function [ x,fx,g,H,out ] = min_f( fin,x0,options )
% AIM:
% 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].
%    f may depend on additional input parameters defined in options.par_f.
% Missing vectors 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, last change: Feb. 9, 2018
% Thanks to Felix Lieder and Markus Lazar for numerous corrections
%
% Test Version with errors - 
% you are free to use it; in return please report bugs to: jarre@hhu.de
%
% calling routine, for example:
%
%   x             = min_f(@sin) % no bounds, just one variable
%
% For more than one variable, a starting point must be provided, e.g.
%
%   [x,y,g,H,out] = min_f(@f_rosen, [-1;1;1]) % no bounds, 3 variables
%
% or, with additional input parameters,
%
%   [x,y]         = min_f(@f_rosen, zeros(n,1), options) 
%
% Mandatory input: 
%   A function handle fin,
%   and, for n > 1, also a starting point x0 in R^n is mandatory.
%   x0 is used to determine n and 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.par_f    If fin depends on additional parameters (not subject 
%                    to optimization), then the field options.par_f is a 
%                    struct containing the input parameters.
%                    For example: function y = f(x,A)
%                                          y = 0.5*x'*A*x;
%                    to be called e.g. as: x = min_f(@f,ones(n,1),options)
%                    where e.g. options.par_f = diag(1:n);
%                    Default:   options.par_f is not provided
%
%   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 [lb,ub] 
%                    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.update - update gradient by central differences and Hessian by:
%                    1:  PSB 
%                    2:  damped BFGS (default)
%   options.pl     - print level (e.g. for additional output if the 
%                    desired accuracy is not reached), default 0
%   options.mss    - Maximum Number of Slow Steps before the algoritm stops
%                    Default 10
%
%
% 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.acc:   a very unreliable (!!!) estimate for the accuracy of the 
%              approximate local minimizer
%   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 Quasi Newton 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.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 y = not_NaN(x)
% For minimization replace NaN with Inf

if isnan(x)
   y = Inf;
else
   y = x;
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: August 2016
% 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 
%                            2:  damped BFGS
%   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).
%   options.pl     - print level (e.g. for additional output if the 
%                    desired accuracy is not reached), default 0
%   options.mss    - Maximum Number of Slow Steps before the algoritm stops
%                    Default 10
%
% Output:
%   x:         an approximate local minimizer 
%   fx:        the associated function value
%   g:         the final gradient approximation
%   H:         the final Hessian approximation
%   out.acc:   a very unreliable (!!!) estimate for the accuracy of the 
%              approximate local minimizer
%   out.iter:  number of iterations needed
%   out.fval:  number of function evaluations needed
%   out.fline: average number of function evaluations in line search
% 
% subroutines used:
%    updategh    - Update gradient and Hessian
%    curvilp     - Evaluate f at a projected point along a certain curve
%    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


% 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);
end
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 = 300*n; 
end
if ~isfield(options,'update')
   options.update = 2; 
end
options.update = round(max(1,min(options.update,2))); 
% eliminate update-options 3-6 of earlier versions of min_f
if ~isfield(options,'pl')
   options.pl = 0; 
end
if ~isfield(options,'mss')
   options.mss = 10; 
end
mss = options.mss;

if ~isfield(options,'par_f')
   f = @(x) fin(x);
else
   f = @(x) fin(x,options.par_f);
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

