global mu_l mu_a c_ea c_el c_pa b mydatafit myDays myb

% Beetles, Cannibalism, and Chaos
% Solution to Problem 3
% Dianne O'Leary 12/2006
%

disp('Read beetledata')
beetledata;

% Part (a). 
% Perform the least squares fit to determine the size
% parameters for each of the four experimental colonies.
% 
% We use lsqnonlin, with lower (lb) and upper (ub) bounds on the
% parameters.  The parameters are ordered

%        c_el c_ea c_pa   b mu_L mu_A
param0 = [0.1; 0.1; 0.1;  6; 0.5; 0.1];
lb     = [0  ;   0;   0; .1;   0;   0];
ub     = [1  ;   1;   1;  9;   1;   1];

% Loop over the colonies.

for mycase=1:4,

% Set mydatafit to the number of larvae, pupae, and adults for this colony.

  mydatafit = [Data(:,1+mycase),Data(:,5+mycase),Data(:,9+mycase)];
  myDays = Data(:,1) + 1;

% Perform the least squares fit.

  [paramlog(:,mycase),rnormlog(mycase)] = lsqnonlin(@evalls2,param0,lb,ub);

% Part (b).
% Compare the predictions of this model to that of Dennis et al.
% using a plot.  A printed summary will be made later.

  [diffs,L,P,A]=evalls2(paramlog(:,mycase));
  [diffs,Ld,Pd,Ad]=evalls2(param_dl(:,mycase));
  figure(mycase)
  subplot(3,1,1)
  plot(1:20,mydatafit(:,1),'r-', 1:20,A, 'rs', 1:20,Ad,'r+')
  legend('Adults')
  title('Predictions by Dennis et al. (+) and our calculations (square)')
  subplot(3,1,2)
  plot(1:20,mydatafit(:,2),'b:', 1:20,P ,'bs', 1:20,Pd,'b+')
  legend('Pupae')
  subplot(3,1,3)
  plot(1:20,mydatafit(:,3),'g--', 1:20,L ,'gs', 1:20,Ld,'g+')
  legend('Larvae')

end

% Part (c).
% For the second colony, explore the sensitivity of the model for
% the second colony.

mycase = 2;
mydatafit = [Data(:,1+mycase),Data(:,5+mycase),Data(:,9+mycase)];

% First: A backward error analysis: how sensitive is the model to
% changes in the b parameter?
% Compute the residual after changing just the b parameter,
% leaving the other parameters fixed at their optimal values.

myparam = paramlog(:,mycase);
myparam(4) = 0;
e4 = [0;0;0;1;0;0];

figure(5)
pert = 4:.25:12;
for j=1:length(pert)
   r(j) = norm(evalls2(myparam+pert(j)*e4));
end
plot(pert,r)
xlabel('b')
ylabel('Residual norm')
axis([min(pert),max(pert) 7 8])
title('How sensitive is the residual to changes in b?')

% Second, a forward error analysis: how sensitive is b to
% small changes in the measured data?
% Repeat the fit, after adding 250 samples of normally
% distributed error to the log of the counts, mean 0, standard deviation 1.

for k=1:250,
   mydatafit = [Data(:,1+mycase),Data(:,5+mycase),Data(:,9+mycase)]; 
   mydatafit = exp(log(mydatafit) + randn(size(mydatafit)));
   [paramlogpert,rnorm(k)] = lsqnonlin(@evalls2,paramlog(:,2),lb,ub);
   b(k) = paramlogpert(4);
end

% Plot the results; we will print a summary sentence later.

figure(6)
plot(b,zeros(1,250),'*') 
xlabel('b estimates')
title('Results of random perturbations on data for the second colony')

% Finally, a homotopy to see how much the residual changes as
% we change b.
% We fix b and do a least squares fit to find optimal values for the
% other 5 parameters.  We plot the resulting residual norm as a function
% of b.  As a side effect, we actually find a better solution than
% we had before.

myb_all = 1:.5:50;
lbs = lb([1:3,5:6]);
ubs = ub([1:3,5:6]);
param = paramlog([1:3,5:6,2]);
mymin = 20;

for k=1:length(myb_all)
   myb = myb_all(k);
   mydatafit = [Data(:,1+mycase),Data(:,5+mycase),Data(:,9+mycase)];
   [param,rnorm(k)] = lsqnonlin(@evallsh,param,lbs,ubs);
   res(k) = norm(evallsh(param));
   if (res(k) < mymin)
      mymin = res(k);
      myminparam= param;
      myminb = myb;
    end
end

figure(7)
plot(myb_all,res,'bo',myb_all,ones(size(myb_all))*mymin*1.05,'r-') 
legend('Best','1.05*Optimal')
xlabel('b')
ylabel('Residual norm')
title('Comparison of optimal residual norm with best found for fixed b')


% Print some summary results 

disp(' ')
disp('Parameter estimates from least squares:')
disp(' ')
disp('colony        c_el     c_ea     c_pa        b     mu_L     mu_A  residual')
for mycase=1:4,
  mydatafit = [Data(:,1+mycase),Data(:,5+mycase),Data(:,9+mycase)];
  disp(sprintf('mine:   %d %f %f %f %f %f %f %f', ...
       mycase, paramlog(:,mycase)',norm(evalls2(paramlog(:,mycase)))))
  disp(sprintf('theirs: %d %f %f %f %f %f %f %f', ...
       mycase, param_dl(:,mycase)',norm(evalls2(param_dl(:,mycase)))))
end

disp(' ')
disp(' ')
disp('Residual norms:')
for mycase=1:4,
  mydatafit = [Data(:,1+mycase),Data(:,5+mycase),Data(:,9+mycase)];
  disp(sprintf('Dataset %d, norm of data vector %f, my residual %f, Dennis et al. residual %f', ...
       mycase, norm(log(mydatafit(2:end,:)),'fro'), ...
       norm(evalls2(paramlog(:,mycase))), ...
       norm(evalls2(param_dl(:,mycase)))   ))
end

disp(' ')
disp('Parameter estimates from homotopy for Colony 2:')
disp(' ')
mycase = 2;
best = [myminparam(1:3),myminb,myminparam(4:5)];
mydatafit = [Data(:,1+mycase),Data(:,5+mycase),Data(:,9+mycase)];
disp(sprintf('mine:   %d %f %f %f %f %f %f %f', ...
     mycase, best,norm(evalls2(paramlog(:,mycase)))))
disp(sprintf('theirs: %d %f %f %f %f %f %f %f', ...
     mycase, param_dl(:,mycase)',norm(evalls2(param_dl(:,mycase)))))

[diffs,L,P,A]=evalls2(best);
  [diffs,Ld,Pd,Ad]=evalls2(param_dl(:,mycase));
  figure(8)
  subplot(3,1,1)
  plot(1:20,mydatafit(:,1),'r-', 1:20,A, 'rs', 1:20,Ad,'r+')
  legend('Adults')
  title('Predictions by Dennis et al. (+) and our calculations (square)')
  subplot(3,1,2)
  plot(1:20,mydatafit(:,2),'b:', 1:20,P ,'bs', 1:20,Pd,'b+')
  legend('Pupae')
  subplot(3,1,3)
  plot(1:20,mydatafit(:,3),'g--', 1:20,L ,'gs', 1:20,Ld,'g+')
  legend('Larvae')



% Print a summary of the results of random perturbations in the data.

disp(' ')
disp(sprintf('When the data is randomly perturbed, the estimate of b'))
disp(sprintf('for the second colony ranges from %f to %f.', min(b),max(b)))