function  [x,histout,complete_history]=imfil(x0,f,budget,bounds,options)
% IMFIL v0.85
% Minimization of noisy functions subject to explicit bound
% constraints + hidden constraints
%
% function  [x,histout,complete_history]=imfil(x0,f,budget,bounds,options)
%
%
% Input: 
%       x0 = initial iterate.
%        f = objective function.
%            The calling sequence for f should be
%            [fout,ifail,icount]=f(x,h)
%
%             h is an optional argument, and should be used only if your
%             function is scale-aware, ie does something useful with 
%             the scale, such as tuning accuracy.
%
%             fout = f(x)
%
%             fail = 0 unless the evaluation of f fails, in which case
%                    set ifail=1 and fout = NaN
%
%             icount = number of calls to the expensive part of f
%                     f should be smart enough to not do anything
%                     if a bound constraint or a cheap hidden constraint
%                     is violated
%
%   budget = max cost, uses icount from f(x)
%            The iteration will terminate after the iteration that
%            exhausts the budget. This argument is required.
%            Suggestion: try 10*N^2. 
%
% WARNING: budget may not be tuned to track expensive part completely
%
% bounds   = [low, high] = N x 2 array of bound constraints on the variables. 
%            This is REQUIRED. If you want to solve an unconstrained
%            problem, you must do so without the scaling we do in
%            f_fixed, and using imfil_core is currently the only way 
%            to do that. The lower bounds are the first column;
%            the upper bounds are the second.
%
%            I plan let you solve unconstrained problems by setting
%            the bounds to Inf or -Inf. That will take some time while
%            I figure out how I want to do the scaling (or make you
%            do it).
%
%  options = options structure
%            Documentation on the way. In the meantime, RTFM or 
%            look at imfil_optset.m.
%
% Output: 
%        x = estimated minimizer
%
%  histout = iteration history, updated after each nonlinear iteration
%          = (number of iterations )x (N+5) array, the rows are
%             [fcount, fval, norm(sgrad), norm(step), iarm, xval]
%             fcount = cumulative function evals
%             fval = current function value
%             norm(sgrad) = current projected stencil grad norm
%             norm(step) = norm of last step
%             iarm=line searches within current iteration
%                 =-1 means first iterate at a new scale
%             xval = transpose of the current iteration
%
% complete_history = complete evaluation history
%
% complete_history is a structure with the data on every evaluation of f.
%
% complete_history.good_points has the successful points for columns.
% complete_history.good_values is a $M x N$ matrix  of the values at
%              the good points. M > 1 for least squares.
%              good_values=f(good_points)
%
% complete_history.failed_points has the unsuccessful points for columns.
%
% You may want to use this to build surrogates, decide to add new points to
% the stencil, or for troubleshooting. The complete history can take up a
% lot of room. I will only return it as output if you ask for it, and will
% not store it at all if you set the complete_history option to 'off'.
% If I don't store it at all, then your functions to add directions to the
% stencil can't use that data.
%
% C. T. Kelley, July 19, 2009.
% This code comes with no guarantee or warranty of any kind.
%

global fname qbounds dbounds imfil_fscale imfil_scale_aware imfil_noise_aware imfil_least_squares
fname=f;
qbounds=bounds;
dbounds=bounds(:,2)-bounds(:,1);
%
% No options? Use the defaults.
%
if nargin < 5
    options=imfil_optset;
end
n=length(x0);
%
% Read the options so f_fixed will know about the globals.
%
imfil_fscale      = imfil_optread('fscale',options);
imfil_scale_aware = imfil_optread('scale_aware',options);
imfil_noise_aware = imfil_optread('noise_aware',options);
imfil_least_squares = imfil_optread('least_squares',options);
imfil_least_squares = imfil_optread('least_squares',options);
imfil_complete_history= imfil_optread('complete_history',options);
%
% Scale the initial iterate to [0,1] before shipping to imfil_core.
% imfil_core uses 0 and 1 as the bounds for optimization. 
%
%fbounds=zeros(n,2); fbounds(:,2)=1;
z0=(x0 - qbounds(:,1))./dbounds;
%if nargout <=2 || imfil_complete_history == 0
%[z,histout]=imfil_core(z0,@f_fixed,budget,options,bounds);
%else
[z,histout,complete_history]=imfil_core(z0,@f_fixed,budget,options,bounds);
%end
%
% Unscale the results to original bounds.
%
x=(dbounds.*z)+qbounds(:,1);
db=diag(dbounds);
qb=diag(qbounds(:,1));
xout=histout(:,6:n+5); 
xlow=ones(size(xout))*qb;
histout(:,6:n+5)=xout*db+xlow;
xout=histout(:,6:n+5); 
%
% Unscale the complete history
%
if nargout == 3 & imfil_complete_history == 1
[mg,ng]= size(complete_history.good_points);
complete_history.good_values=imfil_fscale*complete_history.good_values;
if mg > 0
   clow=qb*ones(size(complete_history.good_points));
   complete_history.good_points=db*complete_history.good_points+clow;
end
[mf,nf]= size(complete_history.failed_points);
if mf > 0
   flow=qb*ones(size(complete_history.failed_points));
   complete_history.failed_points=db*complete_history.good_points+clow;
end
end
%
% Unscale the function and gradients.
%
imfil_fscale=imfil_fscale^(1 + imfil_least_squares);
histout(:,2)=histout(:,2)*imfil_fscale;
histout(:,3)=histout(:,3)*imfil_fscale;

function [fx,iff,icf,tol]=f_fixed(x,h)
% F_FIXED
% Creates an O(1) dummy function with unit bounds.
% So, (I hope) the gradient and stencil are reasonably scaled.
%
% I encode scale and noise awareness in this function, too. If your
% function is not noise aware, I pass a zero noise back to imfil_core.
% imfil_core will do the right thing after that.
%
global fname qbounds dbounds imfil_fscale imfil_scale_aware ...
       imfil_noise_aware imfil_least_squares
[mx,nx]=size(x);
z=x;
for ix=1:nx
   z(:,ix)=(dbounds.*x(:,ix))+qbounds(:,1);
end
func_type=imfil_noise_aware + 10*imfil_scale_aware;
switch func_type
   case 0
     [fx,iff,icf]=feval(fname,z);
     tol=0;
   case 10
     [fx,iff,icf]=feval(fname,z,h);
     tol=0;
   case 1
     [fx,iff,icf,tol]=feval(fname,z);
   case 11
     [fx,iff,icf,tol]=feval(fname,z,h);
end
[mf,nf]=size(fx);
if imfil_least_squares == 1
   val = fx(:,1)'*fx(:,1)/2;
   scale_base=sqrt(val);
