function [x,fx,g,H,out] = mwd(fin,x0,options)
% smooth minimzation with simple bounds without using derivatives
%
% calling routine, for example:
%   [x,y,g,H,out] = mwd(@f_rosen, [-1;1;1])
%
% given a smooth function f=fin: R^n --> R 
%
% minimize  f(x)  for  lb <= x <= ub 
%
% f is not necessarily defined for x \not\in [lb,ub]
%
% Florian Jarre, March 2015
% Test Version with errors -- No guarantee of any kind is given!!!
%
% Mandatory input: 
%   A function handle fin ,
%   A starting point x0, 
%
% Optional input:
%   options.lb    -- a  lower bound on the variable x (Default -infty)
%                    the dimension must be (n,1)
%   options.ub    -- an upper bound on the variable x (Default +infty)
%                    the dimension must be (n,1)
%   options.maxit -- bound on the number of iterations; each iteration
%                    with about 2*n+15 function evaluations (Default 100*n)
%   options.update - update gradient central differences and Hessian by:
%                    default 1:  PSB only
%                            2:  PSB and curvature, unit basis
%                            3:  PSB and curvatrue, randomized bases
%                            4:  curvature only, randomized bases
%                            5:  curvature only, unit basis
%                            6:  zero Hessian (steepest descent)
%                                (2,4,5, and 6 are for test purposes only)
%   options.err   -- if the absolute error for a typical evaluation of
%                    fin is known, this parameter can be set here, else it
%                    is estimated (and used for the finite differences).
%
% Output:
%   x:         an approximate local minimizer 
%   fx:        the associated function value
%   g:         the final gradient approximation
%   H:         the final Hessian approximation
%   out.iter:  number of iterations needed
%   out.fval:  number of function evaluations needed
%   out.fline: average number of function evaluations in line search
%   out.err:   measurement error assumed for the computation of fin ---
%              this is NOT the estimated final accuracy of x but a quantity
%              used to define the step lengths for the finite difference
%              approximations. (Removed in later versions of mwd)
% 
% subroutines used:
%    updategh    - Update gradient and Hessian
%    curvilp     - Evaluate f at a projected point along a certain curve
%    rando       - generate a random orthogonal matrix
%    mwd1        - Line search routine -- in a separate .m-file -- !!!
%
% Algorithm:
% Trust region algorithm with finite difference approximation of the 
% gradient using a PSB update for the Hessian followed by a curvature
% correction based on randomized ON-bases for the central differences.
% Assume that the evaluation cost of f is at least O(n^2). (This would 
% be the case for quadratic functions with dense Hessians!)

%RandStream.setDefaultStream(RandStream('mt19937ar','seed',100));
RandStream.setGlobalStream(RandStream('mt19937ar','seed',100));

f = @(x) fin(x);

% COMPLETE INPUT ARGUMENTS
if (nargin < 2)
   error('starting point must be supplied');
end
n = length(x0);

if (nargin < 3)
   options.lb = -Inf*ones(n,1);
   options.ub =  Inf*ones(n,1);
   options.maxit = 100*n;
   options.update = 1;
else
   if ~isfield(options,'lb')
      options.lb = -Inf*ones(n,1);
   end
   if ~isfield(options,'ub')
      options.ub =  Inf*ones(n,1);
   end
   if ~isfield(options,'maxit')
      options.maxit = 100*n;
   end
   if ~isfield(options,'update')
      options.update = 1;
   end
end
lb = options.lb;
ub = options.ub;
o_up = options.update;

if options.maxit < 2
    disp('enforce at least two iterations');
    options.maxit = 2;
end
maxit = options.maxit;

if min(ub-lb) < 0
   error('bounds are not consistent');
end

if size(x0,2) > size(x0,1)   % make x0 a column vector
   x0 = x0';
end






% PREPARE FOR MINIMIZATION

iter = 0;
dt = (1+norm(x0))*eps^(1/2);
dtt = 2*dt;

xact = max(min(x0,ub-dtt),lb+dtt);
fact = feval(f,xact); iter = iter + 1;     



% estimate gradient and the error of the function evaluation of f near x0 
n = length(xact);
U = rando(n); 
fval = zeros(n,2);

err = eps;
g = zeros(n,1);

for i = 1:n
   fval(i,1) = f(xact+dt*U(:,i));
   fval(i,2) = f(xact+dtt*U(:,i));
   err = max(err, abs(fval(i,1)-0.5*(fact+fval(i,2))));
   g(i) = (fval(i,2)-fact)/dtt;
end
iter = iter+2*n;
if ~isfield(options,'err')
   options.err = err;
