% Laboratory in Oceanography - Data and Methods
% SMAST, UMass Dartmouth
%
% For testing fitting of a Gaussian using equally weighted points
%
% Written by Miles A. Sundermeyer, 7/6/06
% Last modified by Miles A. Sundermeyer, 5/8/09

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create a synthetic data set with noise
x = 1:10;
x = 1:20;
noise = rand(size(x)) - 0.5;

A = 1.1;
xo = 9.9;
sigma = 2.2;
actfitparams = [A xo sigma];

y = noise + A*exp(-(x-xo).^2 / (2*sigma^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
clf
plot(x,y,'bo','linewidth',2)
grid on
xlabel('x')
ylabel('y')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit a Gaussian using either Matlab's lsqnonlin fit or 'by hand'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(1)		% use Matlab's lsqnonlin
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % set an initial guess
  if(1)
    [Aguess, I] = max(y);
    xoguess = x(I(1));
    sigmaguess = 1;
    guess = [Aguess xoguess sigmaguess];
  elseif(1)
    guess = [0 0 0];
  elseif(0)
    guess = [1.15 6.5 2];
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % set upper and lower bounds of search
  LB = [0 -inf 0];
  UB = [inf inf inf];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % set some fitting options
  options = optimset('TolX',1e-5,'MaxFunEvals',500,'Display','off');

  [fitparams1,ResNorm,Res,ExitFlag] =  ...
        lsqnonlin(@gauss_errfn,guess,LB,UB,options,x,y);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % also do fit by taking logarithm of both sides to make a linear problem:
  % log(y) = log(A*exp(-(x-xo)^2 / (2*sigma^2)))
  %	 = log(A) - (x-xo)^2 / (2*sigma^2);
  [fitparams1l,ResNorm,Res,ExitFlag] =  ...
        lsqnonlin(@log_gauss_errfn,guess,LB,UB,options,x,y);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  figure(1)
  hold on
  plot(x,fitparams1(1)*exp(-(x-fitparams1(2)).^2/(2*fitparams1(3)^2)),'b-');
  plot(x,fitparams1l(1)*exp(-(x-fitparams1l(2)).^2/(2*fitparams1l(3)^2)),'r-');

  disp(['lsqnonlin: ',num2str(fitparams1)])
end

if(1)
%else		% span parameter space "by hand"
  Arange = [0 : max(y)/20 : 1.1*max(y)];
  xorange = [-10 : 1.0 : 20];
  sigmarange = [0 : 0.5 : 20];

  [Agrid,xogrid,sigmagrid] = meshgrid(Arange,xorange,sigmarange);

  clear errgrid

  for n=1:length(Agrid(:));
    if (mod(n,1000)==0)
      disp([num2str(n),' of ',num2str(length(Agrid(:)))])
    end
    errgrid(n) = mean((gauss_errfn([Agrid(n) xogrid(n) sigmagrid(n)],x,y)).^2);
    %drawnow
  end

  [minerr, I] = min(errgrid);
  fitparams2 = [Agrid(I(1)) xogrid(I(1)) sigmagrid(I(1))];

  if(0)
    figure(3)
    plot(errgrid)
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  figure(1)
  hold on
  plot(x,fitparams2(1)*exp(-(x-fitparams2(2)).^2/(2*fitparams2(3)^2)),'k-');

  disp(['by hand: ',num2str(fitparams2)])

end

title({['Actual: ',num2str(actfitparams,2)]; ['lsqnonlin: ',num2str(fitparams1,2)];['by hand: ',num2str(fitparams2,2)]});

legend('Data','lsqnonlin','lsqnonlin w/ log','''by hand fit''');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
set(gca,'yscale','log')
%%
set(gca,'yscale','linear')

print -djpeg my_test_gaussfit_1