Solves the either quadratic or linear program.
172346
def _run_opf(self, HH, CC, AA, ll, uu, xmin, xmax, x0, opt):
        
        N = self._nieq

        if HH.nnz > 0:
            solution = qps_pips(HH, CC, AA, ll, uu, xmin, xmax, x0, opt)
        else:
            solution = qps_pips(None, CC, AA, ll, uu, xmin, xmax, x0, opt)

        return solution
845291
def optimize(self, x0, target):
        """"""
        x = x0
        for _ in range(self.maxiter):
            delta = np.linalg.solve(self.h(x, target), -self.g(x, target))
            x = x + delta
            if np.linalg.norm(delta) < self.tol:
                break
        return x
1307818
def quadraticSolver(a, b, c):
    """"""
    if a == 0:
        if b == 0:  return [] # attempting to solve an equation with infinite (0=0) or impossible (0=3) solutions for x
        else:       return [-c / float(b)]
    else:
        d = b**2.0 - 4*a*c # calculate the discriminant
        if d < 0:   return [] # would require math.sqrt a negative number, yielding a complex number
        solvedPart0 = 1/(2.0*a)
        solvedPart1 = -b * solvedPart0
        if d == 0:  return [ solvedPart1 ] # two values at same 'x'
        else: # i.e. d > 0   apply full quadratic formula
            solvedPart2 = math.sqrt(d)*solvedPart0
            return [ solvedPart1+solvedPart2, solvedPart1-solvedPart2 ]
41533
def solve_fba(self, objective):
        
        self._prob.set_objective(self._v_wt[objective])
        return self._solve(lp.ObjectiveSense.Maximize)
852542
public void solve(double[] b, double[] x) {
        if (b.length != qr.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x 1", qr.nrows(), qr.nrows(), b.length));
        }

        if (x.length != qr.ncols()) {
            throw new IllegalArgumentException("A and x dimensions don't match.");
        }

        if (singular) {
            throw new RuntimeException("Matrix is rank deficient.");
        }

        double[] B = b.clone();
        solve(Matrix.newInstance(B));
        System.arraycopy(B, 0, x, 0, x.length);
    }
634709
public static int ConstraintCompare(Solution solution1, Solution solution2)
   {
      double constraints1 = 0;
      int constraint_count = solution1.getConstraints().size();
      for (int i = 0; i < constraint_count; i++)
      {
         constraints1 += Math.abs(solution1.getConstraint(i));
      }

      double constraints2 = 0;
      for (int i = 0; i < constraint_count; i++)
      {
         constraints2 += Math.abs(solution2.getConstraint(i));
      }

      if ((constraints1 != 0.0) || (constraints2 != 0.0))
      {
         if (constraints1 == 0.0)
         {
            return -1;
         }
         else if (constraints2 == 0.0)
         {
            return 1;
         }
         else
         {
            return Double.compare(constraints1, constraints2);
         }
      }
      else
      {
         return 0;
      }
   }
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
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)
469969
def solve(self, S, dimK=None):
        """"""

        # Use dimK specified in __init__ as default
        if dimK is None and self.dimK is not None:
            dimK = self.dimK

        # Start solve timer
        self.timer.start(['solve', 'solve_wo_eval'])

        # Solve CSC problem on S and do dictionary step
        self.init_vars(S, dimK)
        self.xstep(S, self.lmbda, dimK)
        self.dstep()

        # Stop solve timer
        self.timer.stop('solve_wo_eval')

        # Extract and record iteration stats
        self.manage_itstat()

        # Increment iteration count
        self.j += 1

        # Stop solve timer
        self.timer.stop('solve')

        # Return current dictionary
        return self.getdict()
