solve_scalarization: scale objective function to find new point
    Parameters:
        - n: number of cities
        - c,t: alternative edge weights, to compute two objective functions
    Returns list of candidate solutions
765481
def solve_scalarization(n,c,t):
    

    model = base_model(n,c,t)   # base model for minimizing cost or time
    x,u,C,T = model.data

    def explore(C1,T1,C2,T2,front):
        """explore: recursively try to find new non-dominated solutions with a scaled objective
        Parameters:
            - C1,T1: cost and time of leftmost point
            - C1,T1: cost and time of rightmost point
            - front: current set of non-dominated solutions
        Returns the updated front
        """
        alpha = float(C1 - C2)/(T2 - T1)
        # print("%s,%s -- %s,%s  (%s)..." % (C1,T1,C2,T2,alpha)
        init = list(front)
        model.setObjective(quicksum((c[i,j] + alpha*t[i,j])*x[i,j] for (i,j) in x), "minimize")
        optimize(model,front)
        front = pareto_front(front)
        # print("... added %s points" % (len(front)-len(init))
        if front == init:
            # print("no points added, returning"
            return front

        CM = model.getVal(C)
        TM = model.getVal(T)
        # print("will explore %s,%s -- %s,%s and %s,%s -- %s,%s" % (C1,T1,CM,TM,CM,TM,C2,T2)
        if TM > T1:
            front = explore(C1,T1,CM,TM,front)
        if T2 > TM:
            front = explore(CM,TM,C2,T2,front)
        return front


    # store the set of solutions for plotting
    cand = []       # to store the set of solutions for plotting
    model.setObjective(T, "minimize")

    stat = optimize(model,cand)
    if stat != "optimal":
        return []
    C1 = model.getVal(C)
    T1 = model.getVal(T)

    # change the objective function to minimize the travel cost
    model.setObjective(C, "minimize")

    stat = optimize(model,cand)
    if stat != "optimal":
        return []
    C2 = model.getVal(C)
    T2 = model.getVal(T)

    front = pareto_front(cand)
    return explore(C1,T1,C2,T2,front)
565991
def solve(self, lam):
        ''''''
        s = weighted_graphtf(self.nnodes, self.y, self.weights, lam,
                             self.Dk.shape[0], self.Dk.shape[1], self.Dk.nnz,
                             self.Dk.row.astype('int32'), self.Dk.col.astype('int32'), self.Dk.data.astype('double'),
                             self.maxsteps, self.converge,
                             self.beta, self.u)
        self.steps.append(s)
        return self.beta
520511
public static LPResult solve(double[] linearObjectiveFunction, List<LPSolver.LPConstraint> linearConstraintsList, boolean nonNegative, boolean maximize) {
        int m = linearConstraintsList.size();

        List<LinearConstraint> constraints = new ArrayList<>(m);
        for(LPSolver.LPConstraint constraint : linearConstraintsList) {
            String sign = constraint.getSign();
            Relationship relationship = null;
            if(LPSolver.GEQ.equals(sign)) {
                relationship = Relationship.GEQ;
            }
            else if(LPSolver.LEQ.equals(sign)) {
                relationship = Relationship.LEQ;
            }
            else if(LPSolver.EQ.equals(sign)) {
                relationship = Relationship.EQ;
            }
            constraints.add(
                    new LinearConstraint(constraint.getContraintBody(),
                    relationship,
                    constraint.getValue())
            );
        }

        SimplexSolver solver = new SimplexSolver();
        PointValuePair solution = solver.optimize(
                new LinearObjectiveFunction(linearObjectiveFunction, 0.0),
                new LinearConstraintSet(constraints),
                maximize?GoalType.MAXIMIZE:GoalType.MINIMIZE,
                new NonNegativeConstraint(nonNegative),
                PivotSelectionRule.BLAND
        );

        LPResult result = new LPResult();
        result.setObjectiveValue(solution.getValue());
        result.setVariableValues(solution.getPoint());

        return result;
    }
333530
def _adapt_params(self, stepsize, stepsize_bar, h_bar, mu, index_i, alpha, n_alpha=1):
        
        gamma = 0.05  # free parameter that controls the amount of shrinkage towards mu
        t0 = 10.0  # free parameter that stabilizes the initial iterations of the algorithm
        kappa = 0.75
        # See equation (6) section 3.2.1 for details

        estimate = 1.0 / (index_i + t0)
        h_bar = (1 - estimate) * h_bar + estimate * (self.delta - alpha / n_alpha)

        stepsize = np.exp(mu - sqrt(index_i) / gamma * h_bar)
        i_kappa = index_i ** (-kappa)
        stepsize_bar = np.exp(i_kappa * np.log(stepsize) + (1 - i_kappa) * np.log(stepsize_bar))

        return stepsize, stepsize_bar, h_bar
