%% Minimizing a function to find parameters (Basic MATLAB)

sz=8% Size of image
sigma=1.3; %Blob size

%These are the parameters that define the model
Photons=500;
bg=1;
x=sz/2+rand;
y=sz/2+rand;

%make noisy data
data=gaussianblob(newim([sz sz]),[x y],sigma,Photons);
data=data+bg;
data=noise(data,'poisson');

%setup for 'fminsearch'

%the first input to 'BlobMLE' must be a 1D vector, matlab variable type
minfunction=@BlobMLE;
options = optimset('Display','iter');

%setting an initial guess
X0(1)=sz/2;
X0(2)=sz/2;
X0(3)=max(data)*(2*pi*sigma^2);
X0(4)=min(data);

%fminsearch iterates to minimize 'fval'
tic;
[params,fval]=fminsearch(minfunction,X0,options,data,sigma); 
TimeElapsed=toc

fprintf('x: True:%g \t Found: %g\n',x,params(1))
fprintf('y: True:%g \t Found: %g\n',y,params(2))
fprintf('Photons: True:%g \t Found: %g\n',Photons,params(3))
fprintf('bg: True:%g \t Found: %g\n',bg,abs(params(4)))

%show data and model in color image
[fval model]=BlobMLE(params,data,sigma);
joinchannels('RGB',model,data)

%% The above is not fast.  Try the profiler to see what is slow.

profile on
%run above cell
profile viewer

%% Some comments about making things run fast in matlab:

%example of multiplying two vectors 
clear A B C 

N=1000;
A=rand(N,1);
B=rand(N,1);

%printing to the screen can be slow
tic 
for ii=1:N
   C(ii)=A(ii)*B(ii)
end
toc

clear C

tic 
for ii=1:N
   C(ii)=A(ii)*B(ii);
end
toc

% compare pre-allocate 'C' versus growin inside a loop
clear A B C 

N=100000;
A=rand(N,1);
B=rand(N,1);

%grow
tic 
for ii=1:N
   C(ii)=A(ii)*B(ii);
end
toc

%pre-allocate 
C=zeros(N,1);
tic 
for ii=1:N
   C(ii)=A(ii)*B(ii);
end
toc

%AVOID LOOPS!
tic 
C=A.*B;
toc

%% A slightly faster approach
minfunction=@BlobMLE_Fast;
options = optimset('Display','iter');

%setting an initial guess
X0(1)=sz/2;
X0(2)=sz/2;
X0(3)=max(data)*(2*pi*sigma^2);
X0(4)=min(data);

xmatrix=xx(data)-min(xx(data));
ymatrix=yy(data)-min(yy(data));

%fminsearch iterates to minimize 'fval'
tic;
[params,fval]=fminsearch(minfunction,X0,options,data,sigma,xmatrix,ymatrix); 
TimeElapsed=toc

fprintf('x: True:%g \t Found: %g\n',x,params(1))
fprintf('y: True:%g \t Found: %g\n',y,params(2))
fprintf('Photons: True:%g \t Found: %g\n',Photons,params(3))
fprintf('bg: True:%g \t Found: %g\n',bg,abs(params(4)))

%show data and model in color image
[fval model]=BlobMLE_Fast(params,data,sigma,xmatrix,ymatrix);
joinchannels('RGB',model,data)

%%

profile on
%run above cell
profile viewer


%% Convert to MATLAB variables

minfunction=@BlobMLE_Fast;
options = optimset('Display','iter');

%setting an initial guess
X0(1)=sz/2;
X0(2)=sz/2;
X0(3)=max(data)*(2*pi*sigma^2);
X0(4)=min(data);

xmatrix=xx(data)-min(xx(data));
ymatrix=yy(data)-min(yy(data));

xmatrix_M=single(xmatrix);
ymatrix_M=single(ymatrix);
data_M=single(data); 


%fminsearch iterates to minimize 'fval'
tic;
[params,fval]=fminsearch(minfunction,X0,options,data_M,sigma,xmatrix_M,ymatrix_M); 
TimeElapsed=toc

fprintf('x: True:%g \t Found: %g\n',x,params(1))
fprintf('y: True:%g \t Found: %g\n',y,params(2))
fprintf('Photons: True:%g \t Found: %g\n',Photons,params(3))
fprintf('bg: True:%g \t Found: %g\n',bg,abs(params(4)))

%show data and model in color image
[fval model]=BlobMLE_Fast(params,data,sigma,xmatrix,ymatrix);
joinchannels('RGB',model,data_M)