1355444
def solve_dp(self, lam):
        ''''''
        cur_converge = self.converge+1
        step = 0
        # Get an initial estimate using the GFL
        self.solve_gfl(lam)
        beta2 = np.copy(self.beta)
        while cur_converge > self.converge and step < self.max_dp_steps:
            # Weight each edge differently
            u = lam / (1 + np.abs(self.beta[self.trails[::2]] - self.beta[self.trails[1::2]]))
            # Swap the beta buffers
            temp = self.beta
            self.beta = beta2
            beta2 = temp
            # Solve the edge-weighted GFL problem, which updates beta
            self.solve_gfl(u)
            # Check for convergence
            cur_converge = np.sqrt(((self.beta - beta2)**2).sum())
            step += 1
        self.steps.append(step)
        return self.beta
1233595
def solve_full(z, Fval, DPhival, G, A):    
    M, N=G.shape
    P, N=A.shape

    """"""
    m=M    

    """Primal variable"""
    x=z[0:N]

    """Multiplier for equality constraints"""
    nu=z[N:N+P]

    """Multiplier for inequality constraints"""
    l=z[N+P:N+P+M]

    """Slacks"""
    s=z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = np.diag(l/s)

    """Condensed system"""
    if issparse(DPhival):
        if not issparse(A):
            A = csr_matrix(A)        
        H = DPhival + mydot(G.T, mydot(SIG, G))
        J = bmat([[H, A.T], [A, None]])
    else:
        if issparse(A):
            A = A.toarray()
        J = np.zeros((N+P, N+P))
        J[0:N, 0:N] = DPhival + mydot(G.T, mydot(SIG, G))            
        J[0:N, N:] = A.T
        J[N:, 0:N] = A

    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))

    """Prepare iterative solve via MINRES"""
    sign = np.zeros(N+P)
    sign[0:N/2] = 1.0
    sign[N/2:] = -1.0
    S = diags(sign, 0)
    J_new = mydot(S, csr_matrix(J))
    b_new = mydot(S, b)

    dJ_new = np.abs(J_new.diagonal())
    dPc = np.ones(J_new.shape[0])
    ind = (dJ_new > 0.0)
    dPc[ind] = 1.0/dJ_new[ind]
    Pc = diags(dPc, 0)    
    dxnu, info = minres(J_new, b_new, tol=1e-8, M=Pc)
    
    # dxnu = solve(J, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
861047
def solve(self):  # type: () -> SolverResult
        
        start = time.time()
        root_dependency = Dependency(self._root.name, self._root.version)
        root_dependency.is_root = True

        self._add_incompatibility(
            Incompatibility([Term(root_dependency, False)], RootCause())
        )

        try:
            next = self._root.name
            while next is not None:
                self._propagate(next)
                next = self._choose_package_version()

            return self._result()
        except Exception:
            raise
        finally:
            self._log(
                "Version solving took {:.3f} seconds.\n"
                "Tried {} solutions.".format(
                    time.time() - start, self._solution.attempted_solutions
                )
            )
908829
@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        if( B.numRows != numRows )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numRows);
        X.reshape(numCols,B.numCols);

        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.data[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++ ) {
                double []u = QR[n];

                double vv = u[n];
                u[n] = 1;
                QrHelperFunctions_DDRM.rank1UpdateMultR(a, u, gammas[n], 0, n, numRows, temp.data);
                u[n] = vv;
            }

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

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.data[i*X.numCols+colB] = a.data[i];
            }
        }
    }