205233
def solve_mle_rev(C, tol=1e-10, maxiter=100, show_progress=False, full_output=False,
                  return_statdist=True, **kwargs):
    """"""
    M = C.shape[0]

    """Initial guess for primal-point"""
    z0 = np.zeros(2*M)
    z0[0:M] = 1.0

    """Inequality constraints"""
    # G = np.zeros((M, 2*M))
    # G[np.arange(M), np.arange(M)] = -1.0
    G = -1.0*scipy.sparse.eye(M, n=2*M, k=0)
    h = np.zeros(M)

    """Equality constraints"""
    A = np.zeros((1, 2*M))
    A[0, M] = 1.0
    b = np.array([0.0])

    """Scaling"""
    c0 = C.max()
    C = C/c0

    """Symmetric part"""
    Cs = C + C.T

    """Column sum"""
    c = C.sum(axis=0)

    if scipy.sparse.issparse(C):
        Cs = Cs.tocsr()
        c = c.A1
        A = scipy.sparse.csr_matrix(A)
        F = objective_sparse.F
        DF = objective_sparse.DFsym
        convert_solution = objective_sparse.convert_solution
    else:
        F = objective_dense.F
        DF = objective_dense.DF
        convert_solution = objective_dense.convert_solution

    """PDIP iteration"""
    res = primal_dual_solve(F, z0, DF, A, b, G, h,
                            args=(Cs, c),
                            maxiter=maxiter, tol=tol,
                            show_progress=show_progress,
                            full_output=full_output)
    if full_output:
        z, info = res
    else:
        z = res
    pi, P = convert_solution(z, Cs)

    result = [P]
    if return_statdist:
        result.append(pi)
    if full_output:
        result.append(info)

    return tuple(result) if len(result) > 1 else result[0]
540760
public void computeHessianScaling(DMatrixRMaj scaling ) {

		double max = 0;
		for (int i = 0; i < scaling.numRows; i++) {
			// mathematically it should never be negative but...
			double v = scaling.data[i] = sqrt(abs(scaling.data[i]));
			if( v > max )
				max = v;
		}

		// Add this number to avoid divide by zero.
		// Ceres solver just uses 1 no matter what. That's probably not the best approach because if you're doing
		// with very small numbers scaling will be washed out
		max *= 1e-12;
		for (int i = 0; i < scaling.numRows; i++) {
			scaling.data[i] += max;
		}
	}
