Solves a presolved standard form LP problem in the form of
min(c) s.t.
A.x = b
lb <= x <= ub
508377
protected int optimizePresolvedStandardLP() throws Exception {
		log.info("optimizePresolvedStandardLP");
		
		long tStart = System.currentTimeMillis();
		
		LPOptimizationRequest lpRequest = getLPOptimizationRequest();
		if(log.isDebugEnabled() && lpRequest.isDumpProblem()){
			log.debug("LP problem: " + lpRequest.toString());
		}
		
		if(this.dim <= -1){
			if (getLb().size() != getUb().size()) {
				log.error("Lower and upper bounds must have the same dimension");
				throw new IllegalArgumentException("Lower and upper bounds must have the same dimension");
			}
			this.dim = getLb().size();
			double minDeltaBoundsValue = Double.MAX_VALUE;
			int minDeltaBoundsIndex = -1;
			for(int i=0; i<getDim(); i++){
				double deltai = getUb().getQuick(i) - getLb().getQuick(i);
				if(deltai < minDeltaBoundsValue){
					minDeltaBoundsValue = deltai;
					minDeltaBoundsIndex = i;
				}
			}
			log.info("min delta bounds index: " + minDeltaBoundsIndex);
			log.info("min delta bounds value: " + minDeltaBoundsValue);
		}
		//this.boundedLb = new boolean[getDim()];
		//this.boundedUb = new boolean[getDim()];
//		for(int i=0; i<getDim(); i++){
//			if(!isLbUnbounded(getLb().getQuick(i))){
//				boundedLb[i] = true;
//				nOfBoundedLb++;
//			}
//			if(!isUbUnbounded(getUb().getQuick(i))){
//				boundedUb[i] = true;
//				nOfBoundedUb++;
//			}
//		}
		this.meq = (this.meq>-1)? this.meq : ((getA()!=null)? getA().rows() : 0);
		//this.mieq = (this.mieq>-1)? this.mieq : (nOfBoundedLb+nOfBoundedUb);
		this.mieq = (this.mieq>-1)? this.mieq : (2 * getDim());
		
		if(log.isDebugEnabled()){
			log.debug("dim : " + getDim());
			log.debug("meq : " + getMeq());
			log.debug("mieq: " + getMieq());
		}
		
		LPOptimizationResponse lpResponse = new LPOptimizationResponse();
		
		DoubleMatrix1D X0 = getInitialPoint();
		if(X0==null){
			DoubleMatrix1D X0NF = getNotFeasibleInitialPoint();
			if(X0NF!=null){
				double rPriX0NFNorm = Math.sqrt(ALG.norm2(rPri(X0NF)));
				DoubleMatrix1D fiX0NF = getFi(X0NF);
				int maxIndex = Utils.getMaxIndex(fiX0NF);
				double maxValue = fiX0NF.get(maxIndex);
				if (log.isDebugEnabled()) {
					log.debug("rPriX0NFNorm :  " + rPriX0NFNorm);
					log.debug("X0NF         :  " + ArrayUtils.toString(X0NF.toArray()));
					log.debug("fiX0NF       :  " + ArrayUtils.toString(fiX0NF.toArray()));
				}
				if(maxValue<0 && rPriX0NFNorm<=getToleranceFeas()){
					//the provided not-feasible starting point is already feasible
					log.debug("the provided initial point is already feasible");
					X0 = X0NF;
				}
			}
			if(X0 == null){
				BasicPhaseILPPDM bf1 = new BasicPhaseILPPDM(this);
				X0 = bf1.findFeasibleInitialPoint();
			}
		}
		
		//check X0 feasibility
		DoubleMatrix1D fiX0 = getFi(X0);
		int maxIndex = Utils.getMaxIndex(fiX0);
		double maxValue = fiX0.get(maxIndex);
		double rPriX0Norm = Math.sqrt(ALG.norm2(rPri(X0)));
		if(maxValue >= 0. || rPriX0Norm > getToleranceFeas()){//must be fi STRICTLY < 0
			log.warn("rPriX0Norm  : " + rPriX0Norm);
			log.warn("ineqX0      : " + ArrayUtils.toString(fiX0.toArray()));
			log.warn("max ineq index: " + maxIndex);
			log.warn("max ineq value: " + maxValue);
			//the point must be INTERNAL, fi are used as denominators
			throw new Exception("initial point must be strictly feasible");
		}

		DoubleMatrix1D V0 = F1.make(getMeq());
		if(getYlb()!=null && getYub()!=null){
			//NB: the Lagrangian multipliers for eq. constraints used in this interior point method (v)
			//are the opposite of the Lagrangian  multipliers for eq. constraints used in the presolver (y)
			//and so Ylb<=y<=Yub becomes -Yub<=v<=-Ylb  
			for(int i=0; i<getMeq(); i++){
				double v0i = 0;
				if(!isLbUnbounded(getYlb().getQuick(i))){
					if(!isUbUnbounded(getYub().getQuick(i))){
						v0i = -(getYub().getQuick(i)+getYlb().getQuick(i))/2;
					}else{
						v0i = -getYlb().getQuick(i);
					}
				}else{
					if(!isUbUnbounded(getYub().getQuick(i))){
						v0i = -getYub().getQuick(i);
					}else{
						v0i = 0;
					}
				}
				V0.setQuick(i, v0i);
			}
		}
		
		DoubleMatrix1D L0 = getInitialLagrangian();
		if(L0!=null){
			for (int j = 0; j < L0.size(); j++) {
				// must be >0
				if(L0.get(j) <= 0){
					throw new IllegalArgumentException("initial lagrangian must be strictly > 0");
				}
			}
		}else{
			L0 = F1.make(getMieq(), 1.);// must be >0 strictly
			if(getZlb()!=null && getZub()!=null){
				//Zlb<= L <=Zub, meaning that:
				//zlb[i] and zub[i] are the bounds on the Lagrangian of the constraint associated with lb[i]<x[i]<ub[i] 
				//note that zlb.size = zub.size = lb.size = ub.size (and = n of variables of the problem (= getDim())
				//and that L.size = nOfBoundedLb + nOfBoundedUb (and in general < 2*getDim())
				int cntLB = 0;
				int cntUB = 0;
				for(int i=0; i<getDim(); i++){
					double zlbi = (isLbUnbounded(getZlb().getQuick(i)))? 0 : getZlb().getQuick(i);//L must be > 0
					double zubi = (isUbUnbounded(getZub().getQuick(i)))? 1 : getZub().getQuick(i);
					L0.setQuick(cntLB, (zubi- zlbi)/2);
					cntLB++;
					L0.setQuick(getDim() + cntUB, (zubi- zlbi)/2);
					cntUB++;
				}
			}else{
				//inequalities comes from the pairs lower bounds-upper bounds
				//in the calculation of the H matrix fro the KKT system, each pairs gives terms of the form:
				//t = tl + tu
				//tl = -L[i] / fi[i] for the lower bound 
				//tu =  L[dim+i] / fi[dim+i] for the upper bound
				//we want t = 1, and hence
				//L[i] > -cc * fi[i]
				//L[dim+i] = (1 + L[i] / fi[i]) * fi[dim+i]
//				double cc = 10;
//				int nOfLB = getMieq()/2;
//				for (int i = 0; i < nOfLB; i++) {
//					L0.setQuick(i, -cc * fiX0.getQuick(i));
//					L0.setQuick(nOfLB + i, (1 - 10)	* fiX0.getQuick(nOfLB + i));
//					double sum = -L0.getQuick(i)/fiX0.getQuick(i)+L0.getQuick(nOfLB + i)/fiX0.getQuick(nOfLB + i);
//					log.debug("sum["+i+"]:  " + sum);
//				}
			}
		}
		if(log.isDebugEnabled()){
			log.debug("X0:  " + ArrayUtils.toString(X0.toArray()));
			log.debug("V0:  " + ArrayUtils.toString(V0.toArray()));
			log.debug("L0:  " + ArrayUtils.toString(L0.toArray()));
		}
		if(log.isInfoEnabled()){
			log.info("toleranceFeas:  " + getToleranceFeas());
			log.info("tolerance    :  " + getTolerance());
		}

		DoubleMatrix1D X = X0;
		DoubleMatrix1D V = V0;
		DoubleMatrix1D L = L0;
		double previousRPriXNorm = Double.NaN;
		double previousRDualXLVNorm = Double.NaN;
		double previousSurrDG = Double.NaN;
		double t;
		int iteration = 0;
		while (true) {
			
			iteration++;
		    // iteration limit condition
			if (iteration == getMaxIteration()+1) {
				lpResponse.setReturnCode(OptimizationResponse.FAILED);
				log.error("Max iterations limit reached");
				throw new Exception("Max iterations limit reached");
			}
			
			//XList.add(XList.size(), X);
			double F0X = getF0(X);
		  if(log.isInfoEnabled()){
				log.info("iteration: " + iteration);
				log.info("f0(X)=" + F0X);
			}
			if(log.isDebugEnabled()){
				log.debug("X=" + ArrayUtils.toString(X.toArray()));
				log.debug("L=" + ArrayUtils.toString(L.toArray()));
				log.debug("V=" + ArrayUtils.toString(V.toArray()));
			}
			
			
			// determine functions evaluations
			DoubleMatrix1D gradF0X = getGradF0(X);
			DoubleMatrix1D fiX     = getFi(X);
			log.debug("fiX=" + ArrayUtils.toString(fiX.toArray()));
			//DoubleMatrix2D GradFiX = getGradFi(X);
			//DoubleMatrix2D GradFiXOLD = getGradFiOLD(X);
			
			// determine t
			double surrDG = getSurrogateDualityGap(fiX, L);
			t = getMu() * getMieq() / surrDG;
			log.debug("t:  " + t);
						
			// determine residuals
			DoubleMatrix1D rPriX    = rPri(X);
			DoubleMatrix1D rCentXLt = rCent(fiX, L, t);
			DoubleMatrix1D rDualXLV = rDual(gradF0X, L, V);
			//DoubleMatrix1D rDualXLVOLD = rDualOLD(GradFiXOLD, gradF0X, L, V);
			//log.debug("delta: " + ALG.normInfinity(rDualXLVOLD.assign(rDualXLV, Functions.minus)));
			double rPriXNorm    = Math.sqrt(ALG.norm2(rPriX));
			double rCentXLtNorm = Math.sqrt(ALG.norm2(rCentXLt));
			double rDualXLVNorm = Math.sqrt(ALG.norm2(rDualXLV));
			double normRXLVt    = Math.sqrt(Math.pow(rPriXNorm, 2) + Math.pow(rCentXLtNorm, 2) + Math.pow(rDualXLVNorm, 2));
			//@TODO: set log.debug not log.info
			log.info("rPri  norm: " + rPriXNorm);
			log.info("rCent norm: " + rCentXLtNorm);
			log.info("rDual norm: " + rDualXLVNorm);
			log.info("surrDG    : " + surrDG);
			
			// custom exit condition
			if(checkCustomExitConditions(X)){
				lpResponse.setReturnCode(OptimizationResponse.SUCCESS);
				break;
			}
			
			// exit condition
			if (rPriXNorm <= getToleranceFeas() && rDualXLVNorm <= getToleranceFeas() && surrDG <= getTolerance()) {
				lpResponse.setReturnCode(OptimizationResponse.SUCCESS);
				break;
			}
			
		  // progress conditions
			if(isCheckProgressConditions()){
				if (!Double.isNaN(previousRPriXNorm)
					&& !Double.isNaN(previousRDualXLVNorm)
					&& !Double.isNaN(previousSurrDG)) {
					if (  (previousRPriXNorm <= rPriXNorm && rPriXNorm >= getToleranceFeas())
						|| (previousRDualXLVNorm <= rDualXLVNorm && rDualXLVNorm >= getToleranceFeas())) {
						log.error("No progress achieved, exit iterations loop without desired accuracy");
						lpResponse.setReturnCode(OptimizationResponse.FAILED);
						throw new Exception("No progress achieved, exit iterations loop without desired accuracy");
					}
				}
				previousRPriXNorm = rPriXNorm;
				previousRDualXLVNorm = rDualXLVNorm;
				previousSurrDG = surrDG;
			}

			// compute primal-dual search direction
			// a) prepare 11.55 system
			DoubleMatrix2D Hpd = GradLSum(L, fiX);
			//DoubleMatrix2D HpdOLD = GradLSumOLD(GradFiXOLD, L, fiX);
			//log.debug("delta: " + ALG.normInfinity(HpdOLD.assign(Hpd, Functions.minus)));
			DoubleMatrix1D gradSum = gradSum(t, fiX);
			DoubleMatrix1D g = null;
			//if(getAT()==null){
			if(getA()==null){
				g = ColtUtils.add(gradF0X, gradSum);
			}else{
				//g = ColtUtils.add(ColtUtils.add(gradF0X, gradSum), ALG.mult(getAT(), V));
				g = ColtUtils.add(ColtUtils.add(gradF0X, gradSum), ColtUtils.zMultTranspose(getA(), V, F1.make(getDim()), 0));
			}
			
			// b) solving 11.55 system
			if(this.kktSolver==null){
				this.kktSolver = new UpperDiagonalHKKTSolver(getDim(), lpRequest.isRescalingDisabled());
				//this.kktSolver = new DiagonalHKKTSolver(getDim(), lpRequest.isRescalingDisabled());
			}
			if(isCheckKKTSolutionAccuracy()){
				kktSolver.setCheckKKTSolutionAccuracy(true);
				kktSolver.setToleranceKKT(getToleranceKKT());
			}
			kktSolver.setHMatrix(Hpd);
			kktSolver.setGVector(g);
			if(getA()!=null){
				kktSolver.setAMatrix(getA());
				//kktSolver.setATMatrix(getAT());
				kktSolver.setHVector(rPriX);
//				if(rPriXNorm > getToleranceFeas()){
//					kktSolver.setHVector(rPriX);
//				}
			}
			DoubleMatrix1D[] sol = kktSolver.solve();
			DoubleMatrix1D stepX = sol[0];
			//double[] signa = new double[stepX.size()];
//			for(int p=0; p<stepX.size(); p++){
//				signa[p] =	Math.signum(stepX.getQuick(p)); 
//			}
			//SList.add(SList.size(), signa);
			DoubleMatrix1D stepV = (sol[1]!=null)? sol[1] : F1.make(0);
			if(log.isDebugEnabled()){
				log.debug("stepX: " + ArrayUtils.toString(stepX.toArray()));
				log.debug("stepV: " + ArrayUtils.toString(stepV.toArray()));
			}

			// c) solving for L
			DoubleMatrix1D stepL = F1.make(getMieq());
			DoubleMatrix1D gradFiStepX = gradFiStepX(stepX);
			for(int i=0; i<getMieq(); i++){
				stepL.setQuick(i, (- L.getQuick(i) * gradFiStepX.getQuick(i) + rCentXLt.getQuick(i)) / fiX.getQuick(i));
			}
			if(log.isDebugEnabled()){
				log.debug("stepL: " + ArrayUtils.toString(stepL.toArray()));
			}

			// line search and update
			// a) sMax computation 
			double sMax = Double.MAX_VALUE;
			for (int j = 0; j < getMieq(); j++) {
				if (stepL.get(j) < 0) {
					sMax = Math.min(-L.get(j) / stepL.get(j), sMax);
				}
			}
			sMax = Math.min(1, sMax);
			double s = 0.99 * sMax;
			// b) backtracking with f
			DoubleMatrix1D X1 = F1.make(X.size());
			DoubleMatrix1D L1 = F1.make(L.size());
			DoubleMatrix1D V1 = F1.make(V.size());
			DoubleMatrix1D fiX1 = null;
			DoubleMatrix1D gradF0X1 = null;
			//DoubleMatrix2D GradFiX1 = null;
			//DoubleMatrix2D GradFiX1 = null;
			DoubleMatrix1D rPriX1 = null;
			DoubleMatrix1D rCentX1L1t = null;
			DoubleMatrix1D rDualX1L1V1 = null;
			int cnt = 0;
			boolean areAllNegative = true;
			while (cnt < 500) {
				cnt++;
				// X1 = X + s*stepX
				X1 = stepX.copy().assign(Mult.mult(s)).assign(X, Functions.plus);
				DoubleMatrix1D ineqValueX1 = getFi(X1);
				areAllNegative = true;
				for (int j = 0; areAllNegative && j < getMieq(); j++) {
					areAllNegative = (Double.compare(ineqValueX1.get(j), 0.) < 0);
				}
				if (areAllNegative) {
					break;
				}
				s = getBeta() * s;
			}
			
			if(!areAllNegative){
				//exited from the feasible region
				throw new Exception("Optimization failed: impossible to remain within the faesible region");
			}
			
			log.debug("s: " + s);
			// c) backtracking with norm
			double previousNormRX1L1V1t = Double.NaN;
			cnt = 0;
			while (cnt < 500) {
				cnt++;
				X1 = ColtUtils.add(X, stepX, s);
				L1 = ColtUtils.add(L, stepL, s);
				V1 = ColtUtils.add(V, stepV, s);
				
				if (isInDomainF0(X1)) {
					fiX1 = getFi(X1);
					gradF0X1 = getGradF0(X1);
					//GradFiX1 = getGradFi(X1);
					
					rPriX1 = rPri(X1);
					rCentX1L1t = rCent(fiX1, L1, t);
					rDualX1L1V1 = rDual(gradF0X1, L1, V1);
					double normRX1L1V1t = Math.sqrt(ALG.norm2(rPriX1)
							                          + ALG.norm2(rCentX1L1t)
							                          + ALG.norm2(rDualX1L1V1));
					//log.debug("normRX1L1V1t: "+normRX1L1V1t);
					if (normRX1L1V1t <= (1 - getAlpha() * s) * normRXLVt) {
						break;
					}
					
					if (!Double.isNaN(previousNormRX1L1V1t)) {
						if (previousNormRX1L1V1t <= normRX1L1V1t) {
							log.warn("No progress achieved in backtracking with norm");
							break;
						}
					}
					previousNormRX1L1V1t = normRX1L1V1t;
				}
				
				s = getBeta() * s;
				//log.debug("s: " + s);
			}

			// update
			X = X1;
			V = V1;
			L = L1;
		}
		
//		if(lpRequest.isCheckOptimalDualityConditions()){
//			//check duality conditions:
//			if(!checkDualityConditions(X, L, V)){
//				log.error("duality conditions not satisfied");
//				lpResponse.setReturnCode(OptimizationResponse.FAILED);
//				throw new Exception("duality conditions not satisfied");
//			}
//		}
		
		if(lpRequest.isCheckOptimalLagrangianBounds()){
			//check equality constraints Lagrangian bounds
//			if(!checkEqConstraintsLagrangianBounds(V)){
//				log.error("equality constraints Lagrangian multipliers bounds not satisfied");
//				lpResponse.setReturnCode(OptimizationResponse.FAILED);
//				throw new Exception("equality constraints Lagrangian multipliers bounds not satisfied");
//			}
			
			//check inequality constraints Lagrangian bounds
//			if(!checkIneqConstraintsLagrangianBounds(X, L)){
//				log.error("inequality constraints Lagrangian multipliers bounds not satisfied");
//				lpResponse.setReturnCode(OptimizationResponse.FAILED);
//				throw new Exception("inequality constraints Lagrangian multipliers bounds not satisfied");
//			}
		}

		long tStop = System.currentTimeMillis();
		log.debug("time: " + (tStop - tStart));
		log.debug("sol : " + ArrayUtils.toString(X.toArray()));
		log.debug("ret code: " + lpResponse.getReturnCode());
//		log.debug("XList : " + ArrayUtils.toString(XList));
//		for(int s=0; s<SList.size(); s++){
//			log.debug("SList : " + ArrayUtils.toString(SList.get(s)));
//		}
		lpResponse.setSolution(X.toArray());
		setLPOptimizationResponse(lpResponse);
		return lpResponse.getReturnCode();
	}
