Solve nonlinear system using safeguarded Newton iterations


    Parameters
    ----------

    Return
    ------
348999
def newton(f, x, verbose=False, tol=1e-6, maxit=5, jactype='serial'):


    
    if verbose:
        print = lambda txt: old_print(txt)
    else:
        print = lambda txt: None

    it = 0
    error = 10
    converged = False
    maxbacksteps = 30

    x0 = x

    if jactype == 'sparse':
        from scipy.sparse.linalg import spsolve as solve
    elif jactype == 'full':
        from numpy.linalg import solve
    else:
        solve = serial_solve

    while it<maxit and not converged:

        [v,dv] = f(x)

        # TODO: rewrite starting here

#        print("Time to evaluate {}".format(ss-tt)0)

        error_0 = abs(v).max()

        if error_0 < tol:

            if verbose:
                print("> System was solved after iteration {}. Residual={}".format(it,error_0))
            converged = True

        else:

            it += 1

            dx = solve(dv, v)

            # norm_dx = abs(dx).max()

            for bck in range(maxbacksteps):
                xx = x - dx*(2**(-bck))
                vm = f(xx)[0]
                err = abs(vm).max()
                if err < error_0:
                    break

            x = xx

            if verbose:
                print("\t> {} | {} | {}".format(it, err, bck))

    if not converged:
        import warnings
        warnings.warn("Did not converge")
    return [x, it]
