function pid_example_chapter_1
% PID_EXAMPLE_CHAPTER_1
%
% Example script for imfil.m; nonlinear least squares and parameter id.
% This generates the two figures in Chapter 1: "Getting Started"
%
% C. T. Kelley, July 28, 2008
%
% 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 to avoid pain
%
% 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
%
global pid_data time_pts pid_parms pid_y0 pid_tol;
isr1=0;
pflag=0;
hist_optim=cell(4,2);
m=100; t0=0; tf=10;
%
% Set the globals for the objective.
% pid_data is a sampline of the "true" solution.
%
pid_parms=[1,1]'; pid_y0=[10,0]'; pid_tol=1.d-6;
time_pts=(0:m)'*(tf-t0)/m+t0;
%
pid_data=yfunex(time_pts,pid_y0);
%
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);
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);
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);
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);
plot_cases(histso,histsl,histpo,histpl,case_number,'semi-log');
hist_case_4=histso;
xlabel('Calls to integrator');
title('High budget; constraints inactive');
%
figure(2);
p1=subplot(2,2,1);
p2=plot(hist_case_1(:,1),hist_case_1(:,2),'-');
set(p2,'LineWidth',1.5,'Color','black');
set(p1,'FontSize',14,'XTick',[0 50 100 150 200],'XLim',[0 150],'YLim',[20 70]);
ylabel('Value of f');
title('Default scales; constraints active');
p1=subplot(2,2,2);
p2=plot(hist_case_2(:,1),hist_case_2(:,2),'-');
set(p2,'LineWidth',1.5,'Color','black');
set(p1,'FontSize',14,'XTick',[0 50 100 150 200],'XLim',[0 150],'YLim',[20 70]);
title('More scales; constraints active');
p1=subplot(2,2,3);
p2=semilogy(hist_case_3(:,1),hist_case_3(:,2),'-');
set(p1,'FontSize',14,'XTick',[0 50 100 150 200],'XLim',[0 200],...
    'YLim',[1.d-10, 1.d2],'YTick',[1.d-10 1.d-5 1.d-2 1.d2]);
set(p2,'LineWidth',1.5,'Color','black');
xlabel('Calls to integrator'); ylabel('Value of f');
title('Low budget; constraints inactive');
p1=subplot(2,2,4);
p2=semilogy(hist_case_4(:,1),hist_case_4(:,2),'-');
set(p1,'FontSize',14,'XTick',[0 50 100 150 200],'XLim',[0 200],...
    'YLim',[1.d-10, 1.d2],'YTick',[1.d-10 1.d-5 1.d-2 1.d2]);
set(p2,'LineWidth',1.5,'Color','black');
xlabel('Calls to integrator');
title('High budget; constraints inactive');


function [histso,histsl,histpo,histpl]=run_case(x0,budget,bounds,options)
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);
else
[x,histout]=imfil(x0,'pidlsq',budget,bounds,options);
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 150],'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-10, 1.d2],'YTick',[1.d-10 1.d-5 1.d-2 1.d2]);
end
set(p2,'LineWidth',1.5,'Color','black');

       
function yfex=yfunex(range,y0);
% YFUNEX
% Compute the "true" solution at the data points.
%
global pid_parms;
c=pid_parms(1);
k=pid_parms(2);
xip=.5*(-c+sqrt(c*c-4*k));
xim=.5*(-c-sqrt(c*c-4*k));
b=[1,1;xip,xim];
ac=b\y0;
yfc(:,1)=ac(1)*exp(xip*range) + ac(2)*exp(xim*range);
yfc(:,2)=ac(1)*xip*exp(xip*range) + ac(2)*xim*exp(xim*range);
yfex=real(yfc);