else
   if options.err < 0.1*err
      disp('either the function has large second derivative or')
      disp('options.err was too optimistic; err is being increased')
      options.err = err;
   end
   if options.err > 10*err
      disp('options.err was more pessimistic than the estimate generated')
      disp('by the algorithm; the pessimistic estimate is used');
   end
   err = options.err;
end
out.err = err;
dt   = 0.1*err^(1/3);
ub2  = ub-dt; lb2 = lb+dt; % search within the smaller set to allow for
                           % central derivatives
xact = max(min(x0,ub2),lb2);                          

%redo the gradient and the Hessian
H = zeros(n);
if err < 0.001*dtt 
   g = U*g; % with about >= 3 digits accuracy
else
   fval = zeros(n,2);
   for i = 1:n
      fval(i,1) = f(xact-dt*U(:,i));
      fval(i,2) = f(xact+dt*U(:,i));
   end
   iter = iter+2*n;
   g = U*(fval(:,2)-fval(:,1))/(2*dt);
   if o_up > 1 && o_up < 6
      dbeta = (sum(fval,2)-2*fact*ones(n,1))/(dt^2);
      H = U*diag(dbeta)*U'; % curvature correction
      H = 0.5*(H+H');
   end
end

optL.x0   = 0;         % lower bound for curvi-linear search
optL.x1   = 1;         % upper bound for curvi-linear search
optL.err  = err;
optL.tol  = 1.0e-6;
converged = 0;         % stopping the loop
outerit = 1;





while converged == 0  && outerit < maxit % LOOP

   outerit = outerit + 1;
   ng = -g;
   iact = (xact-lb2<1.1*dt & ng<0) + (ub2-xact<1.1*dt & ng>0); %   active 
   inact = ~iact;                                      % inactive indices
   if max(inact)
      Hact = H(inact,inact);
   else
      Hact = H;
   end
   [V,D] = eig(Hact);  % unitary V and diagonal D so that V'*Hact*V = D
   ocl = 1; % curvilp shall return the function value, too
   [tmin,fnew,out1] = ...
      mwd1(@(t) curvilp(f,t,xact,V,D,ng,lb2,ub2,inact,ocl),optL);
   if fnew > fact
      converged = 1;
   end
   if tmin == 0
      converged = 1;
   end

   iter = iter + out1.it;
   
   ocl = 0; % the function value is not needed here
   [~,xnew] = curvilp(f,tmin,xact,V,D,ng,lb2,ub2,inact,ocl);
   fold = fact; fact = fnew;
   
   [g,H] = updategh(f,xnew,xact,g,H,dt,o_up);
   iter  = iter+2*n+1;

   steplength = norm(xnew-xact);
   xact = xnew;
   if converged == 0
      converged = steplength < 0.01*dt && fold - fact < 0.01*dt^2;
   end
   
end % OF LOOP



% DO A FINAL  STEP WITH ORIGINAL BOUNDS

ng = -g;
iact = (xact-lb2<dt & ng<0) + (ub2-xact<dt & ng>0); %   active Indices
inact = ~iact;                                      % inactive indices
if max(inact)
   Hact = H(inact,inact);
else
   Hact = H;
end
[V,D] = eig(Hact);   % unitary V and diagonal D so  that V'*H*V = D

ocl = 1;
[tmin,fnew,out1] =mwd1(@(t) curvilp(f,t,xact,V,D,-g,lb,ub,inact,ocl),optL);
iter      = iter + out1.it;
ocl = 0;
[~,x]     = curvilp(f,tmin,xact,V,D,ng,lb,ub,inact,ocl);
fx        = fnew; 
g         = g+H*(x-xact); % estimate the final gradient
out.iter  = outerit;
out.fval  = iter;
out.fline = iter/outerit-2*n-1; %2*n+1 function evaluations for gradient
end








function [gnew,Hnew] = updategh(f,xnew,xact,g,H,dt,o_up)
% update gradient and Hessian of f
%
% Input:
%        a function handle f
%        xnew: the point at which g and H are searched for
%        xact: the point at which approximations of g and H are given
%        g, H: approximations of gradient and Hessian of f at xact
%        dt: step length for finite difference approximation
%
% Output:
%        gnew: approximation of gradient at xnew
%        Hnew: approximation of Hessian at xnew
%
% Florian Jarre, Feb. 22, 2015
%
% Use PSB and finite differences to update gradient and Hessian
%
% Test Version with errors -- No guarantee of any kind is given!!!
%

n = length(xact);
fnew = f(xnew);
if o_up > 2 && o_up < 5
   U = rando(n); % cost of U is at most O(n) function evaluations
else
   U = eye(n); 
end
fval = zeros(n,2);
for i = 1:n
   fval(i,1) = f(xnew-dt*U(:,i));
   fval(i,2) = f(xnew+dt*U(:,i));
end

gnew = U*(fval(:,2)-fval(:,1))/(2*dt);

dx = xnew-xact; dxdx = sum(dx.^2);
if dxdx > dt^2*n && o_up < 4
                 % PSB update only for sufficiently long steps; else
                 % the finite difference error may falsify the update
   dg = gnew-g;
   ddg = dg-H*dx;
   dH = (ddg*dx'+dx*ddg')*(1.0/dxdx) - ((ddg'*dx/dxdx^2)*dx)*dx';
   dH = 0.5*(dH+dH');
   H = H + dH;
end

if o_up > 1 && o_up < 6
   betat = (sum(fval,2)-2*fnew*ones(n,1))/(dt^2);
   betai = diag(U'*H*U);
   dbeta = betat-betai;

   H = H + U*diag(dbeta)*U'; % curvature correction
end

if o_up > 5
   H = zeros(size(H));
end

Hnew = 0.5*(H+H');


end










function [y,xt] = curvilp(fin,t,x,V,D,ng,lb2,ub2,inact,ocl)
% compute a point on a projected trust region curve and evaluate f
% at this point (used for a curvilinear search where 1 >= t >= 0)
%
% function [y,xt] = curvilp(fin,t,x,V,D,ng,lb2,ub2,inact,ocl)
%
% Input:
%        a function handle f
%        a scalar t in [0,1], a vector x
%        a unitary matrix V, a diagonal matrix D
%        a direction ng ``negative gradient''
%        lower und upper bounds lb2 and ub2 for xt
%        inact a subset of 1:n on which x is being changed
%        ocl = 1 (output, CurviLp) means `evaluate f at xt'
%        ocl = 0  means `return y = Inf' and save the function evaluation
%
% When there is no finite bound then compute
%        f(x+dx)
% where
%        dx = V* (D+(ev_min+(1-t)/(tiny+t))*eye(n))^(-1) *V'*ng
% and    tiny = (1.0e-12)/(1+max(abs(d)))
% and    ev_min = -min(diag(D)) + (1+abs(min(diag(D))))*10*eps.
%
% NOTE: This allows for very long steps even if the estimate
%       H = D*D*V' of the Hessian is positive definite.
%       Given that the gradients are estimated with randomized bases,
%       the ``hard case'' (of More and Sorensen) is not considered here.
%
% In case of finite bounds keep the variables associated with active
% indices fixed and project the remaining variables onto the bounds.
%
% OUTPUT:
%   y        -- the associated function value (if ocl == 1)
%   xt       -- point along a projected curve
%
% Florian Jarre, Feb. 2015
%
% Test Version with errors -- No guarantee of any kind is given!!!
%

if t < 0 || t > 1
   error( ' negative step length tested ' )
end

f = @(x) fin(x);

n = sum(inact);
if n > 0
    e = ones(n,1);
    d = diag(D);

    ev_min = -min(d);
    ev_min = ev_min + (1+abs(ev_min)) * 10*eps;
    tiny = (1.0e-12)/(1+max(abs(d)));

    dtmp = d + (ev_min+(1-t)/(tiny+t))*e;
    xtmp = x(inact)+V*((V'*ng(inact))./dtmp);

    xt = x;
    xt(inact) = xtmp;                      % inactive indices fixed
    iactl = ~inact & ((x-lb2) <  (ub2-x)); % lower active indices
    iactu = ~inact & ((x-lb2) >= (ub2-x)); % upper active indices
    xt(iactl) = lb2(iactl);
    xt(iactu) = ub2(iactu);
    xt = max(min(xt,ub2),lb2); % also project the xtmp-variables
else % all indices are active (this should never happen in a line search)
    xt = x;
    iactl = ((x-lb2) <  (ub2-x));
    iactu = ((x-lb2) >= (ub2-x));
    xt(iactl) = lb2(iactl);
    xt(iactu) = ub2(iactu);
    xt = max(min(xt,ub2),lb2);
end
if ocl == 1
   y = f(xt);
else
   y = Inf;
end
end








function U = rando(n)
% generate a random orthogonal n x n - matrix U
if n <= 0
   error('check input for rando');
end

if n <= 1
   U = randi(2,1,1); % a random number either 1 or 2
   U = 2*U-3;        % a random number either 1 or -1
else
   %a = randn(n);a=a+a';[U,D]=eig(a);
   a = randn(n);[U,~]=qr(a);
end

end
