class QPModel(NLPModel):
    u"""Generic class to represent a quadratic programming (QP) problem.

    minimize    cᵀx + 1/2 xᵀHx
    subject to  L ≤ A x ≤ U
                l ≤ x ≤ u.
    """

    def __init__(self, c, H, A=None, name='GenericQP', **kwargs):
        """Initialize a QP with linear term `c`, Hessian `H` and Jacobian `A`.

        :parameters:
            :c:   Numpy array to represent the linear objective
            :A:   linear operator to represent the constraint matrix.
                  It must be possible to perform the operations `A*x`
                  and `A.T*y` for Numpy arrays `x` and `y` of appropriate size.
                  If `A` is `None`, it will be replaced with an empty linear
                  operator.
            :H:   linear operator to represent the objective Hessian.
                  It must be possible to perform the operation `H*x`
                  for a Numpy array `x` of appropriate size. The operator `H`
                  should be symmetric.

        See the documentation of `NLPModel` for futher information.
        """
        # Basic checks.
        n = c.shape[0]
        if A is None:
            m = 0
            self.A = LinearOperator(n, 0,
                                    lambda x: np.empty((0, 1)),
                                    matvec_transp=lambda y: np.empty((n, 0)),
                                    dtype=np.float)
        else:
            if A.shape[1] != n or H.shape[0] != n or H.shape[1] != n:
                raise ValueError('Shapes are inconsistent')
            m = A.shape[0]
            self.A = A

        super(QPModel, self).__init__(n=n, m=m, name=name, **kwargs)
        self.c = c
        self.H = H

        # Default classification of constraints
        self._lin = range(self.m)             # Linear    constraints
        self._nln = []                        # Nonlinear constraints
        self._net = []                        # Network   constraints
        self._nlin = len(self.lin)            # Number of linear constraints
        self._nnln = len(self.nln)            # Number of nonlinear constraints
        self._nnet = len(self.net)            # Number of network constraints

    def obj(self, x):
        """Evaluate the objective function at x."""
        cHx = self.hprod(x, 0, x)
        cHx *= 0.5
        cHx += self.c
        return np.dot(cHx, x)

    def grad(self, x):
        """Evaluate the objective gradient at x."""
        Hx = self.hprod(x, 0, x)
        Hx += self.c
        return Hx

    def cons(self, x):
        """Evaluate the constraints at x."""
        if isinstance(self.A, np.ndarray):
            return np.dot(self.A, x)
        return self.A * x

    def A(self, x):
        """Evaluate the constraints Jacobian at x."""
        return self.A

    def jac(self, x):
        """Evaluate the constraints Jacobian at x."""
        return self.A

    def jprod(self, x, p):
        """Evaluate Jacobian-vector product at x with p."""
        return self.cons(p)

    def jtprod(self, x, p):
        """Evaluate transposed-Jacobian-vector product at x with p."""
        if isinstance(self.A, np.ndarray):
            return np.dot(self.A.T, p)
        return self.A.T * p

    def hess(self, x, z):
        """Evaluate Lagrangian Hessian at (x, z)."""
        return self.H

    def hprod(self, x, z, p):
        """Hessian-vector product.

        Evaluate matrix-vector product between the Hessian of the Lagrangian at
        (x, z) and p.
        """
        if isinstance(self.H, np.ndarray):
            return np.dot(self.H, p)
        return self.H * p class C1LineModel(NLPModel):
    u"""Restriction of a C¹ model to a line.

    An instance of this class is a model representing the original model
    restricted to the line x + td. More precisely, if the original objective
    is f: ℝⁿ → ℝ, then, given x ∈ ℝⁿ and d ∈ ℝⁿ (d≠0) a fixed direction, the
    objective of the restricted model is ϕ: ℝ → ℝ defined by

        ϕ(t) := f(x + td).

    Similarly, if the original constraints are c: ℝⁿ → ℝᵐ, the constraints of
    the restricted model are γ: ℝ → ℝᵐ defined by

        γ(t) := c(x + td).

    The functions f and c are only assumed to be C¹, i.e., only values and
    first derivatives of ϕ and γ are defined.
    """

    def __init__(self, model, x, d, **kwargs):
        """Instantiate the restriction of a model to the line x + td.

        :parameters:
            :model: ```NLPModel``` whose objective is to be restricted
            :x: Numpy array
            :d: Numpy array assumed to be nonzero (no check is performed).
        """
        name = "line-" + model.name
        kwargs.pop("name", None)
        kwargs.pop("Lvar", None)
        kwargs.pop("Uvar", None)
        kwargs.pop("Lcon", None)
        kwargs.pop("Ucon", None)
        pos = where(d > 0)
        neg = where(d < 0)
        tmax = Min((model.Uvar[pos] - x[pos]) / d[pos])
        tmax = min(tmax, Min((model.Lvar[neg] - x[neg]) / d[neg]))
        tmin = Max((model.Lvar[pos] - x[pos]) / d[pos])
        tmin = max(tmin, Max((model.Uvar[neg] - x[neg]) / d[neg]))
        super(C1LineModel, self).__init__(1,
                                          m=model.ncon,
                                          name=name,
                                          x0=kwargs.get("x0", 0.0),
                                          Lvar=np.array([tmin]),
                                          Uvar=np.array([tmax]),
                                          Lcon=model.Lcon,
                                          Ucon=model.Ucon,
                                          **kwargs)
        if tmin > tmax:
            # no intersection between the original bounds and the line
            self.logger.warn("restricted model is infeasible")
        self.__x = x
        self.__d = d
        self.__f = None  # most recent objective value of `model`
        self.__g = None  # most recent objective gradient of `model`
        self.__c = None  # most recent constraint values of `model`
        self.__model = model

    @property
    def x(self):
        return self.__x

    @property
    def d(self):
        return self.__d

    @property
    def dir(self):
        return self.__d

    @property
    def objval(self):
        return self.__f

    @property
    def gradval(self):
        return self.__g

    @property
    def conval(self):
        return self.__c

    @property
    def model(self):
        return self.__model

    def obj(self, t, x=None):
        u"""Evaluate ϕ(t) = f(x + td).

        :keywords:
            :x: full-space x+td if that vector has already been formed.
        """
        xtd = (self.x + t * self.d) if x is None else x
        self.__f = self.model.obj(xtd)
        return self.objval

    def grad(self, t, x=None):
        u"""Evaluate ϕ'(t) = ∇f(x + td)ᵀ d.

        :keywords:
            :x: full-space x+td if that vector has already been formed.
        """
        xtd = (self.x + t * self.d) if x is None else x
        self.__g = self.model.grad(xtd)
        return np.dot(self.gradval, self.d)

    def cons(self, t, x=None):
        u"""Evaluate γ(t) = c(x + td).

        :keywords:
            :x: full-space x+td if that vector has already been formed.
        """
        xtd = (self.x + t * self.d) if x is None else x
        self.__c = self.model.cons(xtd)
        return self.conval

    def jac(self, t, x=None):
        u"""Evaluate γ'(t) = J(x + td) d.

        :keywords:
            :x: full-space x+td if that vector has already been formed.
        """
        xtd = (self.x + t * self.d) if x is None else x
        return self.model.jprod(xtd, self.d)

    def jprod(self, t, v, x=None):
        u"""Jacobian-vector product, which is just γ'(t)*v with v ∈ ℝ.

        :keywords:
            :x: full-space x+td if that vector has already been formed.
        """
        return self.jac(t, x=x) * v

    def jtprod(self, t, u, x=None):
        u"""Transposed-Jacobian-vector product γ'(t)ᵀ u with u ∈ ℝᵐ.

        :keywords:
            :x: full-space x+td if that vector has already been formed.
        """
        return np.dot(self.jac(t, x=x), u) class ArmijoLineSearch(LineSearch):
    """Armijo backtracking linesearch."""

    def __init__(self, *args, **kwargs):
        u"""Instantiate an Armijo backtracking linesearch.

        The search stops as soon as a step size t is found such that

            ϕ(t) <= ϕ(0) + t * ftol * ϕ'(0)

        where 0 < ftol < 1 and ϕ'(0) is the directional derivative of
        a merit function f in the descent direction d. true.

        :keywords:
            :ftol: constant used in Armijo condition (default: 1.0e-4)
            :bkmax: maximum number of backtracking steps (default: 20)
            :decr: factor by which to reduce the steplength
                   during the backtracking (default: 1.5).
        """
        name = kwargs.pop("name", "Armijo linesearch")
        super(ArmijoLineSearch, self).__init__(*args, name=name, **kwargs)
        self.__ftol = max(min(kwargs.get("ftol", 1.0e-4), 1 - sqeps), sqeps)
        self.__bkmax = max(kwargs.get("bkmax", 20), 0)
        self.__decr = max(min(kwargs.get("decr", 1.5), 100), 1.001)
        self._bk = 0
        return

    @property
    def ftol(self):
        return self.__ftol

    @property
    def bkmax(self):
        return self.__bkmax

    @property
    def decr(self):
        return self.__decr

    @property
    def bk(self):
        return self._bk

    def next(self):
        if self.trial_value <= self.value + self.step * self.ftol * self.slope:
            raise StopIteration()

        self._bk += 1
        if self._bk > self.__bkmax:
            raise LineSearchFailure("backtracking limit exceeded")

        self._step /= self.decr
        if self.step < self.stepmin:
            raise LineSearchFailure("linesearch step too small")

        self._trial_iterate = self.linemodel.x + self.step * self.linemodel.d
        self._trial_value = self.linemodel.obj(self.step, x=self.iterate)

        return self.step class TrustRegionSolver(object):
    """A generic class for implementing solvers for the trust-region subproblem.

    minimize q(d)  subject to  ||d|| <= radius,

    where q(d) is a quadratic function, not necessarily convex.

    The trust-region constraint `||d|| <= radius` can be defined in any
    norm although most derived classes currently implement the Euclidian
    norm only. Note however that any elliptical norm may be used via a
    preconditioner.

    For more information on trust-region methods, see

    A. R. Conn, N. I. M. Gould and Ph. L. Toint, Trust-Region Methods,
    MP01 MPS-SIAM Series on Optimization, 2000.
    """

    def __init__(self, qp, cg_solver, **kwargs):
        """Basic solver for a quadratic trust-region subproblem.

        :keywords:

          :qp: a ``QPModel`` instance
          :cg_solver: a `class` to be instantiated and used as solver for the
                      trust-region problem.
        """
        self._qp = qp
        self._cg_solver = cg_solver(qp, **kwargs)
        self._niter = 0
        self._step_norm = 0.0
        self._step = None
        self._m = None  # Model value at candidate solution
        self.status = ""

    @property
    def qp(self):
        """Return the ``QPModel`` instance."""
        return self._qp

    @property
    def niter(self):
        """Return the number of iterations made by the solver."""
        return self._niter

    @property
    def step_norm(self):
        """Return the norm of the step computed by the solver."""
        return self._step_norm

    @property
    def step(self):
        """Return the step computed by the solver."""
        return self._step

    @property
    def m(self):
        """Return the value of the quadratic model."""
        return self._m

    @property
    def model_value(self):
        """Return the value of the quadratic model."""
        return self._m

    def solve(self, *args, **kwargs):
        """Solve the trust-region subproblem."""
        self._cg_solver.solve(*args, **kwargs)
        self._niter = self._cg_solver.niter
        self._step_norm = self._cg_solver.step_norm
        self._step = self._cg_solver.step
        self._m = self.qp.obj(self.step)        # Compute model reduction.
        self.status = self._cg_solver.status class GeneralizedTrustRegion(TrustRegion):
    """A more general trust-region management class.

    Subclass of `TrustRegion` implementing a more flexible update of the
    trust-region radius.

    See :meth:`update_radius` for more details.
    """

    def __init__(self, **kwargs):
        """Instantiate an object allowing management of a trust region.

        :keywords:
            :radius:        Initial trust-region radius (default: 1.0)
        """
        self.radius = self.radius0 = kwargs.get('radius', 1.0)
        self.radius_max = 1.0e+10
        self.eta0 = 1e-4
        self.eta1 = 0.25
        self.eta2 = 0.75
        self.gamma1 = 0.25
        self.gamma2 = 0.5
        self.gamma3 = 4.0
        self.eps = np.finfo(np.double).eps  # Machine epsilon.

    def update_radius(self, ratio, step_norm, alpha):
        u"""Update the trust-region radius.

        The rule implemented by this method is:

        radius = min(max(α, γ1) * step_norm, γ2*radius)      if ratio <= η0
        radius = max(γ1*radius, min(α*step_norm, γ2*radius)) if ratio ∈ (η0,η1]
        radius = max(γ1*radius,
                     min(α*step_norm, γ3*self.radius))       if ratio ∈ (η1,η2]
        radius = max(radius, min(α*step_norm, γ3*radius))    if ratio > η2.
        """
        if ratio <= self.eta0:
            self.radius = min(max(alpha, self.gamma1) * step_norm,
                              self.gamma2 * self.radius)
        elif ratio <= self.eta1:
            self.radius = max(self.gamma1 * self.radius,
                              min(alpha * step_norm,
                                  self.gamma2 * self.radius))
        elif ratio <= self.eta2:
            self.radius = max(self.gamma1 * self.radius,
                              min(alpha * step_norm,
                                  self.gamma3 * self.radius))
        else:
            self.radius = max(self.radius, min(alpha * step_norm,
                                               self.gamma3 * self.radius)) def norm1(x):