84243
@Override
   public Solution[] computeOptimizationProblem(SuitableOptions cloudOffers, QualityInformation requirements,
         Topology topology, int numPlansToGenerate) {

      // To findSolution method, we pass an Empty solution instead of a null
      // value to or create a new method that does not consider the current one.
      // This way may help for replanning, when even the first attempt for
      // solution will be based on the current deployment

      cloudOffers.sortDescendingPerformance();
      Solution[] bestSols = super.findInitialRandomSolutions(cloudOffers, numPlansToGenerate, topology);
      super.setFitnessOfSolutions(bestSols, requirements, topology, cloudOffers);
      if (logHill.isDebugEnabled()) {
         logHill.debug(
               "Start checking the presence of quality attached to solutions after the first generation HILLCLIMB");
         super.checkQualityAttachedToSolutions(bestSols);
      }
      super.sortSolutionsByFitnessAndReplaceNaN(bestSols);

      int numItersNoImprovement = 0;
      while (numItersNoImprovement < super.getMaxIterNoImprove()) { // each
                                                                    // iteration
                                                                    // finds a
                                                                    // peak

         Solution currentSol = super.findRandomSolution(cloudOffers, topology);
         currentSol.setSolutionFitness(super.fitness(currentSol, requirements, topology, cloudOffers));

         boolean neighborsImprove = true;
         int pathLengthTraversed = 0;
         while (neighborsImprove) {
            Solution[] candidates = super.findNeighbors(currentSol, cloudOffers, topology);
            if (candidates == null) {
               logHill.debug("New candidates not found");
               neighborsImprove = false;
               pathLengthTraversed = 0;
            } else {
               logHill.debug(
                     "Found " + candidates.length + " neighbors of the solution: Are " + Arrays.toString(candidates));

               super.setFitnessOfSolutions(candidates, requirements, topology, cloudOffers);
               Solution bestCandidate = super.getSolutionWithMaximumFitness(candidates);
               if (bestCandidate.getSolutionFitness() > currentSol.getSolutionFitness()) {
                  currentSol = bestCandidate;
                  pathLengthTraversed++;

                  if (pathLengthTraversed > 100) {
                     logHill.debug("Something weird is happening, iteration: " + pathLengthTraversed
                           + "without finding local maxima: Fitness: " + currentSol.getSolutionFitness() + " Solution: "
                           + currentSol.toString());
                  }

               } else {// there is not any better neighbor. See if this solution
                       // can be included among the set of best solutions
                  neighborsImprove = false;
                  pathLengthTraversed = 0;
               }
            }
         }

         // We have a local maxima in currentSol, let's see if it is better than
         // the previous maxima.
         if (super.solutionShouldBeIncluded(currentSol, bestSols)) {

            logHill.debug("After " + numItersNoImprovement
                  + " iterations we have found a local maxima solution that is going to be included: "
                  + currentSol.toString());

            super.insertOrdered(bestSols, currentSol);
            numItersNoImprovement = 0;

            if (logHill.isDebugEnabled()) {
               logHill.debug(
                     "Start checking the presence of quality attached to solutions after adding a solution in HILLCLIMB");
               super.checkQualityAttachedToSolutions(bestSols);
            }

         } else {

            logHill.debug("After " + numItersNoImprovement
                  + " iterations we have found a local maxima, but it will not be included: " + currentSol.toString());

            numItersNoImprovement++;
         }

      }

      if (logHill.isDebugEnabled()) {
         logHill.debug(
               "Start checking the presence of quality attached to solutions after adding a solution in HILLCLIMB");
         super.checkQualityAttachedToSolutions(bestSols);
      }

      return bestSols;

   }
602784
def scalar_projection(v1, v2):
    '''
    '''
    return np.dot(v1, v2) / np.linalg.norm(v2)
1372626
private void computeSpan( List<AssociatedPair> points ) {

		Q.reshape(points.size(), 9);
		int index = 0;

		for( int i = 0; i < points.size(); i++ ) {
			AssociatedPair p = points.get(i);

			Point2D_F64 a = p.p2;
			Point2D_F64 b = p.p1;

			// The points are assumed to be in homogeneous coordinates.  This means z = 1
			Q.data[index++] =  a.x*b.x;
			Q.data[index++] =  a.x*b.y;
			Q.data[index++] =  a.x;
			Q.data[index++] =  a.y*b.x;
			Q.data[index++] =  a.y*b.y;
			Q.data[index++] =  a.y;
			Q.data[index++] =      b.x;
			Q.data[index++] =      b.y;
			Q.data[index++] =  1;
		}

		if( !solverNull.process(Q,4,nullspace) )
			throw new RuntimeException("Nullspace solver should never fail, probably bad input");

		// extract the span of solutions for E from the null space
		for( int i = 0; i < 9; i++ ) {
			X[i] = nullspace.unsafe_get(i,0);
			Y[i] = nullspace.unsafe_get(i,1);
			Z[i] = nullspace.unsafe_get(i,2);
			W[i] = nullspace.unsafe_get(i,3);
		}
	}
810744
private List<Double> knnDistances(S solution, L solutionList) {
    List<Double> listOfDistances = new ArrayList<>() ;
    for (int i = 0 ; i< solutionList.size(); i++) {
      double distanceBetweenSolutions = distance.getDistance(solution, solutionList.get(i)) ;
      if (distanceBetweenSolutions != 0) {
        listOfDistances.add(distanceBetweenSolutions) ;
      }
    }

    return listOfDistances ;
  }