792701
@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        if( B.numRows != numRows )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        X.reshape(numCols,B.numCols);

        U = decomposer.getR(U,true);
        final double gammas[] = decomposer.getGammas();
        final double dataQR[] = QR.data;

        final int BnumCols = B.numCols;

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                a[i] = B.data[i*BnumCols + colB];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                int indexU = n*numRows + n + 1;

                double ub = a[n];
                // U^T*b
                for( int i = n+1; i < numRows; i++ , indexU++ ) {
                    ub += dataQR[indexU]*a[i];
                }

                // gamma*U^T*b
                ub *= gammas[n];

                a[n] -= ub;
                indexU = n*numRows + n + 1;
                for( int i = n+1; i < numRows; i++ , indexU++) {
                    a[i] -= dataQR[indexU]*ub;
                }
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_DDRM.solveU(U.data,a,numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.data[i*X.numCols+colB] = a[i];
            }
        }
    }
682026
def minimise_xyz(xyz):
    """"""
    x, y, z = xyz
    m = max(min(x, y), min(max(x, y), z))
    return (x-m, y-m, z-m)
199118
def solve_kkt(Q_LU, d, G, A, S_LU, rx, rs, rz, ry):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    invQ_rx = rx.btrisolve(*Q_LU)
    if neq > 0:
        h = torch.cat((invQ_rx.unsqueeze(1).bmm(A.transpose(1, 2)).squeeze(1) - ry,
                       invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz), 1)
    else:
        h = invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz

    w = -(h.btrisolve(*S_LU))

    g1 = -rx - w[:, neq:].unsqueeze(1).bmm(G).squeeze(1)
    if neq > 0:
        g1 -= w[:, :neq].unsqueeze(1).bmm(A).squeeze(1)
    g2 = -rs - w[:, neq:]

    dx = g1.btrisolve(*Q_LU)
    ds = g2 / d
    dz = w[:, neq:]
    dy = w[:, :neq] if neq > 0 else None

    return dx, ds, dz, dy