else
   scale_base=max(abs(fx(1)));
end
if scale_base > 0
   if imfil_fscale < 0
      imfil_fscale=abs(imfil_fscale)*scale_base;
   end
else
   imfil_fscale=abs(imfil_fscale);
end
%
% Scale the function value. If you've set imfil_fscale < 0, this will
% also change the sign for the rest of the computation and scale with
% the size of the function at the initial iterate.
%
%if imfil_fscale < 0
%   imfil_fscale=abs(imfil_fscale)*abs(fx(1));
%end
fx=fx/imfil_fscale;
tol=tol/imfil_fscale;
function [x,histout,complete_history] = imfil_core(x0,f,budget,options,bounds)
% IFFCO_CORE v0.8
% Core code for imfil.m
% Minimization of f(x) subject to explict bound constraints
%
%
% function [x, histout,complete_history] ...
%                    = imfil_core(x0,f,budget,options,bounds)
%
% Bound constrained, parallel, implicit filtering code.
% 
% IMPLICIT FILTERING with SR1 and BFGS quasi-Newton methods
%
% Input: 
%       x0 = initial iterate
%        f = objective function,
%
%            the calling sequence for f should be 
%            [fout,ifail,icount]=f(x,h)
%
%            h is an optional argument for f, and should be used only if your
%            function is scale-aware, ie does something useful with 
%            the scale, such as tuning accuracy.
%
%   budget = upper limit on function evaluations. The optimization
%                 will terminate soon after the function evaluation counter
%                 exceeds the budget.
%
%  options = options structure
%                  Documentation on the way. In the meantime, RTFM.
%
%  bounds  = N x 2 array of bound constraints on the variables.
%                   This is now REQUIRED. If you want to solve an unconstrained
%                   problem (1) set the lower/upper bounds to -INF/INF and
%                   hack the scaling within f_fixed (inside imfil.m). I do
%                   not recommend this. It's better to think about your
%                   problem and come up with some bounds on the variables.
%
% Output: 
%        x = estimated minimizer
%  histout = iteration history, updated after each nonlinear iteration 
%          = N+five column array, the rows are
%            [fcount, fval, norm(sgrad), norm(step), iarm, xval]
%            fcount = cumulative function evals
%            fval = current function value
%            norm(sgrad) = current simplex grad norm
%            norm(step) = norm of last step 
%            iarm=line searches within current iteration
%                =-1 means first iterate at a new scale 
%             xval = transpose of the current iteration
%
% complete_history = complete evaluation history
%
% complete_history is a structure with the data on every evaluation of f.
% complete_history.good_points has the successful points for columns.
% complete_history.good_values is a $M x N$ matrix  of the values at
%              the good points. M > 1 for least squares.
%              good_values=f(good_points)
%
% complete_history.failed_points has the unsuccessful points for columns.
%
% You may want to use this to build surrogates, decide to add new points to
% the stencil, or for troubleshooting. The complete history can take up a
% lot of room. I will only return it as output if you ask for it, and will
% not store it at all if you set the complete_history option to 'off'.
% If I don't store it at all, then your functions to add directions to the
% stencil can't use that data.

%
% C. T. Kelley, July 19, 2009
% This code comes with no guarantee or warranty of any kind.
%
global imfil_fname scale_aware imfil_noise_aware imfil_fscale
%
imfil_fname=f;
fcount=0; 
%
% And now for the knobs. Implicit filtering has too many of these and they
% can make a difference. These are headed for the options structure.
%
% The ones we use are:
%
% termtol (default = .01)
%   if norm(proj(difference_grad)) < termtol*h  we terminate at the scale
%   with success
%
% nterm controls termination on stencil failure for centered diffs (defalt = 0)
%       = 0 to terminate on stencil failure before starting the line search
%       = 1 to ignore stencil failure
%
% iquit (default = 3)
%    After iquit consecutive line search or stencil 
%    failures, we terminate the iteration
%
% The knobs are headed for the options structure and should be there for
% V 0.9
%
% set the knobs
%
termtol=.01; iquit=3; nterm=0;
%
% set up the difference scales and options
%
n=length(x0);
%
% In imfil_core the bounds are 0 and 1 on all the variables.
% We use the real bounds to unscale the complete_history array.
%
obounds=zeros(n,2); obounds(:,2)=1;
%
% Get the options we need.
%
imfil_function_error = imfil_optread('function_error',options);
imfil_least_squares  = imfil_optread('least_squares',options);
imfil_maxit          = imfil_optread('maxit',options);
imfil_maxitarm       = imfil_optread('maxitarm',options);
imfil_noise_aware    = imfil_optread('noise_aware',options);
imfil_target         = imfil_optread('target',options);
random_stencil = imfil_optread('random_stencil',options);
scale_aware    = imfil_optread('scale_aware',options);
stencil_wins   = imfil_optread('stencil_wins',options);
verbose        = imfil_optread('verbose',options);
imfil_complete_history = imfil_optread('complete_history',options);
%if imfil_complete_history > 0
    complete_history=struct('good_points',[],'good_values',[],...
                            'failed_points',[]);
%end
%
% Initialize the iteration; create the stencil; set up the scales; ...
%
x=x0; xold=x0; n=length(x0); histout=[];
dscal=imfil_create_scales(options); nscal=length(dscal);
v=imfil_create_stencil(options,n);
xc=x0; hess=eye(n); ns=0; iquitc=0;
stop_now=0;
fval=imfil_target+1;
%
% Sweep through the scales.
%
while (ns < nscal & fcount <= budget & iquitc < iquit & stop_now == 0 ...
       & fval > imfil_target)
    ns=ns+1; h=dscal(ns); 
%
% itc = scale-level iteration counter
%
    itc=0; 
%
% Evaluate the function to test for instant termination and (if
% noise_aware is on) to get the estimate of the noise. Both the noise
% and the value of f  may vary as a function of the scale, so we 
% have to test it every time. 
%
% We can move this outside of the main loop if noise_aware and scale_aware
% are both off. 
%
    [funs,iff,icf,noise_val]=feval(f,x,h);
%
% The first call to f sorts out the scaling. So you can't scale the
% target before that call. 
%
if fcount == 0
    imfil_target=imfil_target/imfil_fscale;
end
%
    fval=f_to_vals(funs,imfil_least_squares);
if imfil_complete_history > 0
   complete_history = ...
            single_point_hist_update(complete_history,x,funs,iff);
end
    fcount=fcount+icf;
    if fval < imfil_target
       stop_now=1;
       icount=0;
    end
    stol=termtol*h; iarm=0; lok=1;
