Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Obtain initial point from Mehrotra's heuristic.
(x, y, z) = self.set_initial_guess(**kwargs)
# Slack variables are the trailing variables in x.
s = x[on:]
ns = self.nSlacks
# Initialize steps in dual variables.
dz = np.zeros(ns)
# Allocate room for right-hand side of linear systems.
rhs = self.initialize_rhs()
finished = False
iter = 0
setup_time = cputime()
# Main loop.
while not finished:
# Display initial header every so often.
if iter % 50 == 0:
self.log.info(self.header)
self.log.info('-' * len(self.header))
# Compute residuals.
pFeas = A * x - b
comp = s * z
sz = sum(comp) # comp = Sz
Qx = Q * x[:on]
dFeas = y * A
dFeas[:on] -= self.c + Qx # dFeas1 = A1'y - c - Qx
# Collect basic info about the problem.
zero = np.zeros(n)
self.b = -qp.cons(zero) # Right-hand side
self.c0 = qp.obj(zero) # Constant term in objective
self.c = qp.grad(zero[:on]) # Cost vector
self.Q = PysparseMatrix(matrix=qp.hess(zero[:on],
np.zeros(qp.original_m)))
# Apply in-place problem scaling if requested.
self.prob_scaled = False
if scale:
self.t_scale = cputime()
self.scale()
self.t_scale = cputime() - self.t_scale
else:
# self.scale() sets self.normQ to the Frobenius norm of Q
# and self.normA to the Frobenius norm of A as a by-product.
# If we're not scaling, set normQ and normA manually.
self.normQ = self.Q.matrix.norm('fro')
self.normA = self.A.matrix.norm('fro')
self.normb = norm_infty(self.b)
self.normc = norm_infty(self.c)
self.normbc = 1 + max(self.normb, self.normc)
# Initialize augmented matrix.
self.H = self.initialize_kkt_matrix()
# It will be more efficient to keep the diagonal of Q around.
self.diagQ = self.Q.take(range(qp.original_n))
def solve(self):
"""Solve."""
n = self.n
x_norm2 = 0.0 # Squared norm of current iterate x, not counting x_feas
# Obtain initial projected residual
self.t_solve = cputime()
if self.qp.A is not None:
if self.factorize and not self.factorized:
self.perform_factorization()
if self.b is not None:
self.rhs[:n] = self.qp.grad(self.x_feasible)
self.rhs[n:] = 0.0
else:
self.rhs[:n] = self.qp.c
self.proj.solve(self.rhs)
r = g = self.proj.x[:n]
self.v = self.proj.x[n:]
# self.CheckAccurate()
else:
g = self.qp.c
# Print out header, say, every 20 iterations.
if self.iter % 20 == 0:
self.log.info(self.header)
self.log.info(self.format, self.iter, self.f, self.pgnorm,
cons_norm_new, al_model.penalty, bc_solver.iter,
bc_solver.status, self.omega, self.eta)
try:
self.post_iteration()
except UserExitRequest:
self.status = "usr"
exitIter = self.niter_total > self.max_iter
self.tsolve = cputime() - tick # Solve time
if slack_model.m != 0:
self.pi_norm = norm(al_model.pi)
self.cons_norm = norm(slack_model.cons(self.x))
else:
self.pi_norm = None
self.cons_norm = None
# Solution output, etc.
if exitOptimal:
self.status = "opt"
self.log.debug('Optimal solution found \n')
elif not exitOptimal and self.status is None:
self.status = "iter"
self.log.debug('Maximum number of iterations reached \n')
self.omega_rel * self.pg0 + self.omega_abs)
self.eta_opt = kwargs.get("eta_opt",
self.eta_rel * cons_norm + self.eta_abs)
self.log.info(u"ηₒₚₜ = %8.1e, ωₒₚₜ = %8.1e",
self.eta_opt, self.omega_opt)
self.iter = 0
self.inner_fail_count = 0
self.niter_total = 0
infeas_iter = 0
exitIter = False
exitOptimal = self.check_convergence(PdL_norm, cons_norm)
tick = cputime()
# Print out header and initial log.
if self.iter % 20 == 0:
self.log.info(self.header)
self.log.info(self.format0, self.iter, self.f, self.pg0,
self.cons0, al_model.penalty, "", "", self.omega,
self.eta)
while not (exitOptimal or exitIter):
self.iter += 1
# Perform bound-constrained minimization
bc_solver = self.setup_bc_solver()
bc_solver.solve()
self.x = bc_solver.x.copy() # may not be useful.
self.niter_total += bc_solver.iter + 1
def find_feasible(self):
"""Obtain `x_feasible` satisfying the constraints.
`rhs` must have been specified.
"""
n = self.n
self.log.debug('Obtaining feasible solution...')
self.t_feasible = cputime()
self.rhs[n:] = self.b
self.proj.solve(self.rhs)
self.x_feasible = self.proj.x[:n].copy()
self.t_feasible = cputime() - self.t_feasible
self.check_accurate()
self.log.debug('... done (%-5.2fs)' % self.t_feasible)
return
def find_feasible(self):
"""Obtain `x_feasible` satisfying the constraints.
`rhs` must have been specified.
"""
n = self.n
self.log.debug('Obtaining feasible solution...')
self.t_feasible = cputime()
self.rhs[n:] = self.b
self.proj.solve(self.rhs)
self.x_feasible = self.proj.x[:n].copy()
self.t_feasible = cputime() - self.t_feasible
self.check_accurate()
self.log.debug('... done (%-5.2fs)' % self.t_feasible)
return
# set stopping tolerance
tol = self.reltol * Fnorm0 + self.abstol
self.tsolve = 0
self.optimal = optimal = Fnorm <= tol
if optimal:
status = 'Optimal solution found'
short_status = 'opt'
tired = self.itn > self.itermax
finished = optimal or tired
self.itn = 0
tick = cputime()
# Main loop.
while not finished:
self.x_old = x.copy()
self.gL_old = gL.copy()
# update penalty parameter
self.merit.penalty = 1.0 / self.new_penalty(Fnorm)
# compute extrapolation step
self.assemble_linear_system(x, y)
rhs = self.assemble_rhs(g, J, y, c)
status, short_status, solved, dx, dy = self.solve_linear_system(
rhs, J=J)
else:
r = range(self.n)
P.put(1, r, r)
# for i in range(self.n):
# P[i,i] = 1
P[self.n:, :self.n] = self.A
# Add regularization if requested.
if self.dreg > 0.0:
r = range(self.n, self.n + self.m)
P.put(-self.dreg, r, r)
msg = 'Factorizing projection matrix '
msg += '(size %-d, nnz = %-d)...' % (P.shape[0], P.nnz)
self.log.debug(msg)
self.t_fact = cputime()
self.proj = LBLContext(P)
self.t_fact = cputime() - self.t_fact
self.log.debug('... done (%-5.2fs)' % self.t_fact)
self.factorized = True
return
P.put(1, r, r)
# for i in range(self.n):
# P[i,i] = 1
P[self.n:, :self.n] = self.A
# Add regularization if requested.
if self.dreg > 0.0:
r = range(self.n, self.n + self.m)
P.put(-self.dreg, r, r)
msg = 'Factorizing projection matrix '
msg += '(size %-d, nnz = %-d)...' % (P.shape[0], P.nnz)
self.log.debug(msg)
self.t_fact = cputime()
self.proj = LBLContext(P)
self.t_fact = cputime() - self.t_fact
self.log.debug('... done (%-5.2fs)' % self.t_fact)
self.factorized = True
return