function [histoutl, histoutg, histoutlx, histoutgx]=driver_lm
% DRIVER_LM
%
% Test for the executive_function option. Solve the PID problem with
% a Levenberg-Marquardt iteration.
%
% C. T. Kelley, Jan 24, 2011.
%
% This code comes with no guarantee or warranty of any kind.
%
% Nonlinear least squares formuation.
%
% This is an overdetermined least squares problem with n=2 and m >> n
%
% pid_parms contains the zero-residual solution to the noise-free problem
% pid_tol is the tolerance given to ode15s
%
m=100; t0=0; tf=10;
%
% Construct the data for the integration.
% pid_data is a sampling of the "true" solution.
%
pid_parms=[1,1]'; pid_y0=[10,0]'; pid_tol=1.d-3;
time_pts=(0:m)'*(tf-t0)/m+t0; 
%
% Find the analytic solution.
%
pid_data=exact_solution(time_pts,pid_y0,pid_parms);
%
% Pack the data into a structure to pass to serial_pidlsq
%
pid_info=struct('pid_y0',pid_y0,'pid_tol',pid_tol,...
                'time_pts',time_pts,'pid_data',pid_data);
%
% Set the bounds, budget, initial iterate for imfil.m.
%
bounds=[2 20; 0 5];
x0=[5,5]'; budget= 50;
%
% Call imfil. Global minimum out of bounds.
%
options=imfil_optset('least_squares',1,...
   'executive_function',@lev_mar_exec,'executive_data',1.0);
[x,histoutl]=imfil(x0,@serial_pidlsq,budget,bounds,options,pid_info);
options=imfil_optset('least_squares',1);
[x,histoutg]=imfil(x0,@serial_pidlsq,budget,bounds,options,pid_info);
%
% Call imfil. Global minimum within bounds.
%
bounds=[0 20; 0 5];
budget=200;
options=imfil_optset('least_squares',1,'scaledepth',20,...
   'executive_function',@lev_mar_exec,'executive_data',1.0);
[x,histoutlx]=imfil(x0,@serial_pidlsq,budget,bounds,options,pid_info);
options=imfil_optset('least_squares',1,'scaledepth',20);
[x,histoutgx]=imfil(x0,@serial_pidlsq,budget,bounds,options,pid_info);
%
% Left side plot.
%
figure(1)
p1=subplot(1,2,1);
p2=plot(histoutg(:,1),histoutg(:,2),'-',histoutl(:,1),histoutl(:,2),'--');
axis('square')
set(p1,'FontSize',14,'XLim',[0 50],'YLim',[20 70]);
set(p2,'LineWidth',1.5,'Color','black');
legend('G-N','L-M');
xlabel('Calls to integrator'); 
ylabel('Value of f');
title('Constraints Active');
%
% Right side plot.
%
p1=subplot(1,2,2);
p2=semilogy(histoutgx(:,1),histoutgx(:,2),'-',...
           histoutlx(:,1),histoutlx(:,2),'--');
set(p1,'FontSize',14,'XTick',[0 50 100 150 200],...
            'XLim',[0 200],'YLim',[1.d-3 100],'YTick',[1.d-3 1.d-1 1.d2]);
axis('square')
set(p2,'LineWidth',1.5,'Color','black');
title('Constraints Inactive');
%legend('Gauss-Newton','Levenberg-Marquardt');
xlabel('Calls to integrator');
%ylabel('Value of f');