161900
def solution_path(self, min_lambda, max_lambda, lambda_bins, verbose=0):
        ''''''
        self.u = np.zeros(self.Dk.shape[0], dtype='double')
        lambda_grid = np.exp(np.linspace(np.log(max_lambda), np.log(min_lambda), lambda_bins))
        aic_trace = np.zeros(lambda_grid.shape) # The AIC score for each lambda value
        aicc_trace = np.zeros(lambda_grid.shape) # The AICc score for each lambda value (correcting for finite sample size)
        bic_trace = np.zeros(lambda_grid.shape) # The BIC score for each lambda value
        dof_trace = np.zeros(lambda_grid.shape) # The degrees of freedom of each final solution
        log_likelihood_trace = np.zeros(lambda_grid.shape)
        beta_trace = []
        best_idx = None
        best_plateaus = None

        if self.edges is None:
            self.edges = defaultdict(list)
            elist = csr_matrix(self.D).indices.reshape((self.D.shape[0], 2))
            for n1, n2 in elist:
                self.edges[n1].append(n2)
                self.edges[n2].append(n1)

        # Solve the series of lambda values with warm starts at each point
        for i, lam in enumerate(lambda_grid):
            if verbose:
                print('#{0} Lambda = {1}'.format(i, lam))

            # Fit to the final values
            beta = self.solve(lam)

            if verbose:
                print('Calculating degrees of freedom')

            # Count the number of free parameters in the grid (dof) -- TODO: the graph trend filtering paper seems to imply we shouldn't multiply by (k+1)?
            dof_vals = self.Dk_minus_one.dot(beta) if self.k > 0 else beta
            plateaus = calc_plateaus(dof_vals, self.edges, rel_tol=0.01) if (self.k % 2) == 0 else nearly_unique(dof_vals, rel_tol=0.03)
            dof_trace[i] = max(1,len(plateaus)) #* (k+1)

            if verbose:
                print('Calculating Information Criteria')

            # Get the negative log-likelihood
            log_likelihood_trace[i] = -0.5 * ((self.y - beta)**2).sum()

            # Calculate AIC = 2k - 2ln(L)
            aic_trace[i] = 2. * dof_trace[i] - 2. * log_likelihood_trace[i]
            
            # Calculate AICc = AIC + 2k * (k+1) / (n - k - 1)
            aicc_trace[i] = aic_trace[i] + 2 * dof_trace[i] * (dof_trace[i]+1) / (len(beta) - dof_trace[i] - 1.)

            # Calculate BIC = -2ln(L) + k * (ln(n) - ln(2pi))
            bic_trace[i] = -2 * log_likelihood_trace[i] + dof_trace[i] * (np.log(len(beta)) - np.log(2 * np.pi))

            # Track the best model thus far
            if best_idx is None or bic_trace[i] < bic_trace[best_idx]:
                best_idx = i
                best_plateaus = plateaus

            # Save the trace of all the resulting parameters
            beta_trace.append(np.array(beta))
            
            if verbose:
                print('DoF: {0} AIC: {1} AICc: {2} BIC: {3}\n'.format(dof_trace[i], aic_trace[i], aicc_trace[i], bic_trace[i]))

        if verbose:
            print('Best setting (by BIC): lambda={0} [DoF: {1}, AIC: {2}, AICc: {3} BIC: {4}]'.format(lambda_grid[best_idx], dof_trace[best_idx], aic_trace[best_idx], aicc_trace[best_idx], bic_trace[best_idx]))

        return {'aic': aic_trace,
                'aicc': aicc_trace,
                'bic': bic_trace,
                'dof': dof_trace,
                'loglikelihood': log_likelihood_trace,
                'beta': np.array(beta_trace),
                'lambda': lambda_grid,
                'best_idx': best_idx,
                'best': beta_trace[best_idx],
                'plateaus': best_plateaus}
968470
public <S extends Solution<?>> void execute(List<S> solutionsList) {
    switch (config.scalarizationType) {
      case ANGLE_UTILITY:
        if (config.extremePoints == null)
          ScalarizationUtils.angleUtility(solutionsList);
        else
          ScalarizationUtils.angleUtility(solutionsList, config.extremePoints);
        break;
      case CHEBYSHEV:
        if (config.idealValues == null)
          ScalarizationUtils.chebyshev(solutionsList);
        else
          ScalarizationUtils.chebyshev(solutionsList, config.idealValues);
        break;
      case NASH:
        if (config.nadirValues == null)
          ScalarizationUtils.nash(solutionsList);
        else
          ScalarizationUtils.nash(solutionsList, config.nadirValues);
        break;
      case SUM_OF_OBJECTIVES:
        ScalarizationUtils.sumOfObjectives(solutionsList);
        break;
      case PRODUCT_OF_OBJECTIVES:
        ScalarizationUtils.productOfObjectives(solutionsList);
        break;
      case TRADEOFF_UTILITY:
        ScalarizationUtils.tradeoffUtility(solutionsList);
        break;
      case UNIFORM:
        ScalarizationUtils.uniform(solutionsList);
        break;
      case WEIGHTED_CHEBYSHEV:
        if (config.idealValues == null)
          ScalarizationUtils.weightedChebyshev(solutionsList, config.weights);
        else
          ScalarizationUtils.weightedChebyshev(solutionsList, config.idealValues, config.weights);
        break;
      case WEIGHTED_PRODUCT:
        ScalarizationUtils.weightedProduct(solutionsList, config.weights);
        break;
      case WEIGHTED_SUM:
        ScalarizationUtils.weightedSum(solutionsList, config.weights);
        break;
    }
  }
