function [ x,fx,g,H,out ] = mwd_bnd( fin,x0,options )
%
% Minimization Without using Derivatives, with BouNDs on the variables.
%
% Given a smooth function f=fin: R^n --> R  (where n >=1)
% approximate a local optimal solution of
%
%    minimize  f(x)  for  lb <= x <= ub 
%
% where f is not necessarily defined for x \not\in [lb,ub].
% Missing values of lb (or ub) are interpreted as -inf (or inf).
%
% The approximate solution generally is accurate to about 6 digits, but 
% in some cases, the final error is considerably larger; it may also 
% happen that the output is far away from any local optimal solution.
%
% Florian Jarre, July 2015
%
% Test Version with errors - you are free to use it; in return please
% report bugs to: jarre@hhu.de
%
% calling routine, for example:
%
%   [x,y]         = mwd_bnd(@sin) % no bounds, just one variable
% or
%   [x,y,g,H,out] = mwd_bnd(@f_rosen, [-1;1;1]) % no bounds, 3 variables
%
% Mandatory input: 
%   A function handle fin,
%   and, for n > 1, also a starting point x0 in R^n is mandatory.
%   x0 does not have to satisfy the bounds,
%   but if it does, the returned value shall be at least as good as f(x0).

%   (When x0 is not provided it is assumed that the input for f is a
%    single (``scalar'') number.)
%
% 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.tol   -- Only relevant for n = 1: 
%                    An approximate tolerance for the minimizer x:
%                    In the nonsmooth case there is a local minimizer 
%                    xm of f satisfying |xm - x| <= tol * int_length where 
%                    int_length is the length of a sub interval of [x0,x1] 
%                    generated by the algorithm. (Default tol = 1e-8)
%                    In the smooth case the above stopping criterion is
%                    based on an estimate - and is not guaranteed!
%                    WARNING: 
%                    Direct search for smooth minimization usually reaches 
%                    at most half the digits of full machine precision
%                    (For n > 1 the stopping criterion is too unreliable
%                    to allow any specification; use maxit instead to limit  
%                    the computational effort.)
%
%   options.smth  -- Only relevant when n=1:
%                    Binary variable indicating whether f is considered
%                    a smooth function. (If yes, spline interpolation
%                    is applied, if not plain golden mean search is used.)
%                    Default is 1.
%   options.update - update gradient by central differences and Hessian by:
%                    1:  PSB only (default)
%                    2:  PSB and curvature correction, unit basis
%                    3:  PSB and curvature correction, randomized bases
%                    4:  curvature correction only, randomized bases
%                    5:  curvature correction only, unit basis
%                    6:  zero Hessian (steepest descent)
%                                (4,5, and 6 are for test purposes only)
%
%
% Output:
%   x:         an approximate local minimizer, 
%   fx:        the associated function value, 
%   g:         the final gradient approximation (unreliable)
%              When n=1 then g = NaN unless the spline interpolation used
%              for minimizing f is deemed reliable; then g approx. f'(x),
%   H:         the final Hessian approximation (very unreliable)
%              When n=1 then H = NaN unless the spline interpolation used
%              for minimizing f is deemed reliable; then H approx. f''(x),
%   out.fval:  number of function evaluations needed,
%   out.iter:  number of iterations needed, (relevant for n>1)
%   out.fline: average number of function eval's in line search. (for n>1)
% 
%
% Algorithm:
% For n=1 use golden mean search and acceleration by spline interpolation. 
% For n>1:
% Trust region algorithm with finite central difference approximation of  
% the gradient using a PSB update for the Hessian.


if (nargin < 2)
   x0 = 0; % no starting point, assume n = 1. 
end
n = length(x0);
   
if n == 1 
   options.xact = x0;
   [x,fx,g,H,out] = mwd11(fin,options);
   out.iter = out.it;
   clear out.it;
   out.fline = out.iter;
   out.fval = out.iter; % for n=1 this is the only meaningful criterion