1171203
def primal_dual_solve(func, x0, Dfunc, A, b, G, h, args=(), tol=1e-10,
                      maxiter=100, show_progress=True, full_output=False):
    """"""
    fcalls, func = wrap_function(func, args)
    Dfcalls, Dfunc = wrap_function(Dfunc, args)

    M, N = G.shape
    P, N = A.shape

    """Total number of inequality constraints"""
    m = M

    def gap(z):
        r"""Gap-function"""
        l = z[N+P:N+P+M]
        s = z[N+P+M:]
        return mydot(l, s)/m

    def centrality(z):
        r"""Centrality function"""
        l = z[N+P:N+P+M]
        s = z[N+P+M:]
        return np.min(l*s)

    def KKT(z, sigma=0.0):
        r"""KKT system (possible perturbed).Primal variable"""
        x = z[0:N]

        """Multiplier for equality constraints"""
        nu = z[N:N+P]

        """Multiplier for inequality constraints"""
        l = z[N+P:N+P+M]

        """Slacks"""
        s = z[N+P+M:]

        """Evaluate objective function"""
        F = func(x)

        """Dual infeasibility"""
        rdual = F+mydot(A.transpose(), nu)+mydot(G.transpose(), l)

        """Primal infeasibilities"""
        rprim1 = mydot(A, x)-b
        rprim2 = mydot(G, x)-h+s

        """Complementary slackness (perturbed)"""
        mu = gap(z)
        rslack = l*s-sigma*mu

        return np.hstack((rdual, rprim1, rprim2, rslack))

    def step_fast(z, KKTval, LU, G, A, mu, beta, gamma, alpha0):
        r"""Affine scaling step."""
        dz = solve_factorized(z, KKTval, LU, G, A)

        """Reduce step length until slacks s and multipliers l are positive"""
        alpha = 1.0*alpha_0
        k = 0
        for k in range(10):
            z_new = z + alpha*dz
            if np.all( z_new[N+P:] > 0.0 ):
                break
            alpha *= 0.5
            k += 1
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")

        """Reduce step length until iterates lie in correct neighborhood"""
        for k in range(10):
            z_new = z + alpha*dz
            KKTval_new = KKT(z_new)
            dual = mynorm(KKTval_new[0:N])
            prim = mynorm(KKTval_new[N:N+P+M])
            mu_new = gap(z_new)
            cent_new = centrality(z_new)

            if (dual <= beta*mu_new and prim <= beta*mu_new and
                cent_new >= gamma*mu_new):
                break
            alpha *= 0.5
            # alpha *= 0.95
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")
        return z_new, mu_new

    def step_safe(z, KKTval, LU, G, A, mu, beta, gamma):
        r"""Centering step."""
        dz = solve_factorized(z, KKTval, LU, G, A)

        """Reduce step length until slacks s and multipliers l are positive"""
        alpha = 1.0
        k = 0
        for k in range(10):
            z_new = z + alpha*dz
            if np.all( z_new[N+P:] > 0.0 ):
                break
            alpha *= 0.5
            k += 1
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction (pos.) reached")

        """Reduce step length until iterates lie in correct neighborhood
        and mu fulfills Armijo condition"""
        k = 0
        for k in range(10):
            z_new = z+alpha*dz
            KKTval_new = KKT(z_new, sigma=SIGMA)
            dual = mynorm(KKTval_new[0:N])
            prim = mynorm(KKTval_new[N:N+P+M])
            mu_new = gap(z_new)
            cent_new = centrality(z_new)

            if (dual <= beta*mu_new and prim <= beta*mu_new and
                cent_new >= gamma*mu_new and
                mu_new<=(1.0-KAPPA*alpha*(1.0-SIGMA))*mu):
                break
            alpha *= 0.5
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")
        return z_new, mu_new

    """INITIALIZATIONInitial Slacks for inequality constraints"""
    s0 = -1.0*(mydot(G, x0)-h)

    """Initial multipliers for inequality constraints"""
    l0 = 1.0*np.ones(M)

    """Initial multipliers for equality constraints"""
    nu0 = np.zeros(P)

    """Initial point"""
    z0 = np.hstack((x0, nu0, l0, s0))

    """Initial KKT-values"""
    KKTval0 = KKT(z0, sigma=0.0)
    mu0 = gap(z0)
    dual0 = mynorm(KKTval0[0:N])
    prim0 = mynorm(KKTval0[N:N+P+M])

    """Initial neighborhood"""
    beta = BETA * np.sqrt(dual0**2 + prim0**2)/mu0
    gamma = 1.0 * GAMMA_MAX

    """Number of fast steps"""
    t = 0

    """Number of iterations"""
    n = 0

    """Dummy variable for step type"""
    step_type = " "

    if show_progress:
        print("%s %s %s %s %s %s %s" %("iter", "gap", "dual", "primal",
                                       "min", "max", "step"))
    """MAIN LOOP"""
    z = z0
    x = z0[0:N]
    KKTval = KKTval0
    dual = dual0
    prim = prim0
    mu = mu0
    Dfunc_val = Dfunc(x)
    LU = factor(z, Dfunc_val, G, A)

    if full_output:
        info = {'z': []}
        info['z'].append(z)

    for n in range(maxiter):
        if show_progress:
            l=z[N+P:N+P+M]
            s=z[N+P+M:]
            print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                     np.min(l*s), np.max(l*s),
                                                     step_type))
        """Attempt fast step"""
        beta_new = (1.0 + GAMMA_BAR**(t+1)) * beta
        gamma_new = GAMMA_MIN + GAMMA_BAR**(t+1)*(GAMMA_MAX - GAMMA_MIN)
        alpha_0 = 1.0 - np.sqrt(mu)/GAMMA_BAR**t

        if alpha_0 > 0.0:
            z_new, mu_new = step_fast(z, KKTval, LU, G, A, mu,
                                      beta_new, gamma_new, alpha_0)
            if mu_new < RHO * mu:
                """Fast successful"""
                z = z_new
                mu = mu_new
                beta = beta_new
                gamma = gamma_new
                t += 1
                step_type = "f"
            else:
                """Perturbed right-had side"""
                KKTval_pert = 1.0*KKTval
                KKTval_pert[N+P+M:] -= SIGMA * mu
                z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
                step_type = "s"
        else:
            """Perturbed right-hand side"""
            KKTval_pert = 1.0*KKTval
            KKTval_pert[N+P+M:] -= SIGMA * mu
            z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
            step_type = "s"

        """Compute new iterates"""
        KKTval = KKT(z, sigma=0.0)
        dual = mynorm(KKTval[0:N])
        prim = mynorm(KKTval[N:N+P+M])
        x = z[0:N]
        Dfunc_val = Dfunc(x)
        LU = factor(z, Dfunc_val, G, A)
        if full_output:
            info['z'].append(z)
        if mu < tol and dual < tol and prim < tol:
            break
    if n == maxiter - 1:
        raise RuntimeError("Maximum number of iterations reached")

    if show_progress:
        l=z[N+P:N+P+M]
        s=z[N+P+M:]
        print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                 np.min(l*s), np.max(l*s),
                                                 step_type))
    if full_output:
        return z[0:N], info
    else:
        return z[0:N]