%
% If you've hit the target then stop_now = 1, and you'll skip most
% of the loop.
%
% Compute the stencil gradient to prepare for the quasi-Newton step.
%
    if stop_now ==0
%        vv=random_augment(v,random_stencil);
vv=imfil_augment_directions(x,v,h,options,bounds);
        [sgrad,fb,xb,icount,sflag,tvar,diff_hist,jac]...
              =stencil_diff(x,f,h*vv,funs,obounds,options,h,complete_history);
        pgrad=x - kk_proj(x-sgrad,obounds);
if imfil_complete_history > 0
   complete_history=many_point_hist_update(complete_history,diff_hist);
end
%
% Scale tvar by fscale.
%
tvar=tvar/imfil_fscale;
%
% If noise_aware = 1 and the scaled variation in f is < than the
% function's estimate of the noise, then I declare stencil failure!
%
        if noise_val > tvar
           sflag=0;
        end
%
% If function_error > 0, then I terminate the entire optimization
% when the scaled variation in f < function_error. I report this as
% stencil failure as well.
% 
       if imfil_function_error > tvar;
          stop_now=1;
          sflag=0;
       end
%
    end
%
% Keep score.
%
    fcount=fcount+icount; 
    histout=[histout', [fcount, fval, norm(pgrad,inf), 0, -1, x']']';
    if norm(pgrad,inf) < stol | (sflag+nterm)==0 | noise_val > tvar ...
       | stop_now==1
%
%   Declare convergence at this scale on stencil failure or tolerance match.
%
       gc=sgrad; 
       if (sflag+nterm) ~= 0 
           iquitc=iquitc+1; 
       else
           iquitc=0;
       end
%
   else
%
%   Take a few quasi-Newton iterates.
%
       iquitc=0;
       while itc < imfil_maxit*n & fval > imfil_target &...
         norm(pgrad,inf) >= stol & lok==1 & fcount < budget & sflag+nterm > 0
           itc=itc+1;
%
%    Compute the difference gradient.
%
       gc=sgrad; 
       if(itc > 1)
%          vv=random_augment(v,random_stencil);
           vv=imfil_augment_directions(x,v,h,options,bounds);
          [sgrad,fb,xb,icount,sflag,tvar,diff_hist,jac]...
               =stencil_diff(x,f,h*vv,funs,obounds,options,h,complete_history); 
           pgrad=x - kk_proj(x-sgrad,obounds);
          fcount=fcount+icount;
          if imfil_complete_history > 0
             complete_history = ...
                 many_point_hist_update(complete_history,diff_hist);
          end
       end
       if sflag ==0 iarm=0; end
%
%     Watch out for stencil failure!
%
       if sflag+nterm > 0
%
%     Update the iterate and the model Hessian.
%
if imfil_least_squares == 0
[hess, xp, fval, fct, iarm, diff_hist] = imfil_qn_update(f, x, fval, ...
    xc, sgrad, gc, hess, obounds, itc, h, options);
funs=fval;
else
[hess, xp, fval, funs, fct,iarm,diff_hist] ...
         = imfil_gauss_newton(f, x, funs, jac, xc, gc, itc, ...
               obounds, h, options);
end
if imfil_complete_history > 0
   complete_history=many_point_hist_update(complete_history,diff_hist);
end
xc=x; x=xp;
%
fcount=fcount+fct;
%
%
%     reduce scale and (?) rotate directions upon failure of line search
%
            if iarm >= imfil_maxitarm
               lok=0;
%
% If the line search fails, but there's a better point in the stencil,
% take that better point. Keep the scale fixed.
%
               if fval > fb
                  lok=1;
               end
               x=xb; fval=fb;
            end
          end
%
%      keep the records
%
if stencil_wins == 1
      if fval > fb
      x=xb; fval=fb;
      end
end
            stepn=norm(xold-x,inf); xold=x; 
    histout=[histout', [fcount, fval, norm(pgrad), stepn, iarm, x']']';
%
      end % end of nonlinear step
%
   end % end of sweep through the scale
%
% At the end of the quasi-Newton iterate, I make sure that the 
% quasi-Newton point is the best I've seen. If it's not, I fix it.
% So, 'stencil_wins' is half-way on.
%
% If you see a negative projected gradient norm in histout, that tells you
% that we took a non-quasi-Newton best point at that scale.
%
if fval > fb 
   x=xb; fval=fb;
   histout=[histout', [fcount, fval, -1, 0, 0, x']']';
end
if verbose == 1 
    [ns,fval,h,norm(pgrad),itc,fcount]
end
end % end of loop over the scales


%
% This function is a wrapper for the options in the function call.
% Hack this at your peril.
%
function [fout,varargout]=f_local(x,varargin)
global imfil_fname scale_aware imfil_noise_aware
%
nv=length(varargin);
if nv > 0
   h=varargin{1};
end
func_type=imfil_noise_aware + 10*scale_aware;
switch func_type
   case 0
     [fout,ifail,icount,tol]=feval(imfil_fname,x);
   case 10
     [fout,ifail,icount,tol]=feval(imfil_fname,x,h);
   case 1
     [fout,ifail,icount,tol]=feval(imfil_fname,x);
   case 11
     [fout,ifail,icount,tol]=feval(imfil_fname,x,h);
end
nout=max(nargout,1)-1;
switch nout
    case 0
    case 1
       if imfil_noise_aware == 1
         varargout{1}=tol;
       else
         error('problem in imfil_core:f_local:tol');
       end
    case 2
       varargout{1}=ifail;
       varargout{2}=icount;
    case 3
       varargout{1}=ifail;
       varargout{2}=icount;
       varargout{3}=tol;
    otherwise
       error('problem in imfil_core:f_local');
end
function vstencil = imfil_create_stencil(options,n)
% IMFIL_CREATE_STENCIL
% Builds the stencil for imfil.m. As we evolve this we will be 
% adding the ability to do random rotations for all or part of
% a stencil and all sorts of other stuff.
%
% function vstencil = imfil_create_stencil(options,n)
%
% Input: 
%        options = imfil.m options structure
%              n = dimension
%
% Output:      
%              v = stencil
%
%
% C. T. Kelley, July 14, 2009
% This code comes with no guarantee or warranty of any kind.
%

stencil=imfil_optread('stencil',options);
vstencil=imfil_optread('vstencil',options);
if nargin ~= 2
   error('imfil_create_stencil requires two arguments');
end
%
% Is vstencil really there?
%
[mv,nv]=size(vstencil);
if mv+nv > 0
   stencil=-1;
