% eda11_07
% example of a grid search that uses a table
% lookup of the function
% f(z) = sin(pi*sin(pi*z*sin(pi*(1-z^2))))

clear all;

% true values of the model parameters
mtrue = [0.5142, 0.7524]';

% the data are a function of time
N = 10001;
tmin = 0;
tmax = 1;
Dt = (tmax-tmin)/(N-1);
t = tmin+Dt*[0:N-1]';

% subsamples indices for plotting purposes
Nsub = 20;
isub = [1:Nsub:N]';

% the true data obey the formula
% dtrue = f(z) with z = m1 + m2*t
z = mtrue(1)+mtrue(2)*t;
dtrue = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
% synthetic estimated data is true data plus random noise
sigmad = 0.1;
dobs = dtrue + random('Normal',0,sigmad,N,1);
etrue = dobs - dtrue;
Etrue = etrue'*etrue;
fprintf('true:     m %f %f    E %f\n', mtrue(1), mtrue(2), Etrue );

% plot the true and observed data
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([tmin,tmax,-1.2,1.2]);
plot(t(isub),dobs(isub),'k.','LineWidth',2);
plot(t(isub),dtrue(isub),'-','LineWidth',7,'Color',[0.9,0.9,0.9]);
xlabel('t');
ylabel('d');

% range of m1 for grid search
M1 = 301;
M1=101;
m1min = 0;
m1max = 1;
Dm1 = (m1max-m1min)/(M1-1);
m1 = m1min+Dm1*[0:M1-1]';

% range of m2 for grid search
M2 = 301;
M2=101;
m2min = 0;
m2max = 1;
Dm2 = (m2max-m2min)/(M1-1);
m2 = m2min+Dm2*[0:M2-1]';

% standard grid search that computes f(z) exactly
% mote use of tic/toc pair to time exectiion
tic;
first = 1;
for m1trial = m1'
for m2trial = m2'
    z = m1trial+m2trial*t;
    dtrial = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
    etrial = dobs - dtrial;
    Etrial = etrial'*etrial;
    if( first==1 )
        Ebest1 = Etrial;
        m1best1 = m1trial;
        m2best1 = m2trial;
        first=0;
    elseif (Etrial < Ebest1 )
        Ebest1 = Etrial;
        m1best1 = m1trial;
        m2best1 = m2trial;
    end
end
end
duration1=toc;

% print and plot results
fprintf('exact:    m %f %f    E %f Dt %f\n', m1best1, m2best1, Ebest1, duration1 );
z = m1best1+m2best1*t;
dest = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
plot(t(isub),dest(isub),'-','LineWidth',4,'Color',[0.6,0.6,0.6]);

% grid search with table lookup of the funtion f(z)
tic;

% compute the table
K = 1001;
zmintable = -0.1;
zmaxtable =  2.1;
Dztable = (zmaxtable-zmintable)/(K-1);
ztable = zmintable+Dztable*[0:K-1]';
dtable=sin(pi*sin(pi*ztable.*sin(pi*(1-ztable.^2))));

% grid search with table lookup
first = 1;
for m1trial = m1'
for m2trial = m2'

    z = m1trial+m2trial*t;
    n = floor((z-zmintable)/Dztable)+1; % table index
    dtrial = dtable(n); % table lookup
    etrial = dobs - dtrial;
    Etrial = etrial'*etrial;
    if( first==1 )
        Ebest2 = Etrial;
        m1best2 = m1trial;
        m2best2 = m2trial;
        first=0;
    elseif (Etrial < Ebest2 )
        Ebest2 = Etrial;
        m1best2 = m1trial;
        m2best2 = m2trial;
    end
    
end
end
duration2=toc;

% print and plot estimated solution
fprintf('lookup:   m %f %f    E %f Dt %f\n', m1best2, m2best2, Ebest2, duration2 );
z = m1best2+m2best2*t;
dest = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
plot(t(isub),dest(isub),'k-','LineWidth',2);