262917
def solve_factorized_aug(z, Fval, LU, 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 = diags(l/s, 0)

    """LU is actually the augmented system W"""
    W = LU

    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
    T = diags(sign, 0)        
   
    """Change rhs"""
    b_new = mydot(T, b)

    dW = np.abs(W.diagonal())
    dPc = np.ones(W.shape[0])
    ind = (dW > 0.0)
    dPc[ind] = 1.0/dW[ind]
    Pc = diags(dPc, 0)    
    dxnu, info = minres(W, b_new, tol=1e-10, 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)
    # ds = s*ds
    # SIG = np.diag(l/s)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
1047141
def xstep(self):
        
        """

        self.X = np.asarray(sl.cho_solve_ATAI(
            self.D, self.rho, self.DTS + self.rho * (self.Y - self.U),
            self.lu, self.piv), dtype=self.dtype)

        if self.opt['LinSolveCheck']:
            b = self.DTS + self.rho * (self.Y - self.U)
            ax = self.D.T.dot(self.D.dot(self.X)) + self.rho*self.X
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
721648
def _get_nonlinear_site_term(self, C, vs30, pga_rock):
        
        v_s = np.copy(vs30)
        v_s[vs30 > 760.] = 760.
        # Nonlinear controlling parameter (equation 8)
        f_2 = C["f4"] * (np.exp(C["f5"] * (v_s - 360.)) -
                         np.exp(C["f5"] * 400.))
        fnl = self.CONSTS["f1"] + f_2 * np.log((pga_rock + self.CONSTS["f3"]) /
                                               self.CONSTS["f3"])
        return fnl
24831
def solve(self):
        
        case = self.om.case
        self._base_mva = case.base_mva
        # TODO: Find an explanation for this value.
        self.opt["cost_mult"] = 1e-4

        # Unpack the OPF model.
        self._bs, self._ln, self._gn, _ = self._unpack_model(self.om)
        # Compute problem dimensions.
        self._ipol, _, self._nb, self._nl, _, self._ny, self._nxyz = \
            self._dimension_data(self._bs, self._ln, self._gn)

        # Compute problem dimensions.
        self._ng = len(self._gn)
#        gpol = [g for g in gn if g.pcost_model == POLYNOMIAL]

        # Linear constraints (l <= A*x <= u).
        A, l, u = self.om.linear_constraints()
#        AA, bb = self._linear_constraints(self.om)

        _, xmin, xmax = self._var_bounds()

        # Select an interior initial point for interior point solver.
        x0 = self._initial_interior_point(self._bs, self._gn, xmin, xmax, self._ny)

        # Build admittance matrices.
        self._Ybus, self._Yf, self._Yt = case.Y

        # Optimisation variables.

        self._Pg = self.om.get_var("Pg")
        self._Qg = self.om.get_var("Qg")
        self._Va = self.om.get_var("Va")
        self._Vm = self.om.get_var("Vm")

        # Adds a constraint on the reference bus angles.
#        xmin, xmax = self._ref_bus_angle_constraint(bs, Va, xmin, xmax)

        # Solve using Python Interior Point Solver (PIPS).
        s = self._solve(x0, A, l, u, xmin, xmax)

        Vang, Vmag, Pgen, Qgen = self._update_solution_data(s)

        self._update_case(self._bs, self._ln, self._gn, self._base_mva,
            self._Yf, self._Yt, Vang, Vmag, Pgen, Qgen, s["lmbda"])

        return s
374302
def solve(self, lam):
        ''''''
        if self.penalty == 'dp':
            return self.solve_dp(lam)
        if self.penalty == 'gfl':
            return self.solve_gfl(lam)
        if self.penalty == 'gamlasso':
            return self.solve_gfl(lam)
        raise Exception('Unknown penalty type: {0}'.format(self.penalty))
1165685
private void solveInternalL() {
        // This takes advantage of the diagonal elements always being real numbers

        // solve L*y=b storing y in x
        TriangularSolver_ZDRM.solveL_diagReal(t, vv, n);

        // solve L^T*x=y
        TriangularSolver_ZDRM.solveConjTranL_diagReal(t, vv, n);
    }
637004
private void solve(int method) throws NullPointerException {
    // solution exists
    if(solved) {
      return;
    }

    // bring in reduced row echelon form
    if(!reducedRowEchelonForm) {
      reducedRowEchelonForm(method);
    }

    if(!isSolvable(method)) {
      if(LOG.isDebugging()) {
        LOG.debugFine("Equation system is not solvable!");
      }
      return;
    }

    // compute one special solution
    final int cols = coeff[0].length;
    int numbound = 0, numfree = 0;
    int[] boundIndices = new int[cols], freeIndices = new int[cols];
    x_0 = new double[cols];
    outer: for(int i = 0; i < coeff.length; i++) {
      for(int j = i; j < coeff[row[i]].length; j++) {
        if(coeff[row[i]][col[j]] == 1) {
          x_0[col[i]] = rhs[row[i]];
          boundIndices[numbound++] = col[i];
          continue outer;
        }
      }
      freeIndices[numfree++] = i;
    }

    StringBuilder msg = new StringBuilder();
    if(LOG.isDebugging()) {
      msg.append("\nSpecial solution x_0 = [").append(FormatUtil.format(x_0, ",", FormatUtil.NF4)).append(']') //
          .append("\nbound Indices ").append(FormatUtil.format(boundIndices, ",")) //
          .append("\nfree Indices ").append(FormatUtil.format(freeIndices, ","));
    }

    // compute solution space of homogeneous linear equation system
    Arrays.sort(boundIndices, 0, numbound);
    int freeIndex = 0;
    int boundIndex = 0;
    u = new double[cols][numfree];

    for(int j = 0; j < u[0].length; j++) {
      for(int i = 0; i < u.length; i++) {
        if(freeIndex < numfree && i == freeIndices[freeIndex]) {
          u[i][j] = 1;
        }
        else if(boundIndex < numbound && i == boundIndices[boundIndex]) {
          u[i][j] = -coeff[row[boundIndex]][freeIndices[freeIndex]];
          boundIndex++;
        }
      }
      freeIndex++;
      boundIndex = 0; // Restart
    }

    if(LOG.isDebugging()) {
      msg.append("\nU");
      for(double[] anU : u) {
        msg.append('\n').append(FormatUtil.format(anU, ",", FormatUtil.NF4));
      }
      LOG.debugFine(msg.toString());
    }

    solved = true;
  }
531567
def optimise_newton(x, a, c, tolerance=0.001):
    
    x_new = x
    x_old = x-1 # dummy value
    while np.abs(x_new - x_old).sum() > tolerance:
        x_old = x_new
        x_new = newton_update(x_old, a, c)
    return x_new
453836
def levenberg_marquardt_update(x, a, c, damping=0.001):
    
    hess = hessian(x, a)
    return x - np.linalg.inv(hess + damping*np.diag(hess)).dot(gradient(x, a, c))
95195
def solve(self, solver_klass=None):
        
        # Start the clock.
        t0 = time()

        # Build an OPF model with variables and constraints.
        om = self._construct_opf_model(self.case)
        if om is None:
            return {"converged": False, "output": {"message": "No Ref Bus."}}

        # Call the specific solver.
#        if self.opt["verbose"]:
#            print '\nPYLON Version %s, %s', "0.4.2", "April 2010"
        if solver_klass is not None:
            result = solver_klass(om, opt=self.opt).solve()
        elif self.dc:
#            if self.opt["verbose"]:
#                print ' -- DC Optimal Power Flow\n'
            result = DCOPFSolver(om, opt=self.opt).solve()
        else:
#            if self.opt["verbose"]:
#                print ' -- AC Optimal Power Flow\n'
            result = PIPSSolver(om, opt=self.opt).solve()

        result["elapsed"] = time() - t0

        if self.opt.has_key("verbose"):
            if self.opt["verbose"]:
                logger.info("OPF completed in %.3fs." % result["elapsed"])

        return result
1245326
def solve_limited(self, assumptions=[]):
        

        if self.glucose:
            if self.use_timer:
                 start_time = time.clock()

            # saving default SIGINT handler
            def_sigint_handler = signal.signal(signal.SIGINT, signal.SIG_DFL)

            self.status = pysolvers.glucose41_solve_lim(self.glucose, assumptions)

            # recovering default SIGINT handler
            def_sigint_handler = signal.signal(signal.SIGINT, def_sigint_handler)

            if self.use_timer:
                self.call_time = time.clock() - start_time
                self.accu_time += self.call_time

            return self.status
1059980
def hessian(self, x, y, kappa_ext, ra_0=0, dec_0=0):
        
        gamma1 = 0
        gamma2 = 0
        kappa = kappa_ext
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy
238344
private DoubleMatrix solveEquationSystem(){
		DoubleMatrix R=new DoubleMatrix(partition.getLength());
		DoubleMatrix M=new DoubleMatrix(partition.getLength(),partition.getLength());
		DoubleMatrix partitionAsVector=new DoubleMatrix(partition.getPoints());
		DoubleMatrix shiftedPartition=new DoubleMatrix(partition.getLength());
		for(int j=1; j<shiftedPartition.length;j++) {
			shiftedPartition.put(j, partition.getPoint(j-1));
		}
		DoubleMatrix partitionIncrements= partitionAsVector.sub(shiftedPartition).put(0,1);
		DoubleMatrix kernelValues=new DoubleMatrix(partition.getLength()-1);
		DoubleMatrix M1_1= new DoubleMatrix(1);
		DoubleMatrix MFirstCol= new DoubleMatrix(partition.getLength()-1);
		DoubleMatrix MSubDiagonal= new DoubleMatrix(partition.getLength()-1);
		DoubleMatrix MSubMatrix= new DoubleMatrix(partition.getLength()-1,partition.getLength()-1);
		DoubleMatrix MSubMatrixSum= new DoubleMatrix(partition.getLength()-1);

		for(int i=0;i<independentValues.length;i++){

			DoubleMatrix oneZeroVector= new DoubleMatrix(partition.getLength());
			DoubleMatrix kernelSum= new DoubleMatrix(partition.getLength());
			DoubleMatrix shiftedKernelVector= new DoubleMatrix(partition.getLength());


			for(int r=0;r<partition.getLength()-1;r++){
				oneZeroVector.put(r, 1);
				kernelValues.put( r,kernel.density((partition.getIntervalReferencePoint(r)-independentValues[i])/bandwidth));
				shiftedKernelVector.put(r+1,kernelValues.get( r) );
				kernelSum=kernelSum.add(oneZeroVector.mmul(kernelValues.get(r)));
			}

			R=R.add(shiftedPartition.neg().add(independentValues[i]).mul(shiftedKernelVector)
					.add(partitionIncrements.mul(kernelSum)).mul(dependentValues[i]));

			M1_1=M1_1.add( kernelSum.get(0));

			MFirstCol=MFirstCol.add(
					partitionAsVector.getRange(0,partitionAsVector.length-1).neg().add(independentValues[i])
					.mul(kernelValues).add(
							partitionIncrements.getRange(1,partitionAsVector.length)
							.mul(kernelSum.getRange(1, kernelSum.length))));

			MSubDiagonal=MSubDiagonal.add(
					partitionAsVector.getRange(0,partitionAsVector.length-1).neg().add(independentValues[i])
					.mul(partitionAsVector.getRange(0,partitionAsVector.length-1).neg().add(independentValues[i]))
					.mul(kernelValues).add(
							partitionIncrements.getRange(1,partitionAsVector.length)
							.mul(partitionIncrements.getRange(1,partitionAsVector.length)
									.mul(kernelSum.getRange(1, kernelSum.length)))));

			MSubMatrixSum=MSubMatrixSum.add(
					partitionAsVector.getRange(0, partitionAsVector.length-1).neg().add(independentValues[i])
					.mul(kernelValues).add(
							partitionIncrements.getRange(1, partitionIncrements.length)
							.mul(kernelSum.getRange(1, kernelSum.length))));
		}

		DoubleMatrix partitionIncrementMatrix= new DoubleMatrix(partition.getLength()-1,partition.getLength()-1);
		DoubleMatrix matrixDefine= DoubleMatrix.ones(partition.getLength()-1);
		for(int m=0;m<matrixDefine.length-1;m++) {
			matrixDefine.put(m, 0);
			partitionIncrementMatrix.putColumn(m, matrixDefine.mul(partitionIncrements.get(m+1)));
		}

		MSubMatrix=partitionIncrementMatrix.mulColumnVector(MSubMatrixSum);
		MSubMatrix=MSubMatrix.add(MSubMatrix.transpose()).add(DoubleMatrix.diag(MSubDiagonal));

		int[] rowColIndex =new int[partition.getLength()-1];
		for(int n=0;n<rowColIndex.length;n++) {
			rowColIndex[n]=n+1;
		}

		M.put(0,0,M1_1.get(0));
		M.put(rowColIndex, 0, MFirstCol);
		M.put(0, rowColIndex, MFirstCol.transpose());
		M.put(rowColIndex, rowColIndex, MSubMatrix);

		return Solve.solve(M, R);
	}
1135882
def solve_sdr(prob, *args, **kwargs):
    

    # lifted variables and semidefinite constraint
    X = cvx.Semidef(prob.n + 1)

    W = prob.f0.homogeneous_form()
    rel_obj = cvx.Minimize(cvx.sum_entries(cvx.mul_elemwise(W, X)))
    rel_constr = [X[-1, -1] == 1]

    for f in prob.fs:
        W = f.homogeneous_form()
        lhs = cvx.sum_entries(cvx.mul_elemwise(W, X))
        if f.relop == '==':
            rel_constr.append(lhs == 0)
        else:
            rel_constr.append(lhs <= 0)

    rel_prob = cvx.Problem(rel_obj, rel_constr)
    rel_prob.solve(*args, **kwargs)

    if rel_prob.status not in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE]:
        raise Exception("Relaxation problem status: %s" % rel_prob.status)

    return X.value, rel_prob.value
875006
def newton_call(self):
        
        # system = self.system
        # exec(system.call.newton)

        system = self.system
        dae = self.system.dae

        system.dae.init_fg()
        system.dae.reset_small_g()
        # evaluate algebraic equation mismatches
        for model, pflow, gcall in zip(system.devman.devices,
                                       system.call.pflow, system.call.gcall):
            if pflow and gcall:
                system.__dict__[model].gcall(dae)

        # eval differential equations
        for model, pflow, fcall in zip(system.devman.devices,
                                       system.call.pflow, system.call.fcall):
            if pflow and fcall:
                system.__dict__[model].fcall(dae)

        # reset islanded buses mismatches
        system.Bus.gisland(dae)

        if system.dae.factorize:
            system.dae.init_jac0()
            # evaluate constant Jacobian elements
            for model, pflow, jac0 in zip(system.devman.devices,
                                          system.call.pflow, system.call.jac0):
                if pflow and jac0:
                    system.__dict__[model].jac0(dae)
            dae.temp_to_spmatrix('jac0')

        dae.setup_FxGy()

        # evaluate Gy
        for model, pflow, gycall in zip(system.devman.devices,
                                        system.call.pflow, system.call.gycall):
            if pflow and gycall:
                system.__dict__[model].gycall(dae)

        # evaluate Fx
        for model, pflow, fxcall in zip(system.devman.devices,
                                        system.call.pflow, system.call.fxcall):
            if pflow and fxcall:
                system.__dict__[model].fxcall(dae)

        # reset islanded buses Jacobians
        system.Bus.gyisland(dae)

        dae.temp_to_spmatrix('jac')
706755
def admm_linearized_simple(x, f, g, L, tau, sigma, niter, **kwargs):
    
    callback = kwargs.pop('callback', None)
    z = L.range.zero()
    u = L.range.zero()
    for _ in range(niter):
        x[:] = f.proximal(tau)(x - tau / sigma * L.adjoint(L(x) + u - z))
        z = g.proximal(sigma)(L(x) + u)
        u = L(x) + u - z
        if callback is not None:
            callback(x)