else
   if (nargin < 3)
      [x,fx,g,H,out] = mwd(fin,x0);
   else
      [x,fx,g,H,out] = mwd(fin,x0,options);
   end
end

end





















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 with n >= 2
%
% minimize  f(x)  for  lb <= x <= ub 
%
% f is not necessarily defined for x \not\in [lb,ub]
%
% last change: July 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
% 
% subroutines used:
%    updategh    - Update gradient and Hessian
%    curvilp     - Evaluate f at a projected point along a certain curve
%    rando       - generate a random orthogonal matrix
%    mwd11       - Line search routine 
%

%RandStream.setDefaultStream(RandStream('mt19937ar','seed',100)); %Matlab
%RandStream.setGlobalStream(RandStream('mt19937ar','seed',100)); %Matlab
%v = 1:625; rand ("state", v); %octave

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^(.65); %XXX Sept.2015: eps^(.5) --> eps^(.65)
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.lb   = 0;         % lower bound for curvi-linear search
optL.ub   = 1;         % upper bound for curvi-linear search
optL.err  = err;
optL.tol  = 1.0e-6;
converged = 0;         % stopping the loop
outerit = 1;

slowsteps = 0;



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] = ...
      mwd11(@(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
   
   if steplength < 10*dt
       slowsteps = slowsteps+1;
   else
       slowsteps = 0;
   end
   if slowsteps >= 10
       converged = 1; % 10 consecutive steps little progress
   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] = ...
    mwd11(@(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
%
% last change: 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
%
% last change: 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











function [x,y,g,H,out] = mwd11(f,options)
% Minimization of a continuous function f: R --> R union {Inf}
%
% last change: July 2015.
%
% Simplest calling routine
%
%   x = mwd11(@f);
%
% or, with more specifications:
%
%   [x,fx,out] = mwd11(@f,options);
%
% Test Version with errors -- No guarantee of any kind is given!!!
%
% Algorithms: Some form of bisection and spline interpolation.
%             Return the lowest point that the algorithm stumbles about
%             (in contrast to the official Matlab routine)
%
%
% INPUT: 
%   Mandatory: A function handle f,
%   Optional:  The structure options as outlined below.
%
%   options.lb   - The minimizer is restricted to the interval [lb,ub]. 
%                  (Default of x0 is -Inf)
%   options.ub   - (Default of x1 is  Inf)
%                  If x0 or x1 are specified f will be evaluated only
%                  within the given bounds
%   options.xact - Reference point for the line search. 
%                  When lb <= xact <= ub the returned value shall be at 
%                  least as good as f(xact).
%
%   options.tol  - An approximate tolerance for the minimizer x:
%                  In the nonsmooth case there exists a local minimizer xm
%                  of f satisfying |xm - x| <= tol * int_length where 
%                  int_length is the length of a sub interval of [x0,x1] 
%                  generated by the algorithm. (Default tol = 1e-8)
%                  In the smooth case the above stopping criterion is
%                  based on an estimate - and is not guaranteed!
%                  WARNING: 
%                  Direct search for smooth minimization can never get more
%                  than about half the digits of full machine precision
%
%   options.smth - Binary variable indicating whether f can be considered
%                  to be a smooth function. (If yes, spline interpolation
%                  is applied, if not plain golden mean search is used.)
%                  Default is 1
%
%
% OUTPUT:
%   x        -- some approximate minimizer
%   y        -- the associated function value
%   g        -- NaN or an estimate of f'  at x, if available from
%               spline interpolation
%   H        -- NaN or an estimate of f'' at x, if available from
%               spline interpolation
%   out.it   -- specifies the number of function evaluations needed.
%
%
% Subroutines used:
%           find_spline.m -- find a least squares spline through xx and ff
%           eval_spline.m -- evaluate the spline at a given point
%           min_spline.m  -- find the minimizer of the spline

done       = 0;    % We are not done yet
spline_int = 0;    % so far no use of spline interpolation

% COMPLETE OPTIONAL INPUT AND GENERATE STARTING POINT xact in (x0,x1):

if (nargin < 2)
   options.lb   = -Inf;
   options.ub   =  Inf;
   options.xact = 0.0;
   options.tol  = 1.0e-8;
   options.smth = 1;
end

% Set x0 and x1: 
if ~isfield(options,'lb') 
   options.lb = -Inf; 
end

if ~isfield(options,'ub')
   options.ub =  Inf;
end

x0 = options.lb;
x1 = options.ub;
if x1 < x0
   error('lower bound larger than upper bound in line search');
end

iterations = 0;

if x0 == x1
   if x0 == -Inf || x0 == Inf
      error('bounds in line search are inconsistent');
   end
   options.xact = x0;
   xact         = x0;
   fact         = f(xact); f0 = fact; f1 = fact;
   iterations   = iterations + 1;
   done         = 1;
   warning('Trivial line search in interval of length zero');
end

% Set xact:
if ~isfield(options,'xact') 
   options.xact = x0;
end
xact = options.xact;
if xact <= x0 || xact >= x1 % Make sure xact is in the interior of [x0,x1]
   if x0  > -Inf
      if x1 < Inf
         xact = 0.5*(x0+x1);
      else
	     xact = x0 + 0.5*abs(x0) + 1;
      end
   else % Case x0 = -Inf and then
      if x1 < Inf
         xact = x1 - 0.5*abs(x1) - 1;
      else
         xact = 0;
      end
   end
end
if done == 0;
   fact = f(xact);
   iterations = iterations + 1; % Number of function evaluations
end
% Now, x0 < xact < x1  but x0 = -Inf of x1 = Inf is possible


% Set the tolerance
if ~isfield(options,'tol') 
   options.tol = 1.0e-8;
end
tol = min(0.1,max(1.0e-8,options.tol)); % Do not allow very high precision 
                                        % or very low precision

if ~isfield(options,'smth') % Assume smooth minimization
   options.smth = 1;
end

                                        
% Complete function values for x0, x1
if done == 0
   f0 = Inf;
   f1 = Inf;
end
if done == 0 && x0 > -Inf
   f0 = f(x0); iterations = iterations+1;
end
if done == 0 && x1 < Inf
   f1 = f(x1); iterations = iterations+1;
end




% GENERATE FINTE BOUNDS

dt = 1 + abs(xact)*1.0E-3; 
% (for large xact an increment by one may be too small)
if x1 < Inf
   dt = max(dt,x1-xact);
end
if x0 > -Inf
   dt = max(dt,xact-x0);
end
dtsave = dt;


% Make one of the bounds finite
if x0 == -Inf && x1 == Inf % Here, it must be ``done = 0''
   xtmp = xact+dt;
   ftmp = f(xtmp); iterations = iterations+1;
   if ftmp < fact
      x0 = xact;
      f0 = fact;
      xact = xtmp;
      fact = ftmp;
   else
      x1 = xtmp;
      f1 = ftmp;
   end
end


% Generate a finite upper bound
if x1 == Inf && f0 >= fact % Here, it must be ``done = 0''
   xtmp = xact + dt;
   ftmp = f(xtmp); iterations = iterations+1;
   itcount = 0;
   while itcount < 15 && ftmp < fact % Line search up to length 10^15
      x0 = xact;
      f0 = fact;
      xact = xtmp;
      fact = ftmp;
      dt = dt*10;
      xtmp = xact + dt;
      ftmp = f(xtmp); iterations = iterations+1; itcount = itcount+1;
   end
   if itcount >= 15
      done = 1;
      if ftmp < fact
         xact = xtmp;
         fact = ftmp;
      end
      out.it = iterations;
      warning('line search may be unbounded (x to Inf)');
   else % itcount < 15 means ftmp >= fact
      x1 = xtmp;
      f1 = ftmp;
   end
end

if x1 == Inf && f0 < fact % Here, it must be ``done = 0''
   x1 = xact;
   f1 = fact;
   xact = 0.5*(xact+x0);
   fact = f(xact); iterations = iterations +1;
end


% Generate a finite lower bound
dt = dtsave;
if x0 == -Inf && f1 >= fact && done == 0
   xtmp = xact - dt; 
   ftmp = f(xtmp); iterations = iterations+1;
   itcount = 0;
   while itcount < 15 && ftmp < fact % line search up to length 10^15
      x1 = xact;
      f1 = fact;
      xact = xtmp;
      fact = ftmp;
      dt = dt*10;
      xtmp = xact - dt;
      ftmp = f(xtmp); iterations = iterations+1; itcount = itcount+1;
   end
   if itcount >= 15
      done = 1;
      if ftmp < fact
         xact = xtmp;
         fact = ftmp;
      end
      out.it = iterations;
      warning('line search may be unbounded (x to -Inf)');
   else % itcount < 15 means ftmp >= fact
      x0 = xtmp;
      f0 = ftmp;
   end
end

if x0 == -Inf && f1 < fact && done == 0
   x0 = xact;
   f0 = fact;
   xact = 0.5*(xact+x0);
   fact = f(xact); iterations = iterations +1;
end
int_length = x1-x0; % Length of the interval containing the minimizer
int_length = max(int_length,10*eps*(abs(x0)+abs(x1))/tol);


% Eliminate Inf-values of f
if min([f0,fact,f1]) == Inf
   error('No finite value of f found in line search')
end
itcount = 0;
if fact == Inf
   if f1 < f0 % look for minimizer near x1
      while fact == Inf && itcount < 15
         itcount = itcount+1;
         x0 = xact; 
         f0 = fact;
         xact = 0.1*xact+0.9*x1;
         fact = f(xact); iterations = iterations+1;
      end
      if fact == Inf
         error('Failure in line search, only finite value may be x1')
      end 
   else % look for minimizer near x0
      while fact == Inf && itcount < 15
         itcount = itcount+1;
         x1 = xact; 
         f1 = fact;
         xact = 0.1*xact+0.9*x0;
         fact = f(xact); iterations = iterations+1;
      end
      if fact == Inf
         error('Failure in line search, only finite value may be x0')
      end 
   end
end % Now, fact is finite
if f1 == Inf
   while f1 == Inf && itcount < 30
      x1old = x1;
      x1 = 0.5*(xact+x1);
      f1 = f(x1); iterations = iterations+1;
   end
   if f1 == Inf
      error('f may be discontinuous at xact')
   else
      if f1 < min(f0,fact)
         x0 = xact;
         f0 = fact;
         xact = x1;
         fact = f1;
         x1 = x1old;
         f1 = Inf;
      end
   end
end
if f0 == Inf
   while f0 == Inf && itcount < 30
      x0old = x0;
      x0 = 0.5*(xact+x0);
      f0 = f(x0); iterations = iterations+1;
   end
   if f0 == Inf
      error('f may be discontinuous at xact')
   else
      if f0 < min(f1,fact)
         x1 = xact;
         f1 = fact;
         xact = x0;
         fact = f0;
         x0 = x0old;
         f0 = Inf;
      end
   end
end      
       





% MAKE SURE f0 AND f1 ARE AT LEAST AS LARGE AS fact 

itcountmax = 9+round(-log(tol)/log(10)); % higher precision near end points
% Number of iterations to identify a minimizer near the boundary

if fact > min(f0,f1) && done == 0
   itcount = 0;
   if f0 < f1
      while fact >= f0 && itcount < itcountmax
         x1 = xact; f1 = fact; itcount = itcount+1;
         xact = 0.9*x0+0.1*x1; fact  = f(xact); 
      end
   else
      while fact >= f1 && itcount < itcountmax
         x0 = xact; f0 = fact; itcount = itcount+1;
         xact = 0.1*x0+0.9*x1; fact  = f(xact); 
      end
   end
   iterations = iterations + itcount;
   if itcount >= itcountmax
      done = 1;
      if f0 < min(f1,fact)
         xact = x0;
         fact = f0;
      end
      if f1 < min(f0,fact)
         xact = x1;
         fact = f1;
      end
      out.it = iterations;
   end
end
% Now, either ``x1-x0 <= tol*int_length'' or ``fact < min(f0,f1)''
if xact <= x0 || x1 <= xact || f0 < fact || f1 < fact
   if done == 0
      error( ' programming error in line search 1' );
   end
end






% NOW, THE ACTUAL LINE SEARCH APPROXIMATING A MINIMIZER IN [x0,x1]

itcount = 0;                      % iteration counter
gm6 = 2/(sqrt(5)+1);              % Golden mean ratio, this is about 0.6
xfval = [x0,xact,x1;f0,fact,f1];  % Record all values of the search
iact = 2;                         % Index of xact in xfval
    
          
   
while done == 0 % *** MAIN LOOP ***
   spline_int = 0;
   itcount = itcount+1;

   % One golden mean search step
   if x1 - xact > xact - x0
      xa = x1 + gm6*(xact-x1); 
      fa = f(xa); iterations  = iterations + 1;  
      xfval = [xfval(:,1:iact),[xa;fa],xfval(:,iact+1:end)];
      if fa <= fact
         x0 = xact; xact = xa; fact = fa; iact = iact+1;
      else
         x1 = xa; 
      end      
   else
      xa = x0 + gm6*(xact-x0); 
      fa = f(xa); iterations  = iterations + 1;
      xfval = [xfval(:,1:iact-1),[xa;fa],xfval(:,iact:end)];
      if fa <= fact
         x1 = xact; xact = xa; fact = fa;
      else
         x0 = xa; iact = iact+1;
      end
   end % of golden mean step     
   if x1 - x0 <= tol*int_length
       done = 1;
   end
      
   
   % Test whether to use spline interpolation
   n = length(xfval);
   if n >=5 && options.smth == 1 && done == 0
      % Check whether the spline would have predicted fact correctly
      % Find points close to xact on both sides (2 <= iact <= end-1)
      if n == 5
         indspl = [1:iact-1,iact+1:5];
      else
         i_set = 0; % indspl not yet set
         if iact <= 2
            if iact <= 1
                error( 'programming error in line search')
            end
            indspl = [1,3,4,5,6]; 
            i_set = 1; % do not change indspl any more 
         end
         n = length(xfval);
         if iact >= n-1
            if iact >= n
                error( 'programming error in line search')
            end
            indspl = [n-5,n-4,n-3,n-2,n];
            i_set = 1; % do not change indspl any more
         end
         if i_set == 0 % now, 3 <= iact <= n-2
            indspl0 = [iact-2,iact-1,iact+1,iact+2]; 
            tmp     = [-Inf,xfval(1,:),Inf]; % Note: tmp(iact+1)=xact
            if tmp(iact+1)-tmp(iact-2) < tmp(iact+4)-tmp(iact+1)
               indspl = [iact-3,indspl0]; 
            else
               indspl = [indspl0,iact+3]; 
            end
         end
      end % Index for spline is set
      xx = xfval(1,indspl);
      ff = xfval(2,indspl);
      ss = find_spline(xx,ff); % Obiges koennte man direkt einsetzen!
      [sact,i] = eval_spline(xact,ss,xx);
      ip1 = min(1+1,length(xx)-1);
      if abs(sact-fact) < 0.1*(abs(ss(3,i))+abs(ss(3,ip1)))*...
                          (xact-xx(i))*(xx(i+1)-xact)
          spline_int = 1; % Estimate of second derivative is 80% correct
      end
   end
      
      
   % Possible spline acceleration (several steps)
   while spline_int == 1 && itcount < 100 && done == 0
      itcount = itcount + 1;
         
      % Recompute the spline including xact
      indspl = iact-2:iact+2; 
      n = length(xfval);
      if iact <= 1 || iact >= n
         error( 'programming error in line search')
      end
      if iact == 2
         indspl = 1:5;
      end
      if iact == n-1
         indspl = n-4:n;
      end
      xx = xfval(1,indspl);
      ff = xfval(2,indspl);
      ss = find_spline(xx,ff); % Obiges koennte man direkt einsetzen!
          
      [t,tval,i] = min_spline(ss,xx);
      if t <= x0 || t >= x1 % x0 = xfval(1,iact-1),  x1 = xfval(1,iact+1)
         spline_int = 0;
      else %%% Spline interpolation may be o.k. %%%
             
         ixf = i+indspl(1)-1; % t in [xfval(1,ixf),xfval(1,ixf+1)]
         if t < xx(i) || t > xx(i+1) || t<xfval(1,ixf) || t>xfval(1,ixf+1)
            error(' programming error in spline-line search ')
         end
         shift_t = 0.1*min(tol*int_length, xx(i+1)-xx(i));
         % less than tol*int_length even in the presence of rounding errors 
         if t < xx(i)+shift_t 
            t = xx(i)+shift_t;
         end
         if t > xx(i+1)-shift_t
            t = xx(i+1)-shift_t;
         end
         ft = f(t); iterations = iterations+1;
         epp = abs(tval-ft); % The approximation error at t
         ip1 = min(1+1,length(xx)-1);
         if epp > 0.2*(abs(ss(3,i))+abs(ss(3,ip1)))*(t-xx(i))*(xx(i+1)-t)
             spline_int = 0; % Prediction not so accurate
             %disp('return to golden mean search')
         end
         
         if epp < 100*eps*(abs(ft)+1)
            done = 1; % Prediction is close to machine precision
            spline_int = 1; % keep spline interpolation for derivatives
         end
         epp = 2*epp/((t-xx(i))*(xx(i+1)-t)); % Second der. of the error
         mss = ((xx(i+1)-t)*ss(3,i)+(t-xx(i))*ss(3,ip1))/(xx(i+1)-xx(i));
         if mss > epp
            tmp = max(0.1*(xx(i+1)-xx(i)), abs(0.5*(xx(i+1)+xx(i))-t) );
            if tmp*epp/(mss-epp) < int_length*tol
               done = 1; % prediction of error in t is small
               spline_int = 1; % keep spline interpolation for derivatives
            end
         end
         xfval = [xfval(:,1:ixf),[t;ft],xfval(:,ixf+1:end)];
         if t < xact
            if ft > fact
               iact = iact+1; % xact remains same but iact increases by 1
            end
         else
            if ft < fact
               iact = iact + 1; % xact changes and also iact increases by 1
            end
         end
            
         x0   = xfval(1,iact-1);
         xact = xfval(1,iact  ); 
         fact = xfval(2,iact  );
         x1   = xfval(1,iact+1);
         if x1-x0 <= tol*int_length
            done = 1;
            spline_int = 1; % keep spline interpolation for derivatives
         end
      end %%% Case where spline interpolation may be o.k. %%%
   end   
   if x1-x0 <= tol*int_length || itcount >= 100
      done = 1;
   end 
end % *** OF MAIN LOOP ***
x = xact; y = fact; out.it = iterations;


if spline_int == 1
   n = length(xfval);
   if xact < xfval(1,1) || xact > xfval(1,n)
      warning(' final point out of range in spline evaluation ');
      g = NaN; % No derivative information available from Spline
      H = NaN; % (Use finite difference instead - not done here)
   else
      i = 1;
      while xact > xx(1,i+1)
         i = i+1;
      end
      t = xact-xx(1,i);
      
      g = ss(2,i)+t*(2*ss(3,i)+t*(3*ss(4,i)));
      H = 2*ss(3,i)+t*(6*ss(4,i));
   end
else
   g = NaN; % No derivative information available from Spline
   H = NaN; % (Use finite difference instead - not done here)
end



end






function s = find_spline(x,f)
% Find the least squares spline interpolant through the points x(i), f(i)
% assuming that there are at least 4 interpolation points.
% On each interval [x(j),x(j+1)] (1<=j<=n-1) the spline is represented as
% s(x) = s(1,j) + s(2,j)*(x-x(j)) + s(3,j)*(x-x(j))^2 + s(4,j)*(x-x(j))^3

% Initialize
[n,m] = size(x);
if n < m % make x a column vector
   x = x'; tmp = n; n = m; m = tmp;
end
if m > 1 || m == 0 || n < 4
  error('input x for spline interpolation is inconsistent');
end

[nn,mm] = size(f);
if nn < mm % make f a column vector
   f = f'; tmp = nn; nn = mm; mm = tmp;
end
if mm > 1 || mm == 0
  error('input f for spline interpolation is inconsistent');
end

if n ~= nn
  error('input x,f for spline interpolation is inconsistent');
end

dx = x(2:n)-x(1:n-1); % vector of increments of x
if min(dx) <= 0
error('vector of spline base points is assumed to be in increasing order');
end

s  = zeros(12,n); % Three splines with n-1 cubic parts each.
% The last column is a dummy to allow assigning f as the first row.
% The first spline interpolates f,
% the second and third interpolate the zero function.
% rows 1 to 4  for the spline interpolating f
% rows 5 to 8  for the spline interpolating zero; nonzero linear first term
% rows 9 to 12 for the spline interpolating zero; nonzero quadr. first term
t1 = [1,5, 9]; % indices for the constant terms
t2 = [2,6,10]; % indices for the linear terms
t3 = [3,7,11]; % indices for the quadratic terms
t4 = [4,8,12]; % indices for the cubic terms

s(1,:) = f; % the constant part of s(1:4,:) is given by f
            % (constant parts of the zero splines are zero)


% Set up the three splines; first the splines on [x(1),x(2)]
s(2,1)  = (f(2)-f(1))/dx(1); % first spline, linear on [x(1),x(2)]
s(6,1)  = 1;                 % second spline, first linear term is 1
s(8,1)  = -1/dx(1)^2;        % interoplate to zero at x(2)
s(11,1) = 1;                 % third spline, first quadratic term is 1
s(12,1) = -1/dx(1);          % interoplate to zero at x(2)

% then the interpolating splines on the remaining subintervals
for i = 2:n-1
   s(t2,i) = s(t2,i-1)+dx(i-1)*(2*s(t3,i-1)+3*dx(i-1)*s(t4,i-1));
   s(t3,i) = s(t3,i-1)+3*dx(i-1)*s(t4,i-1);
   s(t4,i) = (s(t1,i+1)-s(t1,i)-dx(i)*(s(t2,i)+dx(i)*s(t3,i)))/dx(i)^3;
end
s = s(:,1:n-1); % now remove the last column

% The norm is given by || D * \Delta * s(4,:)' ||_2 where, formally,
% D      = Diag(((dx(2:n-1)+dx(1:n-2)).^-.5); and
% \Delta = [-eye(n-2),zeros(n-2,1)]+[zeros(n-2,1),eye(n-2)];
% i.e. a weighted sum of squared jumps of the third derivatives

d  = (dx(2:n-1)+dx(1:n-2)).^-.5; % weights
Ds = (kron(ones(3,1),d').*(s(t4,1:n-2)-s(t4,2:n-1)))';
% Ds contains the weighted jumps in the cubic terms of the splines

% orthogonalize the third with respect to the second spline
alpha     = Ds(:,3)'*Ds(:,2)/(Ds(:,2)'*Ds(:,2));
s(9:12,:) = s(9:12,:)-alpha*s(5:8,:); % updating s(10:12,:) would suffice
Ds(:,3)   = Ds(:,3)  -alpha*Ds(:,2);

% adjust first with the second spline
alpha     = Ds(:,1)'*Ds(:,2)/(Ds(:,2)'*Ds(:,2));
s(1:4,:)  = s(1:4,:)-alpha*s(5:8,:); % updating s(2:4,:) would suffice
Ds(:,1)   = Ds(:,1) -alpha*Ds(:,2);

% adjust first with the third spline
alpha     = Ds(:,1)'*Ds(:,3)/(Ds(:,3)'*Ds(:,3));
s(1:4,:)  = s(1:4,:)-alpha*s(9:12,:); % updating s(2:4,:) would suffice
%Ds(:,1)   = Ds(:,1) -alpha*Ds(:,3); % not needed

% to reduce rounding errors redo the final spline
s = s(1:4,:);
for i = 2:n-1
   s(2,i) = s(2,i-1)+dx(i-1)*(2*s(3,i-1)+3*dx(i-1)*s(4,i-1));
   s(3,i) = s(3,i-1)+3*dx(i-1)*s(4,i-1);
   s(4,i) = (f(i+1)-f(i)-dx(i)*(s(2,i)+dx(i)*s(3,i)))/dx(i)^3;
end
end






function [y,i] = eval_spline(xact,s,x)
% Evaluate the spline given by s and the partition x at the point xact
% Also return the segment in which xact is located
% Assume x in increasing order, length(x) = n and s is 1:4 by 1:n-1

n = length(x);

if xact < x(1) || xact > x(n)
    error(' point out of range in spline evaluation ');
end
i = 1;
while xact > x(i+1)
   i = i+1;
end

t = xact-x(i);
y = s(1,i)+t*(s(2,i)+t*(s(3,i)+t*(s(4,i))));
end







function [t,y,imin] = min_spline(s,x)
% Find the minimum of the spline given by s in the interval [x(1),x(end)]
% The minimizer t with value y in the interval [x(imin),x(imin+1)]
% Assume x in increasing order, length(x) = n and s is 1:4 by 1:n-1

n = length(x);

% Test support points first
[y,imin] = min(s(1,:)); t=x(imin);
dx = x(n)-x(n-1);
y_last = s(1,n-1)+dx*(s(2,n-1)+dx*(s(3,n-1)+dx*(s(4,n-1))));
if y_last < y
   y = y_last; t = x(n); imin = n-1; % Not imin = n, so x(imin+1) exists
end

for i = 1:n-1
   v = s(:,i); % just for convenience
   if v(4) ~= 0
      vv = v(1:3)/(3*v(4));
      if vv(3)^2 >= vv(2)
         root = sqrt(vv(3)^2-vv(2));
         t1 = -vv(3)-root;
         if t1 > 0 && t1 < x(i+1)-x(i)
            y1 = s(1,i)+t1*(s(2,i)+t1*(s(3,i)+t1*(s(4,i))));
            if y1 < y
               y = y1; t = x(i)+t1; imin = i;
            end
         end
         t1 = -vv(3)+root;
         if t1 > 0 && t1 < x(i+1)-x(i)
            y1 = s(1,i)+t1*(s(2,i)+t1*(s(3,i)+t1*(s(4,i))));
            if y1 < y
               y = y1; t = x(i)+t1; imin = i;
            end
         end
      end
   else
      if v(3) ~= 0
         t1 = -0.5*v(2)/v(3);
         if t1 > 0 && t1 < x(i+1)-x(i)
            y1 = s(1,i)+t1*(s(2,i)+t1*(s(3,i)+t1*(s(4,i))));
            if y1 < y
               y = y1; t = x(i)+t1; imin = i;
            end
         end
      end
   end
end
end