602089
def get_scaling_function(y_m, y_M, x_m, x_M):
    
    if x_M <= x_m:
        return lambda _: y_m
    k = (y_M - y_m) / (x_M - x_m)
    b = y_m - k * x_m
    return lambda _: int(k * _ + b)
1002643
def tsp(edges, constraint_generation=False):
    
    edges = populate_edge_weights(edges)
    incoming, outgoing, nodes = node_to_edge(edges)

    nedges, nnodes = len(edges), len(nodes)
    L = LPInstance()

    L.add_objective(edges, objective=MINIMIZE)
    balance = []
    # For each node, select exactly 1 incoming and 1 outgoing edge
    for v in nodes:
        incoming_edges = incoming[v]
        outgoing_edges = outgoing[v]
        icc = summation(incoming_edges)
        occ = summation(outgoing_edges)
        balance.append("{0} = 1".format(icc))
        balance.append("{0} = 1".format(occ))

    # Subtour elimination - Miller-Tucker-Zemlin (MTZ) formulation
    # <http://en.wikipedia.org/wiki/Travelling_salesman_problem>
    # Desrochers and laporte, 1991 (DFJ) has a stronger constraint
    # See also:
    # G. Laporte / The traveling salesman problem: Overview of algorithms
    start_step = nedges + 1
    u0 = nodes[0]
    nodes_to_steps = dict((n, start_step + i) for i, n in enumerate(nodes[1:]))
    edge_store = dict((e[:2], i) for i, e in enumerate(edges))
    mtz = []
    for i, e in enumerate(edges):
        a, b = e[:2]
        if u0 in (a, b):
            continue
        na, nb = nodes_to_steps[a], nodes_to_steps[b]
        con_ab = " x{0} - x{1} + {2}x{3}".format(na, nb, nnodes - 1, i + 1)
        if (b, a) in edge_store:  # This extra term is the stronger DFJ formulation
            j = edge_store[(b, a)]
            con_ab += " + {0}x{1}".format(nnodes - 3, j + 1)
        con_ab += " <= {0}".format(nnodes - 2)
        mtz.append(con_ab)

    # Step variables u_i bound between 1 and n, as additional variables
    bounds = []
    for i in xrange(start_step, nedges + nnodes):
        bounds.append(" 1 <= x{0} <= {1}".format(i, nnodes - 1))

    L.add_vars(nedges)

    """
    Constraint generation seek to find 'cuts' in the LP problem, by solving the
    relaxed form. The subtours were then incrementally added to the constraints.
    """
    if constraint_generation:
        L.constraints = balance
        subtours = []
        while True:
            selected, obj_val = L.lpsolve()
            results = sorted(x for i, x in enumerate(edges) if i in selected) \
                if selected else None
            if not results:
                break
            G = edges_to_graph(results)
            cycles = list(nx.simple_cycles(G))
            if len(cycles) == 1:
                break
            for c in cycles:
                incident = [edge_store[a, b] for a, b in pairwise(c + [c[0]])]
                icc = summation(incident)
                subtours.append("{0} <= {1}".format(icc, len(incident) - 1))
            L.constraints = balance + subtours
    else:
        L.constraints = balance + mtz
        L.add_vars(nnodes - 1, offset=start_step, binary=False)
        L.bounds = bounds
        selected, obj_val = L.lpsolve()
        results = sorted(x for i, x in enumerate(edges) if i in selected) \
            if selected else None

    return results
537272
def main(dimension, iterations):
    
    optimizer = PSOOptimizer()
    solution = optimizer.minimize(sphere, -5.12, 5.12, dimension,
                                  max_iterations(iterations))
    return solution, optimizer
358906
def _scalePoints(points, scale=1, convertToInteger=True):
    
    if convertToInteger:
        points = [
            (int(round(x * scale)), int(round(y * scale)))
            for (x, y) in points
        ]
    else:
        points = [(x * scale, y * scale) for (x, y) in points]
    return points
555727
def scale(self, scalar):
        ''''''
        return self.__class__([self.coefficients[i] * scalar for i in _range(len(self))])