% file: Optimize.m
% description: designed to performs experiments in minimization
% that allow user to test various methods on a user defined problem.
% last rev: 12/15/05
%

% PREPROCESSION SECTION

% PROBLEM DEFINITION:
fn = 'extrosen'; % options: your favorite function to be minimized, or in
                 % case of zero or least squares, function 1/2*||r(x)||^2
fngrad = 'extrosengrad'; % its (possibly approximate) gradient function
fnhess = 'extrosenhess'; % its (possibly approximate) Hessian function
% following functions are used for least squares problems
fnvec = ''; % vector function r(x), for least squares or zero problems
fnjac = ''; % jacobian matrix J(x), for least squares or zero problems
% following two functions are used to define a least squares problem
% generated by a nonlinear least squares problem and are needed by the
% template functions leastsq*.m, which must be entered in all fns above.
datafitfn = ''; % defines a single entry in residual vector 
datafitjac = '';% defines single row in jacobian matrix of residual vector
%
% METHOD DEFINITION:
% options for global method:
% newton, bfgs, steepdesc, trustregion, cg, gaussnewton, broyden, none
globalmtd = 'trustregion';
% 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 = 'cgsteihaug';
% options for stopping method: quadstop, superlinstop, linstop
stopmtd = 'quadstop';
%

% 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:
n = 20; % a dimension parameter, if needed for problem definition
% vectors here should be column vectors
%xnew = [50,10,.1]'; % logistic guess
%xnew = zeros(n,1); % Final exam qu #5
xnew = linspace(-n/2,n/2,n)'; % Final exam qu #6
%xnew = [-0.5 1.4]'; % initial guess
%xnew = (sqrt(1:n)'-sqrt(n)/2).*(rand(n,1)-0.5); % initial guess
%xnew = repmat([-1.2,1],1,n/2)'; % initial guess
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:
outputdata = zeros(parms.maxiter,4);
% user must define the global parms.data here, if needed,
% e.g., run a data generating function as in next line
% parms.data = mydata();
% parms.datasize = size(parms.data,1);
% parms.datafn = datafitfn;
% parms.datajac = datafitjac;

% USER INPUT GLOBAL PARAMETERS
% following variables needed for linear systems via quadfn
%global quadA; 
%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);
%

%
% user data input ends here...but printing and other
% post-processing may require user attention below
%
% BFGS initialization
if (strcmp(globalmtd, 'bfgs'))
  global H; % needed for efficient BFGS
  H = eye(length(xnew));
end
% Broyden initialization
if (strcmp(globalmtd, 'broyden'))
  global H; % needed for efficient Broyden
  H = feval(fnjac,xnew);
end
% 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
%

% PROCESSING SECTION

if (strcmp(globalmtd, 'none'))
  return;
end
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)
      outputdata(countiter,:) = [countiter, alpha, ...
                    norm(xnew-xtrue),abs(feval(fn,xnew)-feval(fn,xtrue))];
    else
      outputdata(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
outputdata = outputdata(1:countiter,:); % trim off excess zeros
disp(outputdata);
format
%