end
%
% Build the stencil.
%
switch stencil
   case -1; % Custom stencil is ok.

   case 0 % central difference stencil
     v = [eye(n),-eye(n)];

   case 1 % one sided stencil = ON THE WAY OUT!
     v = eye(n);

   case 2 % positive basis stencil
     v = [eye(n),-ones(n,1)/sqrt(n)];

   otherwise
     error(' illegal stencil in imfil_create_stencil');
end
%
% If vstencil is a custom job, take it. Otherwise use one of
% the internal choices.
%
if stencil ~= -1
   vstencil=v;
end

function dscal = imfil_create_scales(options)
% IMFIL_CREATE_SCALES
% Creates the scales for imfil.m. Nothing much here right now, but 
% custom scales ... are in the near future.
%
% function dscal = imfil_create_scales(options)
%
%
% C. T. Kelley, September 15, 2008
% This code comes with no guarantee or warranty of any kind.
%
custom_scales=imfil_optread('custom_scales',options);
mcs=length(custom_scales);
if mcs > 0
   dscal=custom_scales;
else
   scalestart = imfil_optread('scalestart',options);
   scaledepth= imfil_optread('scaledepth',options);
%
% Do some error checking. Warn or complain as needed.
%
   if scalestart > scaledepth
      error(' imfil_create_scales: error in scales, scalestart > scaledepth');
   end
   dscal=-(scalestart:scaledepth); dscal=2.^dscal;
end
function [sflag,best_value,best_point,tvar,diff_hist]=...
              collect_stencil_data(good_points,good_values,failed_points,...
                                   best_point_old,best_value_old,fc,options)
least_squares=imfil_optread('least_squares',options);
%
% Who's number one?
%
good_scalars=f_to_vals(good_values,least_squares);
[xbest_value,ibest]=min(good_scalars);
xbest_point=good_points(:,ibest);
if xbest_value < best_value_old
   best_value = xbest_value;
   best_point = xbest_point;
else
   best_value = best_value_old;
   best_point = best_point_old;
end
%
% Find the big loser.
%
worst_value=max(good_scalars);
%
% Assemble the history structure.
%
diff_hist=struct('good_points',good_points,...
'good_values',good_values,'failed_points',failed_points);
%
% What's the total variation?
%
tvar=worst_value-best_value;
%
% Stencil failure?
%
sflag=1;
if abs(fc-best_value) < 1.d-14
   sflag=0;
   jac=[];
   grad=[];
end

function [fct,x,fval,iarm,fres,diff_hist]=...
  armijo_explore(f, sdir, fold, xc, h,...
                 options, obounds)
% ARMIJO_EXPLORE
% Line search for imfil.m.
%
% C. T. Kelley, September 15, 2008
%
% This code comes with no guarantee or warranty of any kind.
%
% function [fct,x,fval,iarm]=...
%   armijo_explore(f, sdir, fold, xc, maxitarm,beta,h,...
%                  options)
%
% This is an internal function, which you are NOT TO HACK! Since
% you may hack it anyway, I will tell you want is going on. If you
% break something, may Alberich's curse be upon you!
%
% Inputs: f = objective function.
%         sdir = quasi-Newton search direction.
%
%         fold = current function value.
%
%         maxitarm = limit on number of step size reductions.
%
%         beta = stepsize reduction factor.
%
%         h = current scale.
%
%         options = imfil options structure.
%
%         obounds = Nx2 array of scaled bound constraints from imfil_core.
%
% Output: fct = cost of line search in function evaluations.
%
%         x = new point.
% 
%         fval = f(x).
%  
%         iarm = number of step length reductions.  
%
%         fres = residual for nonlinear least squares problems.
%
%    diff_hist = history data for this iteration.
%
% Read the options array to get the ones we care about in the line search.
%
imfil_parallel=imfil_optread('parallel',options);
imfil_least_squares=imfil_optread('least_squares',options);
imfil_verbose=imfil_optread('verbose',options);
beta = imfil_optread('armijo_reduction',options);
maxitarm = imfil_optread('maxitarm',options);
%
if imfil_parallel == 0
[fct,x,fval,iarm,aflag,fres,diff_hist]=...
  serial_armijo(f, sdir, fold, xc, h, obounds, options);
else
[fct,x,fval,iarm,aflag,fres,diff_hist]=...
  parallel_armijo(f, sdir, fold, xc, h, obounds, options);
end
%
if iarm == maxitarm & aflag == 1 & imfil_verbose == 1
       disp(' line search failure'); [iarm, h]
end
function [fct,x,fval,iarm,aflag,fres,diff_hist]=...
  parallel_armijo(f, sdir, fold, xc, h, obounds, options)
% PARALLEL_ARMIJO
% Parallel line search for imfil.m.
% Uses extra processors to explore larger stencils.
%
% function [fct,x,fval,iarm,aflag,fres]=...
%    parallel_armijo(f, sdir, fold, xc, h, obounds, options)
%   parallel_armijo(f, sdir, fold, xc, maxitarm,beta,h,...
%                  imfil_verbose, obounds,imfil_least_squares)
%
% Before you modify this (to use a trust region method, say), consider
% writing your own with this as a template. I'll make this very easy
% in a later version of the code. 
%
% Inputs:
%          f = objective function or residual vector, depending
%              on whether this is a least squares problem or not.
%       sdir = search direction.
%       fold = current value of the objective.
%       xc   = current point.
%   maxitarm = limit on number of step size reductions.
%       beta = reduction factor of step size. Currently beta = 1/2.
%          h = current scale.
%     obounds = Nx2 array of scaled bound constraints from imfil_core.
%     options = imfil options structure
% imfil_least_squares = least squares flag.
%
% Output:
%        fct = cost in function evaluations.
%          x = new point.
%       fval = objective function value at new point = fres'*fres/2
%              in the nonlinear least squares case.
%      iarm  = number of step size reductions.
%      aflag = 0 if the search finds a new point, = 1 if not.
%      fres  = residual for nonlinear least squares problems
%  diff_hist = history structure for this iteration
%
%
% C. T. Kelley, July 21, 2009
% This code comes with no guarantee or warranty of any kind.
%
imfil_least_squares=imfil_optread('least_squares',options);
imfil_limit_quasi_newton=imfil_optread('limit_quasi_newton',options);
imfil_verbose=imfil_optread('verbose',options);
beta = imfil_optread('armijo_reduction',options);
maxitarm = imfil_optread('maxitarm',options);
%
% Initialize the line search.
%
lambda=1;
n=length(xc);
fct=0;
aflag=1;
dd=sdir;
fres=[];
frest=fres;
diff_hist=struct('good_points',[],'good_values',[],'failed_points',[]);
%
if imfil_limit_quasi_newton == 1
%
% If the quasi-Newton step is much longer than the scale, shrink it.
%
% I am not convinced that this is the right thing to do, but it
% really helps most of the time.
%
 smax=10*min(h,1); if norm(dd) > smax dd=smax*dd/norm(dd); end