def norm2(x):
def normp(x, p):
def norm_infty(x):
def normest(A, tol=1.0e-6, maxits=100):
        A = np.random.randn(n, m)
        A = np.random.rand(n, n)
        A = .5 * (A.T + A) def where(cond):
    """Bypass Numpy's annoyances.

    Gee does someone need to write a proper Numpy!
    """
    return np.where(cond)[0] def projected_gradient_norm2(x, g, l, u):
    """Compute the Euclidean norm of the projected gradient at x."""
    lower = where(x == l)
    upper = where(x == u)

    pg = g.copy()
    pg[lower] = np.minimum(g[lower], 0)
    pg[upper] = np.maximum(g[upper], 0)

    return norm2(pg[where(l != u)]) def project(x, l, u):
    """Project x into the box [l, u]."""
    return np.maximum(np.minimum(x, u), l) def projected_step(x, d, l, u):
    """Project the step d into the box [l, u].

    The projected step is defined as s := P[x + d] - x.
    """
    return project(x + d, l, u) - x def breakpoints(x, d, l, u):
    """Find the smallest and largest breakpoints on the half line x + t d.

    We assume that x is feasible. Return the smallest and largest t such
    that x + t d lies on the boundary.
    """
    pos = where((d > 0) & (x < u))  # Hit the upper bound.
    neg = where((d < 0) & (x > l))  # Hit the lower bound.
    npos = len(pos)
    nneg = len(neg)

    nbrpt = npos + nneg
    # Handle the exceptional case.
    if nbrpt == 0:
        return (0, 0, 0)

    brptmin = np.inf
    brptmax = 0
    if npos > 0:
        steps = (u[pos] - x[pos]) / d[pos]
        brptmin = min(brptmin, np.min(steps))
        brptmax = max(brptmax, np.max(steps))
    if nneg > 0:
        steps = (l[neg] - x[neg]) / d[neg]
        brptmin = min(brptmin, np.min(steps))
        brptmax = max(brptmax, np.max(steps))

    return (nbrpt, brptmin, brptmax) def cputime():
    """Return the user CPU time since the start of the process."""
    return resource.getrusage(resource.RUSAGE_SELF)[0] class UserExitRequest(Exception):
    """Exception that the caller can use to request clean exit."""

    pass class LineSearchFailure(Exception):
    """Exception raised when a linesearch fails."""

    pass
