Solves a nonconvex problem
      minimize ||x-z||_2^2
      subject to f(x) = x^T P x + q^T x + r ~ 0
      where the relation ~ is given by f.relop (either <= or ==)
490355
def onecons_qcqp(z, f, tol=1e-6):
    

    # if constraint is ineq and z is feasible: z is the solution
    if f.relop == '<=' and f.eval(z) <= 0:
        return z

    if f.eigh is None:
        Psymm = (f.P + f.P.T)/2.
        f.eigh = LA.eigh(np.asarray(Psymm.todense()))

    lmb, Q = f.eigh
    zhat = Q.T.dot(z)
    qhat = Q.T.dot(f.qarray)

    # now solve a transformed problem
    # minimize ||xhat - zhat||_2^2
    # subject to sum(lmb_i xhat_i^2) + qhat^T xhat + r = 0
    # constraint is now equality from
    # complementary slackness
    xhat = lambda nu: -np.divide(nu*qhat-2*zhat, 2*(1+nu*lmb))
    phi = lambda xhat: lmb.dot(np.power(xhat, 2)) + qhat.dot(xhat) + f.r

    s = -np.inf
    e = np.inf
    for l in lmb:
        if l > 0: s = max(s, -1./l)
        if l < 0: e = min(e, -1./l)
    if s == -np.inf:
        s = -1.
        while phi(xhat(s)) <= 0: s *= 2.
    if e == np.inf:
        e = 1.
        while phi(xhat(e)) >= 0: e *= 2.
    while e-s > tol:
        m = (s+e)/2.
        p = phi(xhat(m))
        if p > 0: s = m
        elif p < 0: e = m
        else:
            s = e = m
            break
    nu = (s+e)/2.
    return Q.dot(xhat(nu))
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]
136982
def fd_solve(self):
    
    
    if self.Debug:
      print("qs", self.qs.shape)
      print("Te", self.Te.shape)
      self.calc_max_flexural_wavelength()
      print("maxFlexuralWavelength_ncells', self.maxFlexuralWavelength_ncells")
    
    if self.Solver == "iterative" or self.Solver == "Iterative":
      if self.Debug:
        print("Using generalized minimal residual method for iterative solution")
      if self.Verbose:
        print("Converging to a tolerance of", self.iterative_ConvergenceTolerance, "m between iterations")
      # qs negative so bends down with positive load, bends up with neative load 
      # (i.e. material removed)
      w = isolve.lgmres(self.coeff_matrix, -self.qs, tol=self.iterative_ConvergenceTolerance)  
      self.w = w[0] # Reach into tuple to get my array back
    else:
      if self.Solver == 'direct' or self.Solver == 'Direct':
        if self.Debug:
          print("Using direct solution with UMFpack")
      else:
        print("Solution type not understood:")
        print("Defaulting to direct solution with UMFpack")
      # UMFpack is now the default, but setting true just to be sure in case
      # anything changes
      # qs negative so bends down with positive load, bends up with neative load 
      # (i.e. material removed)
      self.w = spsolve(self.coeff_matrix, -self.qs, use_umfpack=True)
    
    if self.Debug:
      print("w.shape:")
      print(self.w.shape)
      print("w:")
      print(self.w)
767862
def solve(self):
        
        for i in range(4):
            for slot in ["FR", "RB", "BL", "LF"]:
                solver = F2LPairSolver(self.cube, slot)
                if not solver.is_solved():
                    yield tuple([self.cube[slot[i]].colour for i in range(2)]), solver.solve()
                    break
