% file: Optimize.m
% description: performs adaptable experiments in unconstrained minimization.
% This is a driver template.  It is currently set to solve an optimization
% problem which comes from Exercise 9.1 of Nocedal and Wright's text,
% Numerical Optimization.  This is the extended Rosenbrock function.
% 
% programmer: t. shores                     last rev: 5/06/05
%

% PREPROCESSION SECTION

% PROBLEM DEFINITION:
fn = 'extrosenfn'; % options: your favorite function to be minimized, or in
               % case of zero or least squares, function 1/2*||r(x)||^2
fngrad = 'extrosenfngrad'; % its (possibly approximate) gradient function
fnhess = 'extrosenfnhess'; % its (possibly approximate) Hessian function
fnvec = 'leastsqvec'; % vector function r(x), for least squares or zero problems
fnjac = 'leastsqjac'; % jacobian matrix J(x), for least squares or zero problems
n = 20; % a dimension parameter, if needed for problem definition
%

% METHOD DEFINITION:
% options for global method:
% newton, bfgs, steepdesc, trustregion, cg, gaussnewton, broyden, none
globalmtd = 'bfgs';
% options for local method with newton, bfgs, steepdesc, gaussnewton, broyden:
% backtrack, wolfelinesrch, conststep
% options for local method with trustregion: dogleg, cgsteihuag
% options for local method with cg (conjugate gradient): 
% lincg (for linear cg), prpluscg(uses wolfelinesrch)
localmtd = 'wolfelinesrch';
% options for stopping method: quadstop, superlinstop, linstop
stopmtd = 'superlinstop';
%

% GLOBAL PARAMETERS DEFINITIONS:
% since Matlab/Octave don't have protected variables,
% we follow this convention: global parameters in this
% parameter object will be modified in this file and
% *only* in this file, if at all. These parameters
% will be visible to any routine that wants them.
%
global parms
%

% PARAMETERS FOR A LINE SEARCH METHOD: 
parms.wolferho = 1.4; % wolfelinesrch growth factor
parms.backtrackrho = 0.8; % backtrack shrink factor
parms.constrho = 1; % constant step value
parms.c1 = 1e-4; % Wolfe 1 (Armijo)
parms.c2 = 0.9;  % Wolfe 2, use .9 for Newton, lower if line search failure
parms.alpha1 = 1; % initial step for wolfelinesrch and backtrack
parms.alphamax = 2; % maximum allowable alpha for wolfelinesrch
parms.const = 1; % constant stepsize for conststep
%

% PARAMETERS FOR A TRUST REGION METHOD:
parms.deltacurr = 1; % current delta
parms.deltamax = 2; % maximum allowable delta
parms.eta = 1e-3; % accept/reject step measure..must be in [0,1/4)
parms.rhomin = 1/4; % below which region shrinks
parms.deltashrink = 1/2; % delta shrink factor
parms.rhomax = 3/4; % above which region grows
parms.deltagrow = 2; % delta growth factor
parms.maxcgiter = 10; % maximum number of iterations for CG in Steihaug
parms.cgeps = 1e-6; % stopping tolerance for CG in Steihaug
%

% PARAMETERS FOR A CONJUGATE GRADIENT METHOD:
parms.cgwolferho = 1.4; % wolfelinesrch growth factor
parms.cgc1 = 1e-4; % Wolfe 1 (Armijo)
parms.cgc2 = 0.1;  % Wolfe 2
parms.cgalpha1 = 1; % initial step for wolfelinesrch
parms.cgalphamax = 5; % maximum allowable alpha for wolfelinesrch
parms.cgrestart = 0.1; % restart parm for prpluscg--typical angle cos:0.1,period n
%

% PARAMETERS FOR A GAUSS-NEWTON METHOD:
parms.gnmin = 1e-8; % relative lower bound for minimum singular value
%

% PARAMETERS FOR A STOPPING CRITERION:
parms.maxiter = 200; % maximum number of iterations of global routine
parms.releps = 1e-6; % relative stopping tolerance
parms.abseps = 1e-10; % absolute stopping tolerance
parms.ordcoef = .99; % order coefficient for any stopping method
parms.hstep = 1e-5; % for approximate gradients, Jacobians and Hessians
%

% PARAMETERS FOR INITIALIZATION:
% vectors here should be column vectors
xnew = -ones(n,1); % starting point
xtrue = ones(n,1); % exact soln if known
solnflag = 1; % set to true (1) if exact solution known
% add extra locals, e.g., for printout, here:
data = zeros(parms.maxiter,4);
%

% USER INPUT GLOBAL PARAMETERS
%global quadA; % needed for linear systems via quadfn
%global quadb;
% form the matrix
%x = [140 10 10 10 10 10 10 120 10 1 10 10 1 10]';
% smudge it a bit
%x = x.*(1+0.001*rand(n,1));
%z = rand(n,1)-0.5;
%Hm = eye(n) - 2*(z*z')/(z'*z);
%quadA = Hm'*diag(x)*Hm;
%quadb = quadA*xtrue;
%quadA = hilb(n);
%quadb = ones(n,1);
%H = feval(fnjac,xnew); % Broyden initialization
%global parmfitfn; % for definition of one r_i(x), when homogeneous (least sq)
%parmfitfn = 'logisticfn';
%global parmfitjac; % definition of one r_i(x), when homogeneous (least sq)
%parmfitjac = 'logisticfnjac';
%global myP;
%myparms;
%

