function driver_parallel_pid
% DRIVER_PARALLEL_PID
%
% Example script for imfil.m; nonlinear least squares and parameter id.
% This generates the two figures in the section "Looking at the Parallel
% Results".
%
% C. T. Kelley, May 23, 2010
%
% This code comes with no guarantee or warranty of any kind.
%
% Solving the PID problem with implicit filering as simply as possible.
% Compare optimization/nonlinear least squares formulations with parallel
% and serial algorithms. We consider both unconstrained and constrained
% problems and a couple different choices for imfil's options.
%
% We use lots of global variables here, but plan to do this right before
% the final version. 
%
% This is an overdetermined least squares problem with n=2 and m = 100 >> n
%
% pid_parms contains the zero-residual solution to the noise-free problem
% pid_tol is the tolerance given to ode15s
%
hist_optim=cell(4,2);
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);
%
figure(1);
%
% Constrained problem with too few scales.
%
case_number=1;
bounds=[2 20; 0 5];
x0=[5,5]'; budget= 100;
options=imfil_optset('scaledepth',5);
[histso,histsl,histpo,histpl]=run_case(x0,budget,bounds,options,pid_info);
plot_cases(histso,histsl,histpo,histpl,1);
hist_case_1=histso;
legend('serial-opt','serial-lsq','parallel-opt','parallel-lsq');
ylabel('Value of f');
title('Default scales; constraints active');
%
% Constrained problem with more scales.
%
case_number=2;
bounds=[2 20; 0 5];
x0=[5,5]'; budget= 100;
options=imfil_optset;
[histso,histsl,histpo,histpl]=run_case(x0,budget,bounds,options,pid_info);
plot_cases(histso,histsl,histpo,histpl,2);
hist_case_2=histso;
title('More scales; constraints active');
%
% Unconstrained problem with thin budget.
%
case_number=3;
bounds=[0 20; 0 5];
x0=[5,5]';
options=imfil_optset;
[histso,histsl,histpo,histpl]=run_case(x0,budget,bounds,options,pid_info);
plot_cases(histso,histsl,histpo,histpl,case_number,'semi-log');
hist_case_3=histso;
xlabel('Calls to integrator'); ylabel('Value of f');
title('Low budget; constraints inactive');
%
% Unconstrained problem with sufficient budget.
%
case_number=4;
bounds=[0 20; 0 5];
x0=[5,5]';
budget=200;
options=imfil_optset('scaledepth',20);
[histso,histsl,histpo,histpl]=run_case(x0,budget,bounds,options,pid_info);
plot_cases(histso,histsl,histpo,histpl,case_number,'semi-log');
hist_case_4=histso;
xlabel('Calls to integrator');
title('High budget; constraints inactive');
%

function [histso,histsl,histpo,histpl]=run_case(x0,budget,bounds,...
                 options,pid_info)
for pflag=0:1
for lflag=0:1
opt_case=lflag+2*pflag;
%
% Call imfil,
%
options=imfil_optset('least_squares',lflag,'parallel',pflag,options);
if lflag==0
[x,histout]=imfil(x0,@pidobj,budget,bounds,options,pid_info);
else
[x,histout]=imfil(x0,@parallel_pidlsq,budget,bounds,options,pid_info);
end
%
switch opt_case
     case 0 % serial + optimization formulation
         histso=histout;
     case 1 % serial + least squares
         histsl=histout;
     case 2 % parallel + optimization formulation
         histpo=histout;
     case 3 % parallel + least squares
         histpl=histout;
end
%
end
end


function plot_cases(histso,histsl,histpo,histpl,np,plot_form);
p1=subplot(2,2,np);
if nargin == 5
p2=plot(histso(:,1),histso(:,2),'-',histsl(:,1),histsl(:,2),'--',...
        histpo(:,1),histpo(:,2),'-o',histpl(:,1),histpl(:,2),'-x');
set(p1,'FontSize',14,'XTick',[0 50 100 150 200],'XLim',[0 100],'YLim',[20 70]);
else
p2=semilogy(histso(:,1),histso(:,2),'-',histsl(:,1),histsl(:,2),'--',...
        histpo(:,1),histpo(:,2),'-o',histpl(:,1),histpl(:,2),'-x');
set(p1,'FontSize',14,'XTick',[0 50 100 150 200],'XLim',[0 200],...
    'YLim',[1.d-4, 1.d2],'YTick',[1.d-4 1.d-2 1.d2]);
end
set(p2,'LineWidth',1.5,'Color','black');