1327284
function(r, x) {
  var r1 = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
  var r2 = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
  curve25519_sqrmodp(r1, x);									// r1 = x^2
  curve25519_mulasmall(r2, x, 486662);							// r2 = Ax
  curve25519_addmodp(r, r1, r2);								//  r = x^2 + Ax
  curve25519_addmodp(r1, r, curve25519_one());					// r1 = x^2 + Ax + 1
  curve25519_mulmodp(r, r1, x);									//  r = x^3 + Ax^2 + x
}
104461
def bukin(self, x):
        """"""
        s = 0
        for k in xrange((1+len(x)) // 2):
            z = x[2 * k]
            y = x[min((2*k + 1, len(x)-1))]
            s += 100 * np.abs(y - 0.01 * z**2)**0.5 + 0.01 * np.abs(z + 10)
        return s
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)
782191
def _NTU_from_P_solver(P1, R1, NTU_min, NTU_max, function, **kwargs):
    '''
    '''
    P1_max = _NTU_from_P_objective(NTU_max, R1, 0, function, **kwargs)
    P1_min = _NTU_from_P_objective(NTU_min, R1, 0, function, **kwargs)
    if P1 > P1_max:
        raise ValueError('No solution possible gives such a high P1; maximum P1=%f at NTU1=%f' %(P1_max, NTU_max))
    if P1 < P1_min:
        raise ValueError('No solution possible gives such a low P1; minimum P1=%f at NTU1=%f' %(P1_min, NTU_min))
    # Construct the function as a lambda expression as solvers don't support kwargs
    to_solve = lambda NTU1: _NTU_from_P_objective(NTU1, R1, P1, function, **kwargs)
    return ridder(to_solve, NTU_min, NTU_max)
614211
def _D_constraint(self, neg_pairs, w):
    
    diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
    diff_sq = diff * diff
    dist = np.sqrt(diff_sq.dot(w))
    sum_deri1 = np.einsum('ij,i', diff_sq, 0.5 / np.maximum(dist, 1e-6))
    sum_deri2 = np.einsum(
        'ij,ik->jk',
        diff_sq,
        diff_sq / (-4 * np.maximum(1e-6, dist**3))[:,None]
    )
    sum_dist = dist.sum()
    return (
      np.log(sum_dist),
      sum_deri1 / sum_dist,
      sum_deri2 / sum_dist - np.outer(sum_deri1, sum_deri1) / (sum_dist * sum_dist)
    )