%
% user data input ends here...but printing and other
% post-processing may require user attention below
%

% more useful globals which may be modified in other m files:
global countiter; % needed for stopping or restarting count
countiter = 0; % updated by stopping routines
global H; % needed for efficient BFGS or Broyden
H = eye(length(xnew)); % BFGS initialization
%

% PROCESSING SECTION

strng = sprintf('Function: %s',fn);
disp(strng);
strng = sprintf('Global method: %s', globalmtd);
disp(strng);
strng = sprintf('Local method: %s', localmtd);
disp(strng);
strng = sprintf('Stopping criterion: %s', stopmtd);
disp(strng);
strng = sprintf('Size of the system: %d', size(xnew,1));
disp(strng);
if (solnflag)
  disp('Exact solution (transpose): ')
  disp(xtrue')  
else
  disp(' ')
end
doneflag = 0;
gradnew = feval(fngrad,xnew);
fnnew = feval(fn,xnew);
p = -gradnew; % for cg methods
while (~doneflag)
  xold = xnew;
  gradold = gradnew;
  fnold = fnnew;
  acceptflag = 1; % true if this step is accepted, thus adding new data
  switch globalmtd
    case 'trustregion'
      % do a single step
      [xnew,gradnew,alpha,acceptflag] = feval(globalmtd,localmtd,fnhess,fngrad,...
                               fn,gradold,xold);
      parms.deltacurr = alpha; % set the new current delta, called alpha 
      if (acceptflag) 
        % check for done
        doneflag = feval(stopmtd,xold,xnew);
      end
      % set up next step
      fnnew = feval(fn,xnew);
    case {'newton', 'bfgs', 'steepdesc'}
      [xnew,gradnew,alpha,p] = feval(globalmtd,localmtd,fnhess,fngrad, ...
                               fn,gradold,xold);
      % set up next step
      fnnew = feval(fn,xnew);
      %initial step size for linesearch, safeguarding 1, changing if sensible
      alphatmp = min(1, 2.02*(fnnew - fnold)/(gradnew'*p));
      if (alphatmp > 0)
        parms.alpha1 = alphatmp; % change alpha
      end
      % check for done
      doneflag = feval(stopmtd,xold,xnew);
    case {'cg'}
      [xnew,gradnew,alpha,p] = feval(globalmtd,localmtd,fnhess,fngrad, ...
                               fn,gradold,xold,p);
      % check for done
      doneflag = feval(stopmtd,xold,xnew);
      % set up next step
      fnnew = feval(fn,xnew);
    case {'gaussnewton'}
      [xnew,gradnew,alpha,p] = feval(globalmtd,localmtd,fnjac,fnvec,...
                                 fngrad,fn,gradold,xold);
      % set up next step
      fnnew = feval(fn,xnew);
      %initial step size for linesearch, safeguarding 1, changing if sensible
      alphatmp = min(1, 2.02*(fnnew - fnold)/(gradnew'*p));
      if (alphatmp > 0)
        parms.alpha1 = alphatmp; % change alpha
      end
      % check for done
      doneflag = feval(stopmtd,xold,xnew);
    case {'broyden'}
      [xnew,gradnew,alpha,p] = feval(globalmtd,localmtd,fnjac,fnvec,...
                                 fngrad,fn,gradold,xold);
      % set up next step
      fnnew = feval(fn,xnew);
      %initial step size for linesearch, safeguarding 1, changing if sensible
      alphatmp = min(1, 2.02*(fnnew - fnold)/(gradnew'*p));
      if (alphatmp > 0)
        parms.alpha1 = alphatmp; % change alpha
      end
      % check for done
      doneflag = feval(stopmtd,xold,xnew);
    otherwise
      disp('method unknown')
      break;
    end 
  % if you want to modify the output data, put commands here:
  if ((countiter > 0) & (acceptflag))
    if (solnflag)
      data(countiter,:) = [countiter, alpha, ...
                    norm(xnew-xtrue),abs(feval(fn,xnew)-feval(fn,xtrue))];
    else
      data(countiter,:) = [countiter, alpha, ...
                    norm(xnew-xold),abs(feval(fn,xnew)-feval(fn,xold))];
    end
  end
end

% POSTPROCESSING SECTION

% postprocessing commands, e.g., printing commands below, go here:
if (solnflag)
  disp('     k         alpha/delta   ||x(k+1)-xtrue||  |f(x(k+1)-f(xtrue)|');
else
  disp('     k         alpha/delta   ||x(k+1)-x(k)||  |f(x(k+1)-f(x(k))|');
end
format short e
data = data(1:countiter,:); % trim off excess zeros
disp(data);
format