%
end
x=xc;
fval=fold;
%
% Evaluate all steplength choices at once.
%
number_steps=maxitarm+1;
ddm=zeros(n,number_steps);
for i=1:number_steps
   ddm(:,i)=xc-lambda*dd;
   lambda=beta*lambda;
   ddm(:,i)=kk_proj(ddm(:,i),obounds);
end
[fta,iflaga,ictra]=feval(f,ddm);
%
% Put the failed points in the history structure.
%
ibad=(iflaga(1:number_steps)==1);
if sum(ibad) > 0
   diff_hist.failed_points=ddm(:,ibad);
end
%
% Get the good points for the history structure.
%
igood=(iflaga(1:number_steps)==0);
if sum(igood) > 0
    diff_hist.good_points=ddm(:,igood);
    if imfil_least_squares == 1
       diff_hist.good_values=fta(:,igood);
    else
       diff_hist.good_values=fta(igood);
    end
end
%
if imfil_least_squares == 1
     frest=fta;
end
fta=f_to_vals(fta,imfil_least_squares);
fct=fct+sum(ictra);
ilose=(iflaga==1);
fta(ilose) = fold+abs(fold)*.0001+1.d-12;
%
% Traditional duplicates a serial Armijo search. This will take
% the longest possible step and is sometimes the right thing to do.
%
% Non-traditional will take the best point from all the
% ones sampled in the search direction if the full step fails.
%
% I am still playing with this, and will let you chose in the options
% command pretty soon. 
%
traditional = 0;
if traditional == 0
%
[ft,it]=min(fta); xt=ddm(:,it);
iarm=maxitarm;
if ft < fval 
    aflag=0; fval=ft;
    if imfil_least_squares == 1
       fres=frest(:,it); 
    end
    x=xt; iarm=min(it,maxitarm-1); 
end
%
else
%
iarm=-1;
%
while iarm < maxitarm & aflag==1 
    ft=fta(iarm+2);
    xt=ddm(:,iarm+2);
    if ft < fval & aflag==1; aflag=0; 
            if imfil_least_squares == 1
                  fres=frest(:,iarm+2); 
            end
       fval=ft; x=xt; end
    iarm=iarm+1;
end
end
if iarm == maxitarm & aflag == 1 & imfil_verbose == 1
       disp(' line search failure'); [iarm, h]
end
function [fct,x,fval,iarm,aflag,fres,diff_hist]=...
  serial_armijo(f, sdir, fold, xc, h, obounds, options);
% SERIAL_ARMIJO
% Serial line search for imfil.m.
%
% function [fct,x,fval,iarm,aflag,fres]=...
%   serial_armijo(f, sdir, fold, xc, h, obounds, options)
%
% This is a plain vanilla serial line search. Nothing to see here;
% move along.
%
% Before you modify this (to use a trust region method, say), consider
% writing your own with this as a template. I'll make this very easy
% in a later version of the code.
%
% Inputs:
%           f = objective function or residual vector, depending
%              on whether this is a least squares problem or not.
%        sdir = search direction.
%        fold = current value of the objective.
%        xc   = current point.
%           h = current scale.
%     obounds = Nx2 array of scaled bound constraints from imfil_core.
%     options = imfil options structure
%
% Output:
%        fct = cost in function evaluations.
%          x = new point.
%       fval = objective function value at new point = fres'*fres/2
%              in the nonlinear least squares case.
%      iarm  = number of step size reductions.
%      aflag = 0 if the search finds a new point, = 1 if not.
%      fres  = residual for nonlinear least squares problems
%  diff_hist = history structure for this iteration
%
%
% C. T. Kelley, September 15, 2008
% This code comes with no guarantee or warranty of any kind.
%
imfil_least_squares=imfil_optread('least_squares',options);
imfil_limit_quasi_newton=imfil_optread('limit_quasi_newton',options);
imfil_verbose=imfil_optread('verbose',options);
beta = imfil_optread('armijo_reduction',options);
maxitarm = imfil_optread('maxitarm',options);
%
% Initialize the line search.
%
lambda=1;
n=length(xc);
iarm=-1;
fct=0;
aflag=1;
dd=sdir;
fres=[];
frest=fres;
diff_hist=struct('good_points',[],'good_values',[],'failed_points',[]);
%
if imfil_limit_quasi_newton == 1
%
% If the quasi-Newton step is much longer than the scale, shrink it.
%
% I am not convinced that this is the right thing to do, but it
% really helps most of the time.
%
 smax=10*min(h,1); if norm(dd) > smax dd=smax*dd/norm(dd); end
%
%
end
x=xc;
fval=fold;
while iarm < maxitarm & aflag==1 
        d=-lambda*dd;
        xt=x+d;
        xt=kk_proj(xt,obounds);
        [ft,ifl,ict]=feval(f,xt,h); fct=fct+ict;
        diff_hist=single_point_hist_update(diff_hist,xt,ft,ifl);
        if imfil_least_squares == 1
           frest=ft;
           ft=frest'*frest/2;
        end
        if ifl==1
           ft=fold+abs(fold)*.0001+1.d-12;
        end
        if ft < fval & aflag==1; aflag=0; fval=ft; 
                                 fres=frest; x=xt; 
        end
        if aflag==1; lambda=beta*lambda; end
    iarm=iarm+1;
end
if iarm == maxitarm & aflag == 1 & imfil_verbose == 1
       disp(' line search failure'); [iarm, h]
end
function [hess, xp, fvalp, funs, qfct, iarm, diff_hist] ...
         = imfil_gauss_newton(f, x, fun, jac, xc,  gc, itc, ...
               bounds, h, options)
% IMFIL_GAUSS_NEWTON
% Compute a damped Gauss-Newton step.
%
% function [hess, xp, fvalp, funs, qfct, iarm] ...
%         = imfil_gauss_newton(f, x, fun, jac, xc, gc, itc, ...
%               bounds, h, options)
%
% You may elect to modify this routine if you don't want to spend the
% effort computing the Gauss-Newton model Hessian. The current version
% of imfil.m does not use it.
%
% Inputs: 
%         f = objective function; f returns an residual vector.
%         x = current point.
%       fun = current (vector) residual at x.
%       jac = stencil Jacobian at x.
%        xc = previous point; xc is not used in this function, but
%             that may change if we elect to put a quasi-Newton 
%             method in here.
%        gc = stencil gradient at xc.
%       itc = outer iteration counter.
%    bounds = Nx2 array of bound constraints.
%         h = current scale.
%   options = imfil options structure.
%
% Output:
%      hess = Hessian at x; no update done in here.
%        xp = new point.
%     fvalp = least squares error at xp = funs'*funs/2;
%      funs = vector residual at xp.
%      qfct = cost in function evaluations.
%      iarm = number of step size reductions.
% diff_hist = history data for the Gauss-Newton loop
%        
%
% C. T. Kelley, September 15, 2008
% This code comes with no guarantee or warranty of any kind.
%

