#   =================================================================
#   File: algencan.py
#   =================================================================
#
#   =================================================================
#   Module: Python module
#   =================================================================
#
#   Last update of any of the component of this module:
#
#   August 06, 2007.
#
#   Users are encouraged to download periodically updated versions of
#   this code at the TANGO home page:
#
#   www.ime.usp.br/~egbirgin/tango/
#
#   ******************************************************************
#   ******************************************************************

from numpy import *

#   Load the solver wrapper

from pywrapper import solver

def solvers(evalf,evalg,evalh,evalc,evaljac,evalhc,evalhlp,inip,endp):
    """Call the solver."""

    solver(evalf,evalg,evalh,evalc,evaljac,evalhc,evalhlp,inip,endp,param)

#   Parameters of the solver:
#   =========================
#
#   RHOAUTO logical
#   ---------------
#
#   indicates whether the initial penalty parameters will be the
#   ones given by the user (rhoauto = .false.) or will be computed
#   automatically by the Augmented Lagrangian solver (rhoauto =
#   .true.).
#
#   RHOTYPE integer
#   ---------------
#
#   indicates whether the penalty parameters will be updated all
#   toghether (rhotype=1) or individualy (rhotype=2)
#
#   RHOMULT double precision
#   ------------------------
#
#   when a penalty parameter is updated, it is multiplied by rhofrac.
#   So, rhofrac must be greater than 1. Suggested value rhofrac=10.0d0.
#
#   RHOFRAC double precision
#   ------------------------
#
#   this paramater between 0 and 1 is the improvement required in the
#   infeasibility for not to update the penalty parameters.
#
#   When rhotype=1, if the sup-norm of the constraints is smaller than
#   or equal to rhofrac times the sup-norm of the constraints of the
#   previous iteration then the penalty parameters are not updated.
#   Otherwise, they are multiplied by rhomult.
#
#   When rhotype=2, the improvement in feasibility of each constraint
#   is evaluated in separate. If the feasibility of a particular
#   constraint is smaller than or equal to rhofrac times the
#   feasibility of the previous iteration then the penalty parameter
#   is not modified. Otherwise, it is multiplied by rhomult.
#
#   The suggested value for rhofrac is 0.5d0.
#
#   RHO double precision rho(m)
#   ---------------------------
#
#   This the vector of penalty parameters. The user can leave to the
#   method the task to find adequate initial values for the penalty
#   parameters setting rhoauto = .true. Otherwise, he/she can set
#   rhoauto = .false. and atribute a value for each one of the penalty
#   parameters rho(i), i = 1, ..., m,
#
#   GTYPE integer
#   -------------
#
#   Type of first derivatives calculation according to the following
#   convention:
#
#   0 means true first derivatives. In this case, functions evalg
#     and evaljac must be modified by the user to compute the
#     gradient of the objective function and the gradients of the
#     constraints, respectively.
#
#   1 means that a finite difference approximation will be used. In
#     this case, functions evalg and evaljac may have an empty body
#     but must be present. It is also recommended that those
#     empty-body functions set flag = - 1. Last but not least, the
#     option gtype = 1 is not cheap neither safe.
#
#   HPTYPE integer
#   --------------
#
#   Type of Hessian-vector product according to the following
#   convention:
#
#   We will first describe the meaning if the choices in the Augmented
#   Lagrangian framework. The way in which the product of the Hessian
#   matrix of the Augmented Lagrangian by a vector will be done
#   depends on the value of the parameter hptype in the following way:
#
#   9 means that an incremental quotients approximation without any
#     extra consideration will be used. This option requires the
#     evaluation of an extra gradient at each Conjugate Gradient
#     iteration. If gtype = 0 then this gradient evaluation will be
#     done using the user supplied functions evalg and evaljac
#     (evalc will also be used). On the other hand, if gtype = 1, the
#     gradient calculation will be done using just calls to the user
#     provided functions evalf and evalc. nind calls will be done,
#     where nind is the dimension of the current face of the
#     active-set method. This option is not cheap neither safe.
#
#     If you did not code functions evalg and evaljac, to compute
#     the gradient of the objective function and the Jacobian of the
#     constraints then your options finished here.
#
#   0 means that functions to compute the Hessian of the objective
#     function (evalh) and the Hessians of the constraints (evalhc)
#     were provided by the user. So, the product of the Hessian of the
#     Augmented Lagrangian times a vector will be computed using the
#     Hessians provided by these functions and then adding the first
#     order term (for the first order term the user-supplied
#     function to compute the Jacobian of the constraints (evaljac)
#     is also used).
#
#   1 means that, instead of providing individual functions to
#     compute the Hessians of the objective function and the
#     constraints, the user provided a function to compute the
#     product of an arbitrary vector times the Hessian of the
#     Lagrangian.
#
#   2 means that incremental quotients will be used. The difference
#     between hptype = 9 and hptype = 2 is that, in the latter case,
#     the non-differentiability of the Hessian of the Augmented
#     Lagrangian will be taken into account. In particular, when
#     computing the gradient of the Augmented Lagrangian at
#     (x + step p), the constraints that will be considered will be
#     the same constraints that contributed to the computation of the
#     gradient of the Augmented Lagrangian at x.
#
#     If GENCAN is been used to solve a bound-constrained problem
#     which is not the subproblem of an Augmented Lagrangian method,
#     but an isolated bound-constrained problem, then there is no
#     difference between this option and hptype = 9.
#
#     This option also requires the evaluation of an extra gradient at
#     each Conjugate Gradient iteration.
#
#   3 is similar to hptype = 2. The difference is that the
#     contribution of the linear constraints in the Hessian matrix
#     will be computed explicitly. If the problem has not linear
#     constraints then this option is identical to hptype = 2.
#     Moreover, if the problem has no constraints then this option is
#     equal to hptype = 9.
#
#     This option also requires the evaluation of an extra gradient at
#     each Conjugate Gradient iteration.
#
#   4 means that the Hessian matrix will be approximated and then the
#     product of the Hessian approximation by the vector will be
#     computed exactly. In particular, the Hessian matrix will be
#     approximated doing a BFGS correction to the Gauss-Newton
#     approximation of the Hessian. Before the BFGS correction, a
#     structured spectral correction is done to force the Gauss-Newton
#     approximation to be positive definite.
#
#     If the problem has not constraints then the approximation
#     reduces to a BFGS approximation of the Hessian (without memory)
#     and using the spectral approximation (instead of the identity)
#     as initial approximation.
#
#     Numerical experiments suggested that this option is convenient
#     just for constrained problems. This motivated the introduction
#     of the next option.
#
#     This option does NOT require an extra gradient evaluation per
#     iteration and, in this sense, each CG iteration is
#     computationally cheaper than a CG iteration of the previous
#     choices. However, the approximation of the Hessian matrix
#     requires some information (mainly the Jacobian of the
#     constraints) that must be saved during the gradient evaluation.
#     To save this information requires an amount of memory
#     proportional to the number of non-null elements of the Jacobian
#     matrix.
#
#     Quadratic subproblems are convex with this choice.
#
#   5 is an adaptive strategy that choose, at every iteration, between
#     2 and 4. When the gradient of the Augmented Lagrangian is
#     computed, it is verified if at least a constraint contributes to
#     the calculation. If this is the case, 4 is used. Otherwise, 2 is
#     used.
#
#     For problems with equality constraints (that always contributes
#     to the Augmented Lagrangian function) this option is identical
#     to 4.
#
#     For problems without constraints this option is identical to 2.
#
#   6 is identical to 5 but the choice is made between 3 and 4 instead
#     of between 2 and 4.
#
#     For problems with equality constraints (that always contributes
#     to the Augmented Lagrangian function) this option is identical
#     to 4.
#
#     For problems without constraints this option is identical to 3.
#
#   We will now describe the meaning if the choices for unconstrained
#   and bound-constrained problems. In this context the way in which
#   the product of the Hessian matrix by a vector will be done depends
#   on the value of the parameter hptype in the following way:
#
#   0 means that the function to compute the Hessian of the
#     objective function (evalh) was provided by the user. So, the
#     product of the Hessian times a vector will be computed using the
#     Hessian provided by this function.
#
#   1 means that a function (evalhlp) to compute the product of the
#     Hessian of the objective function times an arbitrary vector is
#     being provided by the user.
#
#   9 means that an incremental quotients approximation will be used.
#     This option requires the evaluation of an extra gradient at each
#     Conjugate Gradient iteration. If gtype = 0 then this gradient
#     evaluation will be done using the user supplied function
#     evalg. On the other hand, if gtype = 1, the gradient calculation
#     will be done using just calls to the user provided function
#     evalf. nind calls will be done, where nind is the dimension of
#     the current face of the active-set method.
#
#     If you did not code function evalg to compute the gradient of
#     the objective function then your options finished here.
#
#   4 means that the Hessian matrix will be approximated and then the
#     product of the Hessian approximation by the vector will be
#     computed exactly. The approximation is a BFGS approximation of
#     the Hessian (without memory) and using the spectral
#     approximation (instead of the identity) as initial
#     approximation.
#
#     Numerical experiments suggested that this option is not
#     convenient for unconstrained or just bound-constrained problems.
#     (Note that this option was developed to be used in the Augmented
#     Lagrangian framework.)
#
#     This option does NOT require an extra gradient evaluation per
#     iteration and, in this sense, each CG iteration is
#     computationally cheaper than a CG iteration of the previous
#     choices.
#
#     Quadratic subproblems are convex with this choice.
#
#   In the bound-constrained context, options hptype = 2,3,5 and 6
#   are all identical to hptype = 9.
#
#   INTYPE integer
#   --------------
#
#   This parameter can be used to select the algorithm that will be
#   used to solve the Augmented Lagrangian subproblems, or the problem
#   itself if it an unconstrained problem. In the present implementation
#   just one algorithm is available. So, the value of this variable is
#   ignored.
#
#   PRECOND character * 6
#   ---------------------
#
#   Indicates the type of preconditioning that will be used for
#   Conjugates Gradients according to the following convention:
#
#   'NONE'   means no preconditioner at all.
#
#   'QNAGNC' means Quasi-Newton Correction of the Gauss-Newton
#            approximation of the Hessian. The exact form is this
#            preconditioner is described in:
#
#            E. G. Birgin and J. M. Martinez, "Structured minimal-
#            memory inexact quasi-Newton method and secant
#            preconditioners for Augmented Lagrangian Optimization",
#            submitted, 2005.
#
#   CHECKDER logical
#   ----------------
#
#   If you are using finite differences aproximations for the
#   derivatives (gtype = 1) then checkder must be set to FALSE.
#   On the other hand, if you are using your own coded derivatives
#   you may would like to test them against a finite differences
#   approximation to verify if they are o.k. In this case, set
#   checkder = TRUE.
#
#   EPSFEAS double precision
#   ------------------------
#
#   Feasibility tolerance for the sup-norm of the constraints.
#   (Ignored in the unconstrained and bound-constrained cases.)
#
#   EPSOPT double precision
#   -----------------------
#
#   Optimality tolerance for the sup-norm of the projected gradient of
#   the Augmented Lagrangian in the constrained case and the sup-norm
#   of the projected gradient of the objective function in the
#   unconstrained and the bound-constrained cases.
#
#   MAXOUTIT integer
#   ----------------
#
#   Maximum number of Augmented Lagrangian (outer) iterations.
#   (Ignored in the unconstrained and bound-constrained cases.)
#
#   MAXTOTIT integer
#   ----------------
#
#   Maximum total number of inner iterations in the Augmented
#   Lagrangian context (total means summing up the inner iterations of
#   each outer iteration). In the unconstrained and bound-constrained
#   cases it means just the maximum number of iterations.
#
#   MAXTOTFC integer
#   ----------------
#
#   Idem MAXTOTIT but for number of functional evaluations.
#
#   IPRINT integer
#   --------------
#
#   Controls the ammount of information of the output according to the
#   following convention:
#
#   0 means no output at all.
#
#   1 means information at each outer iteration but without any
#     information of how the subproblems are being solved.
#
#   2 means the same as 1 plus information of each inner iteration.
#
#   3 means the same as 2 plus information of the line searches and
#     the calculation of the truncated Newton direction (using CG) of
#     each inner iteration.
#
#   In all cases, an output file named solution.txt with the final
#   point, Lagrange mutipliers and penalty parameters will be
#   generated. Moreover, the same output of the screen will be saved
#   in a file named algencan.out.
#
#   NCOMP integer
#   -------------
#
#   Every time a vector is printed, just its first ncomp component
#   will be displayed.

param = {
    'gtype' : 0,
    'hptype': 6,

    'intype': 0,

    'precond': 'QNCGNA',

    'rhoauto': True,

    'rhotype': 1,
    'rhomult': 1.0e+01,
    'rhofrac': 0.5,

    'rho': ones(2),

    'checkder': False,

    'epsfeas': 1.0e-04,
    'epsopt' : 1.0e-04,

    'maxoutit': 50,
    'maxtotit': 1000000,
    'maxtotfc': 5000000,

    'iprint': 1,
    'ncomp' : 5
}