1212879
def linsolve(self, A, b):
        
        if self.sparselib == 'umfpack':
            return umfpack.linsolve(A, b)

        elif self.sparselib == 'klu':
            return klu.linsolve(A, b)
352578
def solve(A, b):
    
    """
    A = asarray(A, float)
    b = asarray(b, float)
    if A.shape[0] == 1:

        with errstate(divide="ignore"):
            A_ = array([[1.0 / A[0, 0]]])

        if not isfinite(A_[0, 0]):
            raise LinAlgError("Division error.")

        return dot(A_, b)
    elif A.shape[0] == 2:
        a = A[0, 0]
        b_ = A[0, 1]
        c = A[1, 0]
        d = A[1, 1]
        A_ = array([[d, -b_], [-c, a]])

        with errstate(divide="ignore"):
            A_ /= a * d - b_ * c

        if not npy_all(isfinite(A_)):
            raise LinAlgError("Division error.")

        return dot(A_, b)
    return _solve(A, b)
1039381
def solve(guess_a, guess_b, power, solver='scipy'):
    
    # The problem is 2 dimensional so we need 2 symbols
    x = sp.symbols('x:2', real=True)
    # There is a user specified parameter ``p`` in this problem:
    p = sp.Symbol('p', real=True, negative=False, integer=True)
    # Our system consists of 2-non-linear equations:
    f = [x[0] + (x[0] - x[1])**p/2 - 1,
         (x[1] - x[0])**p/2 + x[1]]
    # We construct our ``SymbolicSys`` instance by passing variables, equations and parameters:
    neqsys = SymbolicSys(x, f, [p])  # (this will derive the Jacobian symbolically)

    # Finally we solve the system using user-specified ``solver`` choice:
    return neqsys.solve([guess_a, guess_b], [power], solver=solver)