%
% Compute stencil gradient (sgrad) and least squares error (fval)
% at current point.
%
funs=fun;
fval=fun'*fun/2;
sgrad=jac'*fun;
%
hess=jac'*jac; % Return the normal equations model Hessian for now
%
imfil_verbose=imfil_optread('verbose',options);
n=length(x);
%
% Get the epsilon-active indices and encode them in a diagonal matrix.
%
epsb=1.d-6;
alist = max((x > bounds(:,1)+epsb) , (x < bounds(:,2)-epsb));
pr=diag(alist);
%
% Give 'em a reduced Hessian to update.
%
hess=eye(n)-pr+pr*hess*pr;
%
% Compute the search direction with a QR factorization.
%
rjac=jac*pr;
[rq,rr]=qr(rjac);
sdir1=(eye(n)-pr)*sgrad;
sdir2=rq'*fun;
sdir2=rr\sdir2;
sdir=sdir1+sdir2;
%
% Bound constrained line search
%
[qfct,xp,fvalp,iarm,fres,diff_hist]=armijo_explore(f, sdir, fval, x, ...
    h, options, bounds);
%
% If the line search fails nothing changes.
%
if length(fres) > 0
   funs=fres;
end
function [hess, xp, fvalp, qfct, iarm, diff_hist] ...
           = imfil_qn_update(f, x, fval, xc, sgrad, gc, hessold, obounds, ...
             itc, h, options)
% IMFIL_QN_UPDATE
% function [hess, xp, fvalp, qfct, iarm, diff_hist] ...
%           = imfil_qn_update(f, x, fval, xc, sgrad, gc, hessold, obounds, ...
%             itc, h, options)
%
% 
% Quasi-Newton update of point and Hessian. This function is never called
% for least squares problems.
%
% Input: 
%        f = objective function.
%        x = current point.
%     fval = f(x). 
%       xc = previous point.
%    sgrad = projected stencil gradient at current point x.
%       gc = projected stencil gradient at previous point xc.
%  hessold = previous model Hessian
%      itc = outer iteration counter.
%  obounds = Nx2 array of scaled bound constraints from imfil_core.
%        h = current scale.
%  options = imfil options structure.
%
% Output:
%     hess = Hessian at x; no update done in here.
%       xp = new point.
%    fvalp = least squares error at xp = funs'*funs/2;
%     qfct = cost in function evaluations.
%     iarm = number of step size reductions.
% diff_hist= history data for the quasi-Newton loop
%
% 
% The steps are (1) update Hessian to obtain hess(current Hessian),
%               (2) use hess, gc, and the bounds to compute the new
%                   search direction,
%               (3) do a line search on that direction,
%               (4) figure out if the search is a success or you should
%                   reduce the scale
%
% C. T. Kelley, September 15, 2008
% This code comes with no guarantee or warranty of any kind.
%
imfil_verbose=imfil_optread('verbose',options);
quasi = imfil_optread('quasi',options);
%
nx=length(xc);
hess=hessold;
%
if itc > 1
%
% Get the epsilon-active indices.
%
epsb=1.d-6;
alist = max((x > obounds(:,1)+epsb) , (x < obounds(:,2)-epsb));
%
switch quasi
   case 1 % BFGS
      hess = bfupdate(x, xc, sgrad, gc, hessold, alist);
   case 2 % SR1
      hess = sr1up(x, xc, sgrad, gc, hessold, alist);
   otherwise % nothing
       hess = eye(nx,nx);
end
end
%
% Search direction
%
sdir = hess\sgrad;
%
% Bound constrained line search
%
[qfct,xp,fvalp,iarm,fres,diff_hist]=armijo_explore(f, sdir, fval, x, ...
    h, options, obounds);
%
%   BFGS update of reduced Hessian, nothing fancy.
%
function hess = bfupdate(x, xc, sgrad, gc, hess,alist)
n=length(x);
pr=diag(alist);
y=sgrad-gc; s=x-xc; z=hess*s;
%
% Turn y into y#.
%
y=pr*y;
if y'*s > 0
   hess = pr*hess*pr + (y*y'/(y'*s)) - pr*(z*z'/(s'*z))*pr;
   hess=eye(n) - pr + hess;
end
if cond(hess) > 1.d6
    hess = eye(n);
end
%
% SR1 update of reduced Hessian.
%
function hess = sr1up(x, xc, sgrad, gc, hess,alist)
n=length(x);
pr=diag(alist);
y=sgrad-gc; s=x-xc; 
z=y - hess*s;
%
% Turn y into y#.
%
y=pr*y; z=pr*z;
if z'*s ~=0
        ptst = z'*(hess*z)+(z'*z)*(z'*z)/(z'*s);
        if ptst > 0 
            hess = pr*hess*pr + (z*z')/(z'*s); 
            hess = eye(n) - pr + hess;
        end
end
if cond(hess) > 1.d6
    hess = eye(n);
end
function [grad,best_value,best_point,icount,sflag,tvar,diff_hist,jac]=...
                 stencil_diff(x,f,dx,fc,bounds,options,h,complete_history)
% STENCIL_DIFF 
% Stencil derivative for optimization and nonlinear least squares. 
% This function also tests for best point in stencil.
%
% function [grad,best_value,best_point,icount,sflag,tvar,diff_hist,jac]=...
%                 stencil_diff(x,f,dx,fc,bounds,options,h)
%
%   Input:  x  = current point
%           f  = objective function (or vector residual for
%                nonlinear least squares)
%           dx = scaled difference directions
%           fc = current function value
%       bounds = bounds (duh!)
%      options = imfil option structure
%            h = current scale
%
%   Output: 
%         grad = stencil gradient; NaN if every point on the stencil is bad
%  best_value  = best value in stencil
%  best_point  = best point in stencil
%       icount = counter of calls to expensive part of function
%                 the function needs to return this.
%        sflag = 0 current point is best in the stencil
%                     This means stencil failure.
%
%              = 1 means there's a better point on the stecil.
%
%         tvar = variation on stencil = max - min
%      
%    diff_hist = every function evaluation for points satisfying the 
%                 bounds stored in a nice struct. I will write some
%                 tools for messing with this stuff, so you don't have
%                 to deal with it directly.
%
%          jac = stencil Jacobian for least squares problems.
%
% C. T. Kelley, July 30, 2009
% This code comes with no guarantee or warranty of any kind.
%
parallel=imfil_optread('parallel',options);
least_squares=imfil_optread('least_squares',options);
%
% Poll the stencil and collect some data.
%
[best_value,best_point,icount,sgood,good_points,good_values,...
         good_dx,good_df,failed_points] = ...
           imfil_poll_stencil(x,f,dx,fc,bounds,options,h,complete_history);