1078217
def solveBinPacking(s,B):
    
    n = len(s)
    U = len(FFD(s,B)) # upper bound of the number of bins
    model = bpp(s,B)
    x,y = model.data

    model.optimize()

    bins = [[] for i in range(U)]
    for (i,j) in x:
        if model.getVal(x[i,j]) > .5:
            bins[j].append(s[i])
    for i in range(bins.count([])):
        bins.remove([])
    for b in bins:
        b.sort()
    bins.sort()
    return bins
1284333
function linComb_(x, y, a, b) {
  var i, c, k, kk;
  k=x.length<y.length ? x.length : y.length;
  kk=x.length;
  for (c=0, i=0; i<k; i++) {
    c+=a*x[i]+b*y[i];
    x[i]=c & mask;
    c>>=bpe;
  }
  for (i=k; i<kk; i++) {
    c+=a*x[i];
    x[i]=c & mask;
    c>>=bpe;
  }
}
190026
private double[][] solveLtransposed(double[][] X) {
    for(int k = L.length - 1; k >= 0; k--) {
      final double[] Lk = L[k], Xk = X[k];
      timesEquals(Xk, 1. / Lk[k]);
      for(int i = 0; i < k; i++) {
        minusTimesEquals(X[i], Xk, Lk[i]);
      }
    }
    return X;
  }