import logging
import numpy as np
from pykrylov.linop import SymmetricallyReducedLinearOperator as ReducedHessian
from nlp.model.nlpmodel import QPModel
from nlp.model.linemodel import C1LineModel
from nlp.ls.linesearch import ArmijoLineSearch
from nlp.tr.trustregion import TrustRegionSolver
from nlp.tr.trustregion import GeneralizedTrustRegion
from nlp.tools import norms
from nlp.tools.utils import where, projected_gradient_norm2, \
    project, projected_step, breakpoints
from nlp.tools.timing import cputime
from nlp.tools.exceptions import UserExitRequest, LineSearchFailure

            # Use a projected search to obtain the next iterate
            (xfree, proj_step) = self.projected_linesearch(x[free_vars],
                                                           l[free_vars],
                                                           u[free_vars],
                                                           gfree, step, ZHZ,
                                                           alpha=1.0)

            # Update the minimizer and the step.
            # Note that s now contains x[k+1] - x[0]
            x[free_vars] = xfree
            s[free_vars] += proj_step

            # Compute the gradient grad q(x[k+1]) = g + H*(x[k+1] - x[0])
            # of q at x[k+1] for the free variables.
            Hs = H * s
            gfree = g[free_vars] + Hs[free_vars]
            gfnormf = norms.norm2(gfree)

            # Convergence and termination test.
            # We terminate if the preconditioned conjugate gradient method
            # encounters a direction of negative curvature, or
            # if the step is at the trust region bound.
            if gfnormf <= cgtol * gfnorm:
                exitOptimal = True
                info = 1
            elif self.solver.status == "trust-region boundary active":
                #  infotr == 3 or infotr == 4:
                exitPCG = True
                info = 2
            elif iters >= itermax:
                exitIter = True
                info = 3

        self.log.debug("leaving projected_newton_step with info=%d", info)
        return (x, s, iters, info)

    def projected_linesearch(self, x, l, u, g, d, H, alpha=1.0):
        u"""Use a projected search to compute a satisfactory step.

        This step must satisfy a sufficient decrease condition for the
        quadratic

            q(s) = gᵀs + ½ sᵀHs,

        where H=Hᵀ and g is a vector. Given the parameter α, the step is

           s[α] = P[x + α d] - x,

        where d is the search direction and P the projection into the
        box [l, u]. The final step s = s[α] satisfies the sufficient decrease
        condition

           q(s) ≤ μ₀ gᵀs,

        where μ₀ ∈ (0, 1).

        The search direction d must be a descent direction for the quadratic q
        at x such that the quadratic is decreasing along the ray  x + α d
        for 0 ≤ α ≤ 1.
        """
        self.log.debug("performing projected linesearch")
        mu0 = 0.01
        interpf = 0.5
        nsteps = 0

        # Find the smallest break-point on the ray x + α d.
        (_, brptmin, _) = breakpoints(x, d, l, u)
        self.log.debug("nearest  breakpoint: %7.1e", brptmin)

        # Reduce alpha until the sufficient decrease condition is
        # satisfied or x + α w is feasible.

        search = True
        while search and alpha > brptmin:

            # Calculate P[x + alpha*w] - x and check the sufficient
            # decrease condition.
            nsteps += 1

            s = projected_step(x, alpha * d, l, u)
            Hs = H * s
            gts = np.dot(g, s)
            q = .5 * np.dot(Hs, s) + gts
            if q <= mu0 * gts:
                search = False
            else:
                alpha *= interpf

        # Force at least one more constraint to be added to the active
        # set if alpha < brptmin and the full step is not successful.
        # There is sufficient decrease because the quadratic function
        # is decreasing in the ray x + alpha*w for 0 <= alpha <= 1.
        if alpha < 1 and alpha < brptmin:
            alpha = brptmin

        # Compute the final iterate and step.
        s = projected_step(x, alpha * d, l, u)
        x = project(x + alpha * d, l, u)
        return (x, s)

    def solve(self):
        """Solve method.

        All keyword arguments are passed directly to the constructor of the
        trust-region solver.
        """
        self.log.debug("entering solve")
        model = self.model
        ls_fmt = "%7.1e  %8.1e"

        # Project the initial point into [l,u].
        self.x = project(self.x, model.Lvar, model.Uvar)

        # Gather initial information.
        self.f = model.obj(self.x)
        self.f0 = self.f
        self.g = model.grad(self.x)  # Current  gradient
        self.g_old = self.g.copy()
        self.x_old = self.x.copy()

        pgnorm = projected_gradient_norm2(self.x, self.g,