%
% Use the points and values to get the stencil statisitcs.
%
[sflag,best_value,best_point,tvar,diff_hist]=...
              collect_stencil_data(good_points,good_values,failed_points,...
                                   best_point,best_value,fc,options);
%
% sgood = 0 means there are no good points. Shrink time!
%
if sgood > 0
%
% Estimate the derivative.
% The idea is that fprime*dx approx df, where
% fprime   = m x n is the gradient-transpose or the Jacobian
% good_dx  = n x pnew (pnew <= vsize) is the matrix of good steps
% good_df  = m x pnew is the collection of good function values a x + dx
% 
% So the least squares problem to be solved is 
%         min || fprime * good_dx - good_df ||
% for the columns of fprime. The minimum norm solution is
% fprime = good_df*pinv(good_dx).
%
   pdx=pinv(good_dx);
   fprime=good_df*pdx;
   if least_squares == 0
      grad=fprime';
      jac=[];
   else
      jac=fprime;
      grad=fprime'*fc;
   end
else
   grad=NaN*x;
   jac=[];
end
function px = kk_proj(x,bounds)
% KK_PROJ
% Projection onto the feasible hyperrectangle.
%
% function px = kk_proj(x,bounds)
%
% Not exciting stuff.
%
% C. T. Kelley, September 15, 2008
% This code comes with no guarantee or warranty of any kind.
%
ndim=length(x);
px=zeros(ndim,1);
px=min(bounds(:,2),x);
px=max(bounds(:,1),px);
function fvals=f_to_vals(funs,least_squares)
%
% F_TO_VALS
% Evaluate fvals = f^T f/2 when f is a vector least squares residual.
%
% function fvals=f_to_vals(funs,least_squares)
%
% There is no reason you'd want to mess with this.
%
% C. T. Kelley, September 15, 2008
% This code comes with no guarantee or warranty of any kind.
%
[m,n]=size(funs);
if m == 0 && n == 0
   fvals=[];
end
fvals=zeros(n,1);
if least_squares == 1
  for i=1:n
    fvals(i)=funs(:,i)'*funs(:,i)/2;
  end
else
    fvals=funs;
end

function vout=imfil_augment_directions(x,vin,h,options,bounds)
% IMFIL_AUGMENT_DIRECTIONS
% Enrich the stencil with random directions and a the user-supplied
% new_direcions function.
%
% C. T. Kelley, August 13, 2009.
% This code comes with no guarantee or warranty of any kind.
%
% Add the random directions.
%
random_stencil = imfil_optread('random_stencil',options);
vout = random_augment(vin,random_stencil);
%
% See if there's a new_directions function.
%
new_directions = imfil_optread('add_new_directions',options);
lnew=length(new_directions);
if lnew > 0
    dbv=bounds(:,2)-bounds(:,1); db=diag(dbv);
    unscaled_x = db*x + bounds(:,1);
    unscaled_v = db*vin;
    unscaled_vnew=feval(new_directions, unscaled_x, h, unscaled_v);
    [mv,nv]=size(unscaled_vnew);
    if mv > 0
        vnew=inv(db)*unscaled_vnew;
        [mv,nv]=size(vnew);
        for i=1:nv
           vnew(:,i)=vnew(:,i)/norm(vnew(:,i));
        end
        vout=[vout, vnew];
    end
end

function vout=random_augment(vin,k)
% RANDOM_AUGMENT
% Add k random unit vectors to the stencil.
%
% This makes some of the theory work, but it is not always
% a good idea to do this. Adding more vectors makes stencil failure
% less likely, which is not good if you have a full stencil and you're
% wasting lots of time in the line search.
%
% Adding random vectors makes more sense if the minimizer is 
% on a hidden or explicit constraint boundary.
%
% C. T. Kelley, August 2, 2009.
% This code comes with no guarantee or warranty of any kind.
%
[n,nn]=size(vin);
if k==0
   vout=vin;
end
rv=rand(n,k);
for ir=1:k
   rv(:,ir)=rv(:,ir)/norm(rv(:,ir));
end
vout=[vin,rv];

function [best_value,best_point,icount,sgood,good_points,good_values,...
         good_dx,good_df,failed_points] = ...
              imfil_poll_stencil(x,f,dx,fc,bounds,options,h,complete_history)
% IMFIL_POLL_STENCIL
% Poll the stencil and organize the results.
% I do not think you will want to modify this code.
%
% function [best_value,best_point,icount,sgood,good_points,good_values,...
%           good_dx,good_df,failed_points] = ...
%                          imfil_poll_stencil(x,f,dx,fc,bounds,options,h)
%
% Input:
%        x = center of stencil
%        f = objective function
%       dx = h*V = array of scaled directions
%       fc = f(x)
%   bounds = N x 2 array of lower/upper bounds
%  options = imfil options structure
%        h = current scale
%
% Output:
%     best_value    = lowest value of f on the stencil
%     best_point    = best_value = f(best_point)
%     icount        = cost in units of calls to f
%     sgood         = number of successful evaluations
%     good_points   = list of points at which f returned a value
%     good_values   = f(good_points)
%     good dx       = vector of good_points - x
%     good df       = vector of f(good_points) - fc
%     failed_points = list of points at which f did not return a value
%
% C. T. Kelley, Aug 13, 2009
% This code comes with no guarantee or warranty of any kind.
%
parallel=imfil_optread('parallel',options);
least_squares=imfil_optread('least_squares',options);
[n,vsize]=size(dx);
%
% Record the best point and best function value in the stencil.
% Get started by using the center point.
%
best_point=x;
if least_squares == 1
best_value=fc'*fc/2; % best point in the stencil
    else