converged = 0; % stopping criterion for main loop (avoid the word `stop')
iter = 0;
dt = (1+norm(x0))*eps^(.65); 
dtt = 2*dt;

xact = max(min(x0,ub-dtt),lb+dtt);
if min(xact-lb)<0 || max(xact-ub) > 0
   disp('lb and ub are close together, consider replacing with fixed var')
end
xact = max(min(xact,ub),lb); 
% reproject in case that dtt is large and the above is no longer feasible
fact = feval(f,xact); iter = iter + 1;     

if not_NaN(fact) == Inf
   disp('function to be minimized in min_f not defined at starting point');
   converged = 1;
end
if fact == -Inf
   disp('function to be minimized in min_f is -Inf at starting point');
   converged = 1;
end
% NOTE: When converged == 1 the gradient/Hessian Evaluation below should 
%       be omitted (correction of the code postponed since this should not 
%       happen in the first place.)


% estimate gradient and the error of the function evaluation of f near x0 
n = length(xact);
U = eye(n); % to have a deterministic algorithm 
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

rescalef = 0;
%largeg   = max(100,10*sqrt(n));
largeg   = 1; %XXX
if norm(g) > largeg % consider the problem as badly scaled
   rescalef     = 1;
   scaleffactor = largeg/norm(g);
   err  = err*scaleffactor;
   %fact = fact*scaleffactor; % not used any more
   %g    = g*scaleffactor;    % not used any more
   if ~isfield(options,'par_f')
      f = @(x) scaleffactor*fin(x);
   else
      f = @(x) scaleffactor*fin(x,options.par_f);
   end
else
   scaleffactor = 1;
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
err = min(err,1.0e-6)*(1+norm(x0)); 
%XXX ad hoc choice - some smoothness needed
dt   = 0.1*err^(1/3);
ub2  = ub-1.05*dt; lb2 = lb+1.05*dt; % search within the smaller set to 
                                     % allow for central differences
if min(ub2-lb2) < 0
    ub2 = ub; lb2 = lb; % hope that function values are defined outside
end

xact = max(min(x0,ub2),lb2); 
fact = f(xact);

%redo the gradient and do the Hessian
if o_up == 1
   H = zeros(n);
end
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
   dbeta = (sum(fval,2)-2*fact*ones(n,1))/(dt^2);
   dbeta = max(dbeta,1.0e-4*norm(dbeta)/sqrt(n)); % new YYY
   H = diag(dbeta); % here, U = I (identity)
end
   
if not_NaN(norm(g)) == Inf
   disp( ' Finite difference stencil at initial point in min_f failed ')
   converged = 1;
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-10;   % Todo: adaptive tolerance depending on norm of 
                       %       the Newton step (if defined)
outerit = 1;

slowsteps = 0;


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

   outerit = outerit + 1; 
   ng = -g;
   
   % PROJECT THE HESSIAN
   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 % all variables are at lower or upper bounds 
      Hact = H;   % Hact as if all variables were inactive
      converged = 1;
      if options.pl > 0
         disp('stop since all bounds are active') 
      end
   end
   if not_NaN(norm(Hact,'fro')) == Inf
      Hact = eye(size(Hact));
   end
   
   % LINE SEARCH
   if converged == 0
      [V,D] = eig(Hact);  % unitary V and diagonal D so that V'*Hact*V = D
      pars = preppars(D,V,ng,inact); % prepare parameters for search step
      t0 =0.99*(pars.opd+pars.ev_min*pars.tiny)/(pars.opd-pars.ev_min);
      %above corresponds to 0.99 * the plain (Quasi-) Newton step
      t0 = max(0,min(t0,0.99)); 
      optL.xact = t0; 
   
      ocl = 1; % curvilp shall return the function value, too 
      [tmin,fnew,~,~,out1] = ...
            mwd11(@(t) curvilp(f,t,xact,V,lb2,ub2,pars,ocl),optL);
      if fnew >= fact
         converged = 1;
         if options.pl > 0
           disp('stop since the line search did not yield any improvement') 
         end
      end
      if tmin == 0
         converged = 1;
         % this case should be redundant since then fnew >= fact
         if options.pl > 0
            disp('stop since the line search retruned step length zero') 
         end
      end
      if fnew == -Inf
        converged = 1;
         if options.pl > 0
            disp('stop since a point with function value -Inf was found') 
         end
      end
      
      iter = iter + out1.iter;
   
      ocl = 0; % the function value is not needed here
      [~,xnew] = curvilp(f,tmin,xact,V,lb2,ub2,pars,ocl);
      
      if converged == 1
         xnew = xact;
         fnew = fact;
      end
   else
      xnew = xact;
      fnew = fact;
   end
   
   % UPDATE GRADIENT AND HESSIAN
   steplength = norm(xnew-xact); 
   if steplength == 0
      converged = 1;
   else
      if outerit == 2 
         if o_up == 2 % for BFGS, initial H must be positive definite
            lbd    = max(0.001,abs(g'*(xnew-xact))/norm(xnew-xact)^2);
            H = lbd*eye(n);
         end
      end
      [g,H] = updategh(f,xnew,xact,g,H,dt,o_up);
      iter  = iter+2*n+1;
   end
   if not_NaN(norm(g)) == Inf
      disp('Finite difference stencil at some iteration in min_f failed')
      converged = 1;
   end

   fold = fact; fact = fnew; 
   xact = xnew; %xold is not used
   
   % STOPPING TEST
   if converged == 0
      converged = steplength < 0.01*dt && ...
                  fold - fact < 0.01*dt^2/scaleffactor;
      if options.pl > 0 && converged == 1
         disp('stop since step length and function reduction were tiny') 
      end
   end
   
   if steplength < 10*dt
       slowsteps = slowsteps+1;
   else
       slowsteps = 0;
   end
   if slowsteps >= mss 
      converged = 1; % mss consecutive steps little progress
      if options.pl > 0 
       fprintf('stop since %0.5g consecutive steps little progress\n',mss); 
         fprintf('dt is : %0.5g \n',dt);
         fprintf('last step length is : %0.5g \n',steplength);
      end
   end
      
end % OF MAIN LOOP



% DO A FINAL  STEP WITH ORIGINAL BOUNDS
if fact > -Inf && fact < Inf && norm(g) < Inf
   ng = -g;
   iact = (xact-lb<0.9*dt & ng<0) + (ub-xact<0.9*dt & ng>0); % act. Indices
   inact = ~iact;                                        % inactive indices
   if max(inact)
      Hact = H(inact,inact);
   else
      Hact = H;
   end

   if not_NaN(norm(Hact,'fro')) == Inf
      Hact = eye(size(Hact));
   end
   [V,D] = eig(Hact);   % unitary V and diagonal D so  that V'*H*V = D
   pars = preppars(D,V,ng,inact); %parameters used for search step
   ocl = 1;
   [tmin,fnew,~,~,out1] = ...
       mwd11(@(t) curvilp(f,t,xact,V,lb,ub,pars,ocl),optL);
   iter      = iter + out1.iter;
   ocl = 0;
   [~,x]     = curvilp(f,tmin,xact,V,lb,ub,pars,ocl);
   if fnew >= fact
      x = xact;
      fnew = fact;
   end
   fx        = fnew; 
   g         = g+H*(x-xact); % estimate the final gradient
else
   x  = xact;
   fx = fact;
   inact = 1:n;
end

if rescalef == 1
   fx = fx/scaleffactor;
   g  = g /scaleffactor;
   H  = H /scaleffactor;
end

out.iter  = outerit;
out.fval  = iter;
out.fline = iter/outerit-2*n-1;   % 2*n+1 function evaluations for gradient
out.acc   = not_NaN(norm(g(inact))); % assuming that inactive indices did 
                                     % not change in the last step
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
%
%
% Use finite differences to update gradient and PSB/BFGS for the Hessian
%
% Test Version with errors -- No guarantee of any kind is given!!!
%

n = length(xact);
%fnew = f(xnew);

U = eye(n); 
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 > 0.01*dt^2*n && o_up == 1
                 % 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 dxdx > 0.01*dt^2*n && o_up == 2
                 % damped BFGS update (for sufficiently long steps)
   y     = gnew-g;
   s     = dx;
   ys    = y'*s;
   hs    = H*s;
   shs   = s'*hs; %XXX H must be positive definite
   if ys < 0.01*shs
      theta = 0.99*shs/(shs-ys);
      y = theta*y+(1-theta)*hs;
   end
   dH = y*(y'/ys) - hs*(hs'/shs);
   dH = 0.5*(dH+dH.');
   H = H + dH;
end


Hnew = 0.5*(H+H.');


end








function [y,xt] = curvilp(fin,t,x,V,lb2,ub2,pars,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-10)/(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.
%       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!!!

f = @(x) fin(x);
al = 1; % a plane search modifying al and t simultaneously does not 
        % seem to pay off

if t < 0 || t > 1
   xt = x;
   y = Inf;
else
   xt = xpdx(t,al,x,V,lb2,ub2,pars);

   if ocl == 1
      y = f(xt);
   else
      y = Inf;
   end
end
end








function xt = xpdx(t,al,x,V,lb2,ub2,pars)
% compute x+dx(t,alpha) = x+alpha*dx(t/alpha) with alpha in [0,1]

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

n = sum(pars.inact);
if n > 0
    e = ones(n,1);
    d = pars.d;
    inact = pars.inact;
    ev_min = pars.ev_min;
    tiny = pars.tiny;
    
    t = t/al;
    dtmp = d + (ev_min+pars.opd*(1-t)/(tiny+t))*e;
    xtmp = x(inact)+al*(V*((pars.Vtng)./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

end








function pars = preppars(D,V,ng,inact)
% compute parameters needed for the search step

n = sum(inact);
pars.inact = inact;

if n > 0
   d = diag(D);
   pars.d = d;
%   pars.opd = 1; 
   pars.opd = 0.5*(1+norm(d)); 
   
   ev_min = -min(d);
%   ev_min = max(0,-min(d)); % for pos.def. Hessian do not allow more 
                            % descent along small curvature than Newton
   pars.ev_min = ev_min + (1+abs(ev_min)) * 10.0e-10;
   pars.tiny = (1.0e-10)/(1+max(abs(d)));
   pars.Vtng = V'*ng(inact);
   
else % all indices are active (this should never happen in a line search)
   pars.ev_min = eps;
   pars.tiny = eps;
   pars.d = zeros(0,1);
end

end








function [x,y,g,H,out] = mwd11(fin,options)
% Minimization of a continuous function f: R --> R union {Inf}
%
% last change: August 2016.
%
% 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 lb is -Inf)
%   options.ub   - (Default of ub is  Inf)
%                  If lb or ub 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 [lb,ub] 
%                  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
%
%
%
% 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.iter -- specifies the number of function evaluations needed.
%   out.acc  -- estimated length of final interval containing x
%
%
% 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
iterations = 0;    % number of function evaluations

% COMPLETE OPTIONAL INPUT AND GENERATE STARTING POINT xact in (lb,ub):

if (nargin < 2)
   options.lb   = -Inf;
   options.ub   =  Inf;
   options.xact = 0.0;
   options.tol  = 1.0e-8;
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
if ~isfield(options,'par_f')
   f = @(x) not_NaN(fin(x));
else
   f = @(x) not_NaN(fin(x,options.par_f));
end


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 or x1 = Inf is possible


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


                                        
% 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
      if ftmp < fact
         xact = xtmp;
         fact = ftmp;
      end
      out.iter = iterations;
      warning('line search may be unbounded (x to Inf)');
      done = 1;
   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
      if ftmp < fact
         xact = xtmp;
         fact = ftmp;
      end
      out.iter = iterations;
      warning('line search may be unbounded (x to -Inf)');
      done = 1;
   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);
tol = max(tol, 10*eps*(1+abs(x0)+abs(x1))/int_length);


% Eliminate Inf-values of f
if min([f0,fact,f1]) == Inf
   warning('No finite value of f found in line search')
   done = 1;
end
itcount = 0;
if fact == Inf && done == 0
   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
         warning('Only infinite objective values found in line search')
         done = 1;
      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
         warning('Only infinite objective values found in line search')
         done = 1;
      end 
   end
end % Now, fact is finite
if f1 == Inf && done == 0
   while f1 == Inf && itcount < 30
      x1old = x1;
      x1 = 0.5*(xact+x1);
      f1 = f(x1); iterations = iterations+1;
   end
   if f1 == Inf
      warning('Only infinite objective values found in line search')
      done = 1;
   else
      if f1 < min(f0,fact)
         x0 = xact;
         f0 = fact;
         xact = x1;
         fact = f1;
         x1 = x1old;
         f1 = Inf;
      end
   end
end
if f0 == Inf && done == 0
   while f0 == Inf && itcount < 30
      x0old = x0;
      x0 = 0.5*(xact+x0);
      f0 = f(x0); iterations = iterations+1;
   end
   if f0 == Inf
      error('Only infinite objective values found in line search')
   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 LARGER THAN 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 || x0 == xact || x1 == xact
      done = 1;
      if f0 < min(f1,fact)
         xact = x0;
         fact = f0;
      end
      if f1 < min(f0,fact)
         xact = x1;
         fact = f1;
      end
      out.iter = 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
lref = 0;                         % last refinement not used so far
    
          
   
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; f0 = fact; xact = xa; fact = fa; iact = iact+1;
      else
         x1 = xa; f1 = fa;
      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; f1 = fact; xact = xa; fact = fa;
      else
         x0 = xa; f0 = fa; 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 && 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]; % "Spline" with 4 interpolat. points
      else
%     Strategy: Choose 5 points, (if possible) two larger than xact, two
%     smaller, and the last one the closest to xact (among the rest).
         i_set = 0; % indspl not yet set
         if iact <= 2
            if iact <= 1
                error( 'programming error in line search 2')
            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 3')
            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); 
      [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
      
      
   while spline_int == 1 && itcount < 100 && done == 0
   %          use minimizer of the spline function (possibly several steps)
      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 4')
      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); 
          
      [t,tval,i] = min_spline(ss,xx);  
      % SPLINE MINIMIZER t in [xx(i),xx(i+1)]
      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 5')
         end
         shift_t = 0.1*min(tol*int_length, xx(i+1)-xx(i));
         % the shift is (much) less than tol*int_length 
         t = max(t,xx(i  )+shift_t);
         t = min(t,xx(i+1)-shift_t);
         tnewton = xact-0.5*ss(2,3)/ss(3,3); % Newton step for spline
         safetygap = min(abs(t-tnewton),0.1*min(x1-xact,xact-x0));
         % If xact is close to x0 or x1 and t is close to xact on the other
         % side, convergence may be slow ==> move a bit away from xact
         if xact-x0 < 0.1*(x1-xact) && t > xact
            t = t+safetygap; 
         end 
         if x1-xact < 0.1*(xact-x0) && t < xact
            t = t-safetygap; 
         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);
         f0   = xfval(2,iact-1);
         xact = xfval(1,iact  ); 
         fact = xfval(2,iact  );
         x1   = xfval(1,iact+1);
         f1   = xfval(2,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. %%%
      dddx = xfval(1,2:end)-xfval(1,1:end-1);
      if min(dddx) < 10*eps*(1+abs(x0)+abs(x1))
         done = 1;
      end   
   end   
   if x1-x0 <= tol*int_length || itcount >= 100
      done = 1;
   end 
   dddx = xfval(1,2:end)-xfval(1,1:end-1);
   if min(dddx) < 10*eps*(1+abs(x0)+abs(x1))
      done = 1;
   end
   
   if done == 1 % do a (final?) check
      check = 0;
      acc = max(x1-xact,xact-x0);
      if acc > int_length*tol
         check = 1;
      end
      if fact >= min(f0,f1)
         check = 0;
      end
      if min(x1-xact,xact-x0) <= 0
         check = 0;
      end
      if check > 0
         slope0 = (f0-fact)/(x0-xact); 
         slope1 = (f1-fact)/(x1-xact); 
         t0 = 0.5*(x0+xact);
         t1 = 0.5*(x1+xact);
         topt = t0-slope0*(t1-t0)/(slope1-slope0);
         if abs(xact-topt) < int_length*tol
             check = 0;
             %acc = abs(xact-topt);
         end
      end
      if check > 0
         if x1-xact > xact-x0
            t = xact + 0.9*int_length*tol;
            ft = f(t); iterations = iterations+1;
            xfval = [xfval(:,1:iact),[t;ft],xfval(:,iact+1:end)];
            if ft < fact
               done = 0;
               lref = 1;
               iact = iact + 1; % xact changes and also iact increases by 1
               x0   = xact;
               f0   = fact;
               xact = t;
               fact = ft;
            else
               x1 = t;
               f1 = ft;
            end
         else % we have x1-xact <= xact-x0
            t = xact - 0.9*int_length*tol;
            ft = f(t); iterations = iterations+1;
            xfval = [xfval(:,1:iact-1),[t;ft],xfval(:,iact:end)];
            iact = iact+1;
            if ft < fact
               done = 0;
               lref = 1;
               iact = iact - 1; % xact changes and also iact increases by 1
               x1   = xact;
               f1   = fact;
               xact = t;
               fact = ft;
            else
               x0 = t;
               f0 = ft;
            end
         end
         acc = max(x1-xact,xact-x0);
         if done == 1 && acc > int_length*tol
         % repeat the above (automatically this is the other side)
            if x1-xact > xact-x0
               t = xact + 0.9*int_length*tol;
               ft = f(t); iterations = iterations+1;
               xfval = [xfval(:,1:iact),[t;ft],xfval(:,iact+1:end)];
               if ft < fact
                  done = 0;
                  lref = 1;
                  iact = iact + 1; % xact changes, also iact decreases by 1
                  x0   = xact;
                  f0   = fact;
                  xact = t;
                  fact = ft;
               else
                  x1 = t;
                  f1 = ft;
               end
            else % we have x1-xact <= xact-x0
               t = xact - 0.9*int_length*tol;
               ft = f(t); iterations = iterations+1;
               xfval = [xfval(:,1:iact-1),[t;ft],xfval(:,iact:end)];
               iact = iact+1;
               if ft < fact
                  done = 0;
                  lref = 1;
                  iact = iact - 1; % xact changes, also iact decreases by 1
                  x1   = xact;
                  f1   = fact;
                  xact = t;
                  fact = ft;
               else
                  x0 = t;
                  f0 = ft;
               end
            end
         end
      end
   end
end % *** OF MAIN LOOP ***
x = xact; y = fact; out.iter = iterations;
spline_int = max(spline_int, lref);


if spline_int == 1 && y > -Inf
   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
      % Recompute the spline including xact
      indspl = iact-2:iact+2; 
      if iact <= 1 || iact >= n
         error( 'programming error in line search 6')
      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); 
      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

acc = max(x1-xact,xact-x0);
if spline_int == 1
   if isnan(g)
      out.acc = acc;
   else
      out.acc = min(abs(g)/(eps+abs(H)),acc); 
   end
else
   out.acc = acc;
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

% y is the smallest value so far
for i = 1:n-1
   v = s(:,i); % just for convenience
   if v(4) ~= 0 
      discr = v(3)^2-3*v(4)*v(2);
      if discr >= 0
         tmp = v(3)+sign(v(3))*sqrt(discr);
         t1 = -tmp/(3*v(4));
         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 = -v(2)/tmp;
         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 % when this is instable then the minimizer is outside 
                   % the current interval
         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

