function [stats] = sweepZn(m, inpOpts, simOpts, sweepOpts)
% function [stats] = sweepZn(m, inpOpts, simOpts, sweepOpts)
% 
% stats has 6 fields
% ymax1     -   first pulse maximum height
% ymin1     -   first pulse minimum height
% tmax1     -   first pulse maximum timing
% tmin1     -   first pulse minimum timing
% ySS       -   steady state value at long sim. time
% ySS_std   -   stdev of data for above's steady state (in case of limit cyc)
%
% At command line, run:
% 
% load simSweepKxKy0
% load modelNutlin
% 
% stats = sweepParams(m, inpOpts, simOpts, sweepOpts)


if ~exist('getPulseOpts', 'var')
    getPulseOpts = struct('tau', 1, 'threshold', 0.5, 'flooor', 8);
end


% **************************
% Important parameters
% **************************
tAve = 150;      % times beyond which the system has relaxed to its w-limit set (for averaging to 'steady-state')


% extend to make it possible to load either detailed or abstract models
m = m.update(m.p, m.k);
p0 = m.p;

% simulate model
[t, y] = simulate3y(m, inpOpts, simOpts);

h = figure;
hold on
plot(t{1}, y{1}, '--')
legend(m.yNames)

p1V = logspace(sweepOpts.LB, sweepOpts.UB, sweepOpts.N);
p2V = logspace(sweepOpts.LB, sweepOpts.UB, sweepOpts.N);

for i = 1:length(p1V)
    disp(sprintf('Simulating i = %d', i));
    for j = 1:length(p2V)
        % update model
        m.p(sweepOpts.p1) = p0(sweepOpts.p1)*p1V(i);
        m.p(sweepOpts.p2) = p0(sweepOpts.p2)*p1V(j);
        m = m.update(m.p, m.k);
        % simulate model
        [t, y] = simulate3y(m, inpOpts, simOpts);
        % get pulse stats
        [pP, pM, mP, mM] = getPulseStats(t{1}, y{1}, getPulseOpts);
        try
            % compute the pulse maxima
            ymax1{i,j} = pP.height;
            ymax2{i,j} = pM.height;
            % compute the pulse minima
            ymin1{i,j} = -mP.height;
            ymin2{i,j} = -mM.height;
        end
        try
            % compute the pulse maxima
            tmax1{i,j} = pP.time;
            tmax2{i,j} = pM.time;
            % compute the pulse minima
            tmin1{i,j} = -mP.time;
            tmin2{i,j} = -mM.time;
        end
        try
            % compute the eventual steady-state as well as the standard
            % deviation in case the trajectory is in a limit cycle
            iAve = find(t > tAve);
            ySS{i,j}     = mean(y(:, iAve),    2);
            ySS_std{i,j} =  std(y(:,iAve), [], 2);
        end
        figure(h);
        hold on
        plot(t{1}, y{1}, pP.time, pP.height, 'x', mP.time, mP.height, 'x')
        legend(m.yNames)
        hold off
        drawnow
%         pause
    end
end
save tmp
stats = struct('ymax1', ymax1, 'ymax2',   ymax2, ...
               'tmax1', tmax1, 'tmax2',   tmax2, ...
               'ySS',   ySS,   'ySS_std', ySS_std);