best_value=fc;
end
worst_value=best_value;
%
iflag=zeros(vsize,1);
m=length(fc); % m > 1 tells me it's a least squares problem
fp=zeros(m,vsize);
failed_points=[];
good_points=[];
good_values=[];
sgood=0;
fval=zeros(vsize,1);
icount=0;
%
% fp(:,i) and fval(i) are not defined outside of the bounds
% or if iflag(i)=1.
%
% First cull the points which violate the bound constraints.
%
% Collect the feasible points, differences, and functions in xp1 and dx1.
%
%
pold=0;
dx1=[];
xp1=[];
for i=1:vsize
  xp(:,i)=x+dx(:,i);
  if isok(xp(:,i),bounds)
     pold=pold+1;
     dx1=[dx1,dx(:,i)];
     xp1=[xp1,xp(:,i)];
  end
end
fp1=fp(:,1:pold);
xp=xp1; dx=dx1; fp=fp1;
%
% Query the complete_history structure to see if you've evaluated f
% at any of these points before.
%
[oldindex,oldpoints,oldvalues,oldflags]= ...
               scan_history(complete_history,xp,fp,dx);
newindex=~oldindex;
xp=xp1(:,newindex);
iflago=zeros(1,vsize);
if sum(oldindex) > 0
fp(:,oldindex)=oldvalues;
iflago(oldindex)=oldflags;
end
%xp=xpout; dx=dxout; fp=fpout;
%
% Evaluate f, in parallel if possible. Flag the failed points.
%
pnew=sum(newindex);
fp1=[];
iflag=[];
if parallel == 0
    for i=1:pnew
        [fpx,iflagx,ict]=feval(f,xp(:,i),h);
        fp1=[fp1,fpx];
        iflag=[iflag,iflagx];
        icount=icount+ict;
    end
else
    if pnew > 0
        [fp1,iflag,ictrp]=feval(f,xp,h);
        icount=icount+sum(ictrp);
    end
end
if pnew > 0
   fp(:,newindex)=fp1;
   iflago(newindex)=iflag;
end
fp1=fp;
iflag=iflago;
%
% Identify the failed points.
%
ibad=(iflag(1:pold)==1);
if sum(ibad) > 0
   failed_points=xp1(:,ibad);
end
%
% Store the good points and their function values.
%
igood=(iflag(1:pold)==0);
sgood=sum(igood);
good_dx=[];
good_df=[];
if sgood > 0
   good_points=xp1(:,igood);
   if least_squares == 1
      good_fp=fp1(:,igood);
   else
      good_fp=fp1(igood);
   end
   good_dx=dx1(:,igood);
   for ig=1:sgood
      if least_squares == 1
         good_df(:,ig)=good_fp(:,ig)-fc;
      else
         good_df(ig)=good_fp(ig)-fc;
      end
   end
%   good_values=f_to_vals(good_fp,least_squares);
   good_values=good_fp;
end

%
% test x for feasibility re bounds
%
function lisok=isok(x,bounds)
il=min(x >= bounds(:,1));
iu=min(x <= bounds(:,2));
isok=min(il,iu);
lisok = (isok == 1);

%
% Update the complete_history structure after a many calls to f.
%
function new_hist=many_point_hist_update(old_hist,diff_hist,bounds)
new_hist=old_hist;
if nargin == 3
   diff_hist.good_points=mscale(diff_hist.good_points,bounds);
   diff_hist.failed_points=mscale(diff_hist.failed_points,bounds);
end
new_hist.good_points=[new_hist.good_points,diff_hist.good_points];
new_hist.good_values=[new_hist.good_values,diff_hist.good_values];
new_hist.failed_points=[new_hist.failed_points,diff_hist.failed_points];

function ds=mscale(d,bounds)
%
% Unscale an array of points from imfil_core to the original bounds.
%
[md,nd]=size(d);
dbounds=bounds(:,2)-bounds(:,1);
db=diag(dbounds);
qbounds=diag(bounds(:,1));
ds=d;
if md > 0
   dlow=qbounds*ones(size(d));
   ds=db*d+dlow;
end

%
% Update a history structure after a single call to f.
%
function new_hist=single_point_hist_update(old_hist,x,fout,ifail,bounds)
new_hist=old_hist;
if nargin == 5
   x=bscale(x,bounds);
end
if ifail == 1
    new_hist.failed_points=[old_hist.failed_points,x];
else
    new_hist.good_points=[old_hist.good_points,x];
    new_hist.good_values=[old_hist.good_values,fout];
end

function sx=bscale(x,bounds)
%
% Unscale a vector from imfil_core to the original bounds.
%
dbounds=bounds(:,2)-bounds(:,1);
db=diag(dbounds);
sx=db*x+bounds(:,1);
function [oldindex,oldpoints,oldvalues,oldflags] = ...
               scan_history(complete_history,xp,fp,dx)
[mp,np]=size(xp);
imold=zeros(1,np);
oldpoints=[];
oldvalues=[];
oldflags=[];
for i=1:np
    [fpt,ift]=scal_complete_history(complete_history,xp(:,i));
    if ift > -1
       oldpoints=[oldpoints,xp(:,i)];
       oldvalues=[oldvalues,fpt];
       oldflags=[oldflags,ift];
       imold(i)=1;
    end
end
oldindex=(imold==1);


function [fhist,iflaghist]=scal_complete_history(complete_history,x)
fhist=[];
iflaghist=-1;
losers=complete_history.failed_points;
[ml,nl]=size(losers);
winners=complete_history.good_points;
[mw,nw]=size(winners);
iquit=0;
if nw > 0
for i=1:nw
  d=norm(x-winners(:,i),inf);
  if d < 1.d-12
     iquit=1;
     fhist=complete_history.good_values(:,i);
     iflaghist=0;
     break;
  end
end
end
if iquit == 0 & nl > 0
for i=1:nl
   d=norm(x-losers(:,i),inf);
   if d < 1.d-12
     fhist=NaN;
     iquit=1;
     iflaghist=1;
     break;
  end
end
end
function valout=imfil_optread(str,optin)
% IMFIL_OPTREAD
%
% Reads the options one at a time. This is an internal
% function for imfil.m. 
%
% Input: 
%      str = name of option
%    optin = options structure
%
% Output: 
%   valout = value of option
%
% Example: fscale=imfil_optread('fscale',optin);
%
% The options array is fully documented in the comment lines to
% imfil_optset. Type `help imfil_optset' for the details.
%
% There is no error checking in imfil_optread. It's coming.
% 
%
% C. T. Kelley, September 15, 2008
% This code comes with no guarantee or warranty of any kind.
%

str=lower(str);
if nargin <= 2
   n=-1;
end
parms=fieldnames(optin);
nopt=length(parms);
did_ok=0;
for p=1:nopt
    if strcmp(str,parms(p));
       valout=getfield(optin,str);
       did_ok=1;
    end
end

if did_ok == 0;
   error('Error in imfil_optread');
end