555844
def  dBinaryRochedz(r, D, q, F):
    
    return -r[2]*(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])**-1.5 -q*r[2]*((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**-1.5
474478
function(f, x0, options) {
    options = options || {}

    this.f = f
    this.n = x0.length
    // cache the init guess
    this.xInitGuess = x0
    this.x = MathUtils.clone(x0)

    this.df = options.df || MathUtils.getApproximateGradientFunction(this.f)
    this.maxIterator = options.maxIterator || 200
    this.err = options.err || 1E-6

    // the gradient of the function evaluated at x[k]: g[k]  (x[0] = x0)
    this.g = this.df(this.x)

    // the inverse of approximate Hessian matrix: B[k]  (B[0] = I)
    this.B = MathUtils.getIdentityMatrix(this.n)

    // direction: p[k]
    this.p = []

    this.stepsize = 0
    this.convergence = Infinity
    this.isConverges = false

    this.iterator = -1
}
1171341
def _nln_constraints(self, nb, nl):
        
        Pmis = NonLinearConstraint("Pmis", nb)
        Qmis = NonLinearConstraint("Qmis", nb)
        Sf = NonLinearConstraint("Sf", nl)
        St = NonLinearConstraint("St", nl)

        return Pmis, Qmis, Sf, St
123141
def constant_cost_withGradients(x):
    
    return np.ones(x.shape[0])[:,None], np.zeros(x.shape)
232905
def opt(self, x_init, f_fp=None, f=None, fp=None):
        
        tnc_rcstrings = ['Local minimum', 'Converged', 'XConverged', 'Maximum number of f evaluations reached',
             'Line search failed', 'Function is constant']

        assert f_fp != None, "TNC requires f_fp"

        opt_dict = {}
        if self.xtol is not None:
            opt_dict['xtol'] = self.xtol
        if self.ftol is not None:
            opt_dict['ftol'] = self.ftol
        if self.gtol is not None:
            opt_dict['pgtol'] = self.gtol

        opt_result = optimize.fmin_tnc(f_fp, x_init, messages=self.messages,
                       maxfun=self.max_f_eval, **opt_dict)
        self.x_opt = opt_result[0]
        self.f_opt = f_fp(self.x_opt)[0]
        self.funct_eval = opt_result[1]
        self.status = tnc_rcstrings[opt_result[2]]
1363597
def forward(Q, p, G, h, A, b, Q_LU, S_LU, R, eps=1e-12, verbose=0, notImprovedLim=3,
            maxIter=20, solver=KKTSolvers.LU_PARTIAL):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    # Find initial values
    if solver == KKTSolvers.LU_FULL:
        D = torch.eye(nineq).repeat(nBatch, 1, 1).type_as(Q)
        x, s, z, y = factor_solve_kkt(
            Q, D, G, A, p,
            torch.zeros(nBatch, nineq).type_as(Q),
            -h, -b if b is not None else None)
    elif solver == KKTSolvers.LU_PARTIAL:
        d = torch.ones(nBatch, nineq).type_as(Q)
        factor_kkt(S_LU, R, d)
        x, s, z, y = solve_kkt(
            Q_LU, d, G, A, S_LU,
            p, torch.zeros(nBatch, nineq).type_as(Q),
            -h, -b if neq > 0 else None)
    elif solver == KKTSolvers.IR_UNOPT:
        D = torch.eye(nineq).repeat(nBatch, 1, 1).type_as(Q)
        x, s, z, y = solve_kkt_ir(
            Q, D, G, A, p,
            torch.zeros(nBatch, nineq).type_as(Q),
            -h, -b if b is not None else None)
    else:
        assert False

    # Make all of the slack variables >= 1.
    M = torch.min(s, 1)[0]
    M = M.view(M.size(0), 1).repeat(1, nineq)
    I = M < 0
    s[I] -= M[I] - 1

    # Make all of the inequality dual variables >= 1.
    M = torch.min(z, 1)[0]
    M = M.view(M.size(0), 1).repeat(1, nineq)
    I = M < 0
    z[I] -= M[I] - 1

    best = {'resids': None, 'x': None, 'z': None, 's': None, 'y': None}
    nNotImproved = 0

    for i in range(maxIter):
        # affine scaling direction
        rx = (torch.bmm(y.unsqueeze(1), A).squeeze(1) if neq > 0 else 0.) + \
            torch.bmm(z.unsqueeze(1), G).squeeze(1) + \
            torch.bmm(x.unsqueeze(1), Q.transpose(1, 2)).squeeze(1) + \
            p
        rs = z
        rz = torch.bmm(x.unsqueeze(1), G.transpose(1, 2)).squeeze(1) + s - h
        ry = torch.bmm(x.unsqueeze(1), A.transpose(
            1, 2)).squeeze(1) - b if neq > 0 else 0.0
        mu = torch.abs((s * z).sum(1).squeeze() / nineq)
        z_resid = torch.norm(rz, 2, 1).squeeze()
        y_resid = torch.norm(ry, 2, 1).squeeze() if neq > 0 else 0
        pri_resid = y_resid + z_resid
        dual_resid = torch.norm(rx, 2, 1).squeeze()
        resids = pri_resid + dual_resid + nineq * mu

        d = z / s
        try:
            factor_kkt(S_LU, R, d)
        except:
            return best['x'], best['y'], best['z'], best['s']

        if verbose == 1:
            print('iter: {}, pri_resid: {:.5e}, dual_resid: {:.5e}, mu: {:.5e}'.format(
                i, pri_resid.mean(), dual_resid.mean(), mu.mean()))
        if best['resids'] is None:
            best['resids'] = resids
            best['x'] = x.clone()
            best['z'] = z.clone()
            best['s'] = s.clone()
            best['y'] = y.clone() if y is not None else None
            nNotImproved = 0
        else:
            I = resids < best['resids']
            if I.sum() > 0:
                nNotImproved = 0
            else:
                nNotImproved += 1
            I_nz = I.repeat(nz, 1).t()
            I_nineq = I.repeat(nineq, 1).t()
            best['resids'][I] = resids[I]
            best['x'][I_nz] = x[I_nz]
            best['z'][I_nineq] = z[I_nineq]
            best['s'][I_nineq] = s[I_nineq]
            if neq > 0:
                I_neq = I.repeat(neq, 1).t()
                best['y'][I_neq] = y[I_neq]
        if nNotImproved == notImprovedLim or best['resids'].max() < eps or mu.min() > 1e32:
            if best['resids'].max() > 1. and verbose >= 0:
                print(INACC_ERR)
            return best['x'], best['y'], best['z'], best['s']

        if solver == KKTSolvers.LU_FULL:
            D = bdiag(d)
            dx_aff, ds_aff, dz_aff, dy_aff = factor_solve_kkt(
                Q, D, G, A, rx, rs, rz, ry)
        elif solver == KKTSolvers.LU_PARTIAL:
            dx_aff, ds_aff, dz_aff, dy_aff = solve_kkt(
                Q_LU, d, G, A, S_LU, rx, rs, rz, ry)
        elif solver == KKTSolvers.IR_UNOPT:
            D = bdiag(d)
            dx_aff, ds_aff, dz_aff, dy_aff = solve_kkt_ir(
                Q, D, G, A, rx, rs, rz, ry)
        else:
            assert False

        # compute centering directions
        alpha = torch.min(torch.min(get_step(z, dz_aff),
                                    get_step(s, ds_aff)),
                          torch.ones(nBatch).type_as(Q))
        alpha_nineq = alpha.repeat(nineq, 1).t()
        t1 = s + alpha_nineq * ds_aff
        t2 = z + alpha_nineq * dz_aff
        t3 = torch.sum(t1 * t2, 1).squeeze()
        t4 = torch.sum(s * z, 1).squeeze()
        sig = (t3 / t4)**3

        rx = torch.zeros(nBatch, nz).type_as(Q)
        rs = ((-mu * sig).repeat(nineq, 1).t() + ds_aff * dz_aff) / s
        rz = torch.zeros(nBatch, nineq).type_as(Q)
        ry = torch.zeros(nBatch, neq).type_as(Q) if neq > 0 else torch.Tensor()

        if solver == KKTSolvers.LU_FULL:
            D = bdiag(d)
            dx_cor, ds_cor, dz_cor, dy_cor = factor_solve_kkt(
                Q, D, G, A, rx, rs, rz, ry)
        elif solver == KKTSolvers.LU_PARTIAL:
            dx_cor, ds_cor, dz_cor, dy_cor = solve_kkt(
                Q_LU, d, G, A, S_LU, rx, rs, rz, ry)
        elif solver == KKTSolvers.IR_UNOPT:
            D = bdiag(d)
            dx_cor, ds_cor, dz_cor, dy_cor = solve_kkt_ir(
                Q, D, G, A, rx, rs, rz, ry)
        else:
            assert False

        dx = dx_aff + dx_cor
        ds = ds_aff + ds_cor
        dz = dz_aff + dz_cor
        dy = dy_aff + dy_cor if neq > 0 else None
        alpha = torch.min(0.999 * torch.min(get_step(z, dz),
                                            get_step(s, ds)),
                          torch.ones(nBatch).type_as(Q))
        alpha_nineq = alpha.repeat(nineq, 1).t()
        alpha_neq = alpha.repeat(neq, 1).t() if neq > 0 else None
        alpha_nz = alpha.repeat(nz, 1).t()

        x += alpha_nz * dx
        s += alpha_nineq * ds
        z += alpha_nineq * dz
        y = y + alpha_neq * dy if neq > 0 else None

    if best['resids'].max() > 1. and verbose >= 0:
        print(INACC_ERR)
    return best['x'], best['y'], best['z'], best['s']
158939
def solve(self, guess):
        
        guess = np.array(guess, dtype=np.float, ndmin=1)
        f = self.func
        args = self._args

        def lmfunc(params, vec):
            vec[:] = f(params, *args).flatten()

        self.lm_prob.set_residual_func(self.data.flatten(),
                                       self.invsigma.flatten(),
                                       lmfunc, None)
        self.lm_soln = soln = self.lm_prob.solve(guess)

        self.params = soln.params
        self.puncerts = soln.perror
        self.covar = soln.covar
        self.mfunc = self.make_frozen_func(soln.params)

        # fvec = resids * invsigma = (data - mdata) * invsigma
        self.resids = soln.fvec.reshape(self.data.shape) / self.invsigma
        self.mdata = self.data - self.resids

        # lm_soln.fnorm can be unreliable ("max(fnorm, fnorm1)" branch)
        self.chisq = (self.lm_soln.fvec**2).sum()
        if soln.ndof > 0:
            self.rchisq = self.chisq / soln.ndof

        return self
1144996
def xors_constraint(v=1):
    
    assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__()
    model, x, y, z = _init()
    r = True
    model.addConsXor([x,y,z], r)
    model.addCons(x==v)
    _optimize("Standard XOR (as boolean)", model)