% [X,LE,LI,F,CE,CI,MU,ERR,Oit,condH,G] =
%   qpenalty(f,paramf,E,paramE,I,paramI,O,paramO,x0[,mu0,tol,maxit])
% Constrained optimisation by quadratic-penalty method
%
% Reference: ch. 17 in:
%   Nocedal and Wright: "Numerical Optimization", Springer, 1999.
%
% Example usage: the following will minimise Rosenbrock's function subject
% to 1 quadratic equality constraint and 2 quadratic inequality constraints
% (with respective parameters), using the newtoncg unconstrained optimiser,
% starting at the origin and running 50 iterations at most:
%
%   [X,LE,LI,F,CE,CI,MU] = qpenalty(@rosenbrock,{},...
%      {@quadf},{{A1,b1,c1}},...
%      {@quadf,@quadf},{{A2,b2,c2},{A3,b3,c3}},...
%      newtoncg,{'quadratic'},...
%      [0;0],1,[],50);
%   figure(1);
%   fcontours(@rosenbrock,{},...
%      {@quadf},{{A1,b1,c1}},...
%      {@quadf,@quadf},{{A2,b2,c2},{A3,b3,c3}},...
%      [-2 2;-1 3],50);
%   plotseq(X);
%
% In:
%   f: a handle to a function that takes as input a list of vectors and
%      returns 3 lists: the function values, the gradients and the Hessians
%      (see rosenbrock.m for an example in 2D).
%   paramf: cell array containing the parameters for f. Use {} to pass no
%      parameters, or to use f's default parameters.
%   E: a cell array {e1,...} of equality constraints. Each ei is a handle to
%      a function like f. Use {} if there are no equality constraints.
%   paramE: a cell array {p1,...} containing the parameters for e1,... (pi is
%      itself a cell array, like paramf). Use {} for pi to pass no parameters,
%      or to use pi's default parameters.
%   I, paramI: like E, paramE but for inequality constraints.
%   O: a handle to a function that implements an unconstrained optimisation
%      method, with the same argument format as newtoncg.m.
%   paramO: cell array containing the parameters for O. Use {} to pass no
%      parameters, or to use O's default parameters. Do not include x0, tol,
%      maxit (they are set by qpenalty).
%   x0: n x 1 vector, the starting point.
%   mu0: (positive) initial value of the penalty parameter. Default: 1.
%   tol: tolerance to stop iterating. Default: 1e-4.
%   maxit: maximum number of iterations. Default: 1000.
% Out:
%   X: ? x n list of iterates for the primal variables.
%   LE: ? x le list of estimates for the Lagrange multipliers of the equality
%      constraints.
%   LI: ? x li list of estimates for the Lagrange multipliers of the inequality
%      constraints.
%   F: ? x 1 list of the function value at the iterates.
%   CE: ? x le list of the equality constraints value at the iterates.
%   CI: ? x li list of the inequality constraints value at the iterates.
%   MU: ? x 1 list of the penalty parameter value at the iterates.
%   ERR: ? x 1 list of the error in the KKT system at the iterates.
%   Oit: ? x 1 list of the number of iterations in the unconstrained
%      optimisation method at each iterate.
%   condH: ? x 1 list of the condition number of the Hessian of the
%      quadratic-penalty function at the iterates.
%   G: ? x 1 list of the norm of the objective function gradient at the
%      iterates.
%
% Any non-mandatory argument can be given the value [] to force it to take
% its default value.

% Copyright (c) 2005 by Miguel A. Carreira-Perpinan

function [X,LE,LI,F,CE,CI,MU,ERR,Oit,condH,G] = ...
      qpenalty(f,paramf,E,paramE,I,paramI,O,paramO,x0,mu0,tol,maxit)

% ---------- Argument defaults ----------
if ~exist('mu0','var') | isempty(mu0) mu0 = 1; end;
if ~exist('tol','var') | isempty(tol) tol = 1e-4; end;
if ~exist('maxit','var') | isempty(maxit) maxit = 1000; end;
% ---------- End of "argument defaults" ----------

% STUFF HERE

  % INSIDE THE LOOP FOR mu_k: the call to the unconstrained solver of Q, with
  % - x: the starting point
  % - tau: your tau_k
  % - maxit2: maximum number of iterations
  [X1,F1] = O(@Q,{f,paramf,E,paramE,I,paramI,mu},x',tau,maxit2,paramO{:});
  x = X1(end,:);
  
% STUFF HERE

