function [xmin,fmin,nmin] = myfminL(myf,a,b,L,nsamp)

% function [xmin,fmin,nmin] = myfminL(myf,a,b,L,nsamp)
%
% This function minimizes the function myf (of a single variable)
% over the interval [a,b]
% using a set of random starting points sent to Matlab's fmincon,
% making use of a Lipschitz constant for the function.
%
% As more information is accumulated about the function, certain
% intervals are excluded from consideration, reducing the 
% number of minimizations.
%
% Inputs:
%   myf    a "handle" for a Matlab function
%   a      left  endpoint of the interval
%   b      right endpoint of the interval
%   L      a Lipschitz constant for myf:
%            | myf(x)-myf(y)| <= L (x-y) for all x, y in [a,b].
%          (Note that the absolute value of a bound on the 1st derivative of
%           myf is a valid L.)
%   nsamp  the number of random samples to generate
%
%  Outputs:
%    xmin  the computed minimizer
%    fmin  the computed minimum myf(xmin)
%    nmin  the number of calls made to fmincon
%
%  The function could be improved by only generating random numbers in
%  intervals that are not already excluded, and by iterating until all 
%  intervals are excluded or else reporting what intervals have not  
%  been excluded.
%
%  myfminL  Dianne P. O'Leary 10/2006
%                             11/2008 update
%                             10/2009 update


% In the array named excluded_intervals, we maintain a list of 
% intervals that do not need to be searched for the minimizer.
% This list is initially empty, but for bookkeeping purposes,
% we set it to two intervals of 0 length.

excluded_intervals = [a,a;b,b];

% Initialize our best guess at the minimizer to be the left endpoint.

xmin = a;
fmin = feval(myf,a); % 10/2009

% Initialize an array of Nsamp random starting points.
% Initialize a counter for the number of minimizations performed.

x0 =  a + (b-a)*rand(nsamp,1);

nmin = 0;

% Consider each of the starting points.

for i=1:nsamp,  % 11/2008

   x = x0(i);

% If x is not in an excluded interval, ...

   t1 = x >= excluded_intervals(:,1);
   t2 = x <= excluded_intervals(:,2);

   if (sum(t1.*t2) == 0)

% ... use it as a starting point for minimization.
% The point x is in the interval beginning with 
% excluded_intervals(i1(1)-1,2) and ending with excluded_intervals(i1(1)  ,1).
% Search this interval for the minimizer.

     nmin = nmin + 1;
     i1 = find(t1);
     lo = excluded_intervals(i1(end),  2);
     up = excluded_intervals(i1(end)+1,1);
     [xopt,fopt] = fmincon(myf,x,[],[],[],[],lo,up); % 10/2009

% See if xopt is the best point seen yet.

     if (fopt < fmin)
        xmin = xopt;
        fmin = fopt;
     else

% See if any interval can be excluded.
% If the new interval is long enough, insert it in the 
% list of excluded intervals.

        deltax = (fopt-fmin)/L;

        if (deltax > eps*a)
           excluded_intervals = ...
             update_interval_list(max(xopt-deltax,a), ...
                                  min(xopt+deltax,b), ...
                                  excluded_intervals);  % 11/2008

        end
     end
   end
end

[nmin, xmin, fmin]