%--------------------------------------
% example of 1D cost function minisation
%--------------------------------------
% 6/12/13	N. Teanby
%--------------------------------------
close all
clear all

%--------------------------------------
% make a fake data set
% (NB usually this would be read in from a datafile)
%--------------------------------------
lambda = 60;
a0     = 10;

t = 0:10:300;
n = length(t);

% calculate a decay curve and add some random noise
a     = zeros(n,1);
noise = 0.5*rand(n,1);
for j=1:n
  a(j) = a0*exp(-t(j)/lambda) + noise(j);
end

% define the errorbars
e = ones(1,n)*0.5;

% first point has zero error
a(1) = a0;
e(1) = 0;

%--------------------------------------


%--------------------------------------
% (1) plot fake data
%--------------------------------------
figure(1)
errorbar(t,a,e,'ok')
xlabel('Time (s)')
ylabel('Activity')
xlim([0 305])
grid on


%--------------------------------------
% grid search for minimum phi
%--------------------------------------
% (2) Function: a = a0*exp(-t/lambda) 
%               based on the fact that it is radiactive decay

% (3) Range:    50 - 100 seconds seems sensible based on the plot

% (4) Grid:     make an evenly spaced grid of parameters
lambda_trial = 50:2:100;
n_trial      = length(lambda_trial);

% (5) Search:   loop around the paramters
% make a blank vector to store phi in
phi          = zeros(n_trial,1);
for i=1:n_trial
  
  % (5a) Evaluate the function, ie the trial model fit
  a_trial = zeros(n,1);
  for j=1:n
    a_trial(j) = a0*exp(-t(j)/lambda_trial(i));
  end
  
  % (5b) Calculate the cost function
  % calc the misfit function for this value of lambda_trial
  % [dont include the first point as it has no error bar
  % so would give a divide by zero error. Usually you
  % would include the first data point though.]
  phi(i)  = 0;
  for j=2:n
    phi(i) = phi(i) + ( (a(j)-a_trial(j)) / e(j) )^2;
  end
  
end

% (6) Minimise
% find minimum phi and its index imin
[phi_min,imin]=min(phi);

%--------------------------------------
% (7) plot phi variation
%--------------------------------------
figure(2)
plot(lambda_trial,phi,'-ok')
xlabel('\lambda (s)')
ylabel('\phi(\lambda)')
grid on

%--------------------------------------
% (8) plot best fit and a few of the trial fits
%--------------------------------------
figure(3)
% data
errorbar(t,a,e,'ok')
xlabel('Time (s)')
ylabel('Activity')
xlim([0 305])

% overlay fits
hold on
% best fit
a_trial = a0*exp(-t/lambda_trial(imin));
plot(t,a_trial,'-b')

% fits for end points
a_trial = a0*exp(-t/lambda_trial(1));
plot(t,a_trial,'-r')
a_trial = a0*exp(-t/lambda_trial(n_trial));
plot(t,a_trial,'-g')

hold off
grid on

% legend
label_best  = strcat('Best fit, \lambda=',num2str(lambda_trial(imin)));
label_start = strcat('Start, \lambda=',num2str(lambda_trial(1)));
label_end   = strcat('End, \lambda=',num2str(lambda_trial(n_trial)));
legend('Data',label_best,label_start,label_end)