697645
function constrainedMDSGradient(x, fxprime, distances, constraints) {
        var loss = 0, i;
        for (i = 0; i < fxprime.length; ++i) {
            fxprime[i] = 0;
        }

        for (i = 0; i < distances.length; ++i) {
            var xi = x[2 * i], yi = x[2 * i + 1];
            for (var j = i + 1; j < distances.length; ++j) {
                var xj = x[2 * j], yj = x[2 * j + 1],
                    dij = distances[i][j],
                    constraint = constraints[i][j];

                var squaredDistance = (xj - xi) * (xj - xi) + (yj - yi) * (yj - yi),
                    distance$$1 = Math.sqrt(squaredDistance),
                    delta = squaredDistance - dij * dij;

                if (((constraint > 0) && (distance$$1 <= dij)) ||
                    ((constraint < 0) && (distance$$1 >= dij))) {
                    continue;
                }

                loss += 2 * delta * delta;

                fxprime[2*i]     += 4 * delta * (xi - xj);
                fxprime[2*i + 1] += 4 * delta * (yi - yj);

                fxprime[2*j]     += 4 * delta * (xj - xi);
                fxprime[2*j + 1] += 4 * delta * (yj - yi);
            }
        }
        return loss;
    }
600200
def solve(self):
        '''
        
        '''
        self.defUtilityFuncs()
        self.DiscFacEff = self.DiscFac*self.LivPrb
        self.makePFcFunc()
        self.defValueFuncs()
        solution = ConsumerSolution(cFunc=self.cFunc, vFunc=self.vFunc, vPfunc=self.vPfunc,
                                    mNrmMin=self.mNrmMin, hNrm=self.hNrmNow,
                                    MPCmin=self.MPC, MPCmax=self.MPC)
        #solution = self.addSSmNrm(solution)
        return solution
