function [x,histout,chistout]= history_test
% HISTORY_TEST
%
% Example script for imfil.m; nonlinear least squares and parameter id.
% Plot some things from the complete_history structure.
% This generates the final figure in Chapter 1: "Getting Started"
%
% C. T. Kelley, July 21, 2009
%
% 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);
bounds=[2 20; 0 5];
bounds=[0 5; 0 5];
x0=[5,5]'; budget= 400;
for il=0:1
least_squares=1;
options=imfil_optset('least_squares',least_squares,'parallel',il);
if least_squares==0
    [x,histout,chistout]=imfil(x0,'pidobj',budget,bounds,options);
else
    [x,histout,chistout]=imfil(x0,'pidlsq',budget,bounds,options);
end
figure(1)
subplot(1,2,il+1)
if least_squares == 0
   vals=chistout.good_values;
else
   [mv,nv]=size(chistout.good_values);
   vals=zeros(1,nv);
   for i=1:nv
      vals(i)=chistout.good_values(:,i)'*chistout.good_values(:,i)/2;
   end
end
plot(vals)
figure(2)
p1=subplot(1,2,il+1);
p2=plot(chistout.good_points(1,:),chistout.good_points(2,:),'*');
%[mg,ng]=size(chistout.good_points)
%p2=plot(chistout.good_points(1,1),chistout.good_points(2,1),'*');
%hold on;
%for i=2:ng
%plot(chistout.good_points(1,i),chistout.good_points(2,i),'*');
%end
axis('square');
if il==0
    title('Serial');
else
    title('Parallel');
end
set(p1,'FontSize',14,'XLim',bounds(1,:),'YLim',bounds(2,:));
end

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);