778486
def bpp(s,B):
    
    n = len(s)
    U = len(FFD(s,B)) # upper bound of the number of bins
    model = Model("bpp")
    # setParam("MIPFocus",1)
    x,y = {},{}
    for i in range(n):
        for j in range(U):
            x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
    for j in range(U):
        y[j] = model.addVar(vtype="B", name="y(%s)"%j)

    # assignment constraints
    for i in range(n):
        model.addCons(quicksum(x[i,j] for j in range(U)) == 1, "Assign(%s)"%i)

    # bin capacity constraints
    for j in range(U):
        model.addCons(quicksum(s[i]*x[i,j] for i in range(n)) <= B*y[j], "Capac(%s)"%j)

    # tighten assignment constraints
    for j in range(U):
        for i in range(n):
            model.addCons(x[i,j] <= y[j], "Strong(%s,%s)"%(i,j))

    # tie breaking constraints
    for j in range(U-1):
        model.addCons(y[j] >= y[j+1],"TieBrk(%s)"%j)

    # SOS constraints
    for i in range(n):
        model.addConsSOS1([x[i,j] for j in range(U)])

    model.setObjective(quicksum(y[j] for j in range(U)), "minimize")
    model.data = x,y

    return model
1210682
def lcp(s1, s2):
    '''
    '''
    i = 0
    for i, (c1, c2) in enumerate(zip(s1, s2)):
        if c1 != c2:
            return i
    return min(len(s1), len(s2))
360705
def _minimize_constraints_fun_summation(x):
    '''
    
    '''
    summation = sum([x[i] for i in CONSTRAINT_PARAMS_IDX])
    return CONSTRAINT_UPPERBOUND >= summation >= CONSTRAINT_LOWERBOUND
826744
function(a, b){
    return new BB(
            min(a.l, b.l),
            min(a.b, b.b),
            max(a.r, b.r),
            max(a.t, b.t)
        );
}
95755
def _expon_solve_lam_from_mu(mu, b):
    

    def lam_eq(lam, mu, b):
        # Small offset added to denominator to avoid 0/0 erors
        lam, mu, b = Decimal(lam), Decimal(mu), Decimal(b)
        return ( (1 - (lam*b + 1) * np.exp(-lam*b)) /
                 (lam - lam * np.exp(-lam*b) + Decimal(1e-32)) - mu )

    return optim.brentq(lam_eq, -100, 100, args=(mu, b), disp=True)
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;
    }
1165072
def restricted_cover(l, succsOf):
    
    fzl = frozenset(l)
    lut = dict()
    for i in l:
        lut[i] = fzl.intersection(succsOf(i))
    return lambda x: lut[x]
1200752
def _cho_solve(c_and_lower, b, check_finite=True):
    """"""

    L = c_and_lower[0]
    y = cpxl.solve_triangular(L, b, trans=0, lower=True,
                              check_finite=check_finite)
    return cpxl.solve_triangular(L, y, trans=1, lower=True,
                                 check_finite=check_finite)