%
% This tests the use of multilevel Monte Carlo for a
% continuous-time Markov process with adaptive time-stepping
%
% Developed with Chris Lester, summer 2014
%

function ctmp_adapted

close all; clear all;

%
% call mlmc_test to perform MLMC tests
%

nvert = 2;      % plotting option (1 for slides, 2 for papers, 3 for full set)
M     = 4;      % refinement cost factor (2^gamma in general MLMC Thm)

N     = 1000;   % samples for convergence tests
L     = 6;      % levels for convergence tests 

N0    = 20;     % initial number of samples on coarsest levels
Eps   = [ 1 2 5 10 20 ];  % desired accuracies for MLMC calcs

figs = mlmc_test(@sde_l, M, N,L, N0,Eps, nvert);

subplot(nvert,2,2*nvert-1)
legend('\epsilon=1','\epsilon=2','\epsilon=5','\epsilon=10','\epsilon=20',1)

subplot(nvert,2,2*nvert)
axis([1 20 2e5 1e8])

print('-deps2','ctmp_adapted.eps')


%-------------------------------------------------------
%
% level l estimator for continuous-time Markov process with factor 4 refinement
%

function [sum1 sum2] = sde_l(l,N)

coeffs = [  1  0.04 .002  0.5 ];  % propensity coefficients
nu     = [ -1    0   -2    2; ...  % stoichiometric matrix
            0   -1    1   -1; ...
            0    1    0    0 ];

M = 4;    % refinement factor in adaptive timesteps on each level
T = 30;   % simulation time interval

sum1(1:4) = 0;
sum2(1:2) = 0;

for N1 = 1:10000:N
  N2 = min(10000,N-N1+1);

  t  = zeros(1,N2);
  tf = t;
  tc = t;

  qf = repmat( [10^5; 0; 0], 1,N2);
  qc = qf;

  lamf = zeros(4,N2);
  lamc = lamf;

  if l==0

    while any(t<T)
      tm = t;
      t  = tf;

      p1 = poissrnd(lamf.*repmat(t-tm,4,1));

      qf = qf + nu*p1;

      qf   = max(qf,0);
      lamf = propensity(qf,coeffs);
      tf   = min(tf+dt(qf,lamf,nu,M,l), T);
    end

  else

    while any(t<T)
      tm = t;
      t  = min(tf,tc);

      lam_min = min(lamf,lamc);
      lam_dif = abs(lamf-lamc);

      p1 = poissrnd(lam_min.*repmat(t-tm,4,1));
      p2 = poissrnd(lam_dif.*repmat(t-tm,4,1));

      qf = qf + nu*(p1 + (lamf>lamc).*p2);
      qc = qc + nu*(p1 + (lamc>lamf).*p2);

      ind = find(tf==t);
      if size(ind)>0
        qf(:,ind)   = max(qf(:,ind),0);
        lamf(:,ind) = propensity(qf(:,ind),coeffs);
        tf(ind)     = min(tf(ind)+dt(qf(:,ind),lamf(:,ind),nu,M,l), T);
      end

      ind = find(tc==t);
      if size(ind)>0
        qc(:,ind)   = max(qc(:,ind),0);
        lamc(:,ind) = propensity(qc(:,ind),coeffs);
        tc(ind)     = min(tc(ind)+dt(qc(:,ind),lamc(:,ind),nu,M,l-1), T);
      end
    end
  end

  Pf = qf(3,:);
  if (l>0)
    Pc = qc(3,:);
  else
    Pc = zeros(1,N2);
  end

  sum1(1) = sum1(1) + sum(Pf-Pc);
  sum1(2) = sum1(2) + sum((Pf-Pc).^2);
  sum1(3) = sum1(3) + sum((Pf-Pc).^3);
  sum1(4) = sum1(4) + sum((Pf-Pc).^4);
  sum2(1) = sum2(1) + sum(Pf);
  sum2(2) = sum2(2) + sum(Pf.^2);
end

%
%------ propensity function ---------
%

function lam = propensity(q,coeffs)

lam = [ coeffs(1)*q(1,:); ...
        coeffs(2)*q(2,:); ...
        coeffs(3)*q(1,:).*(q(1,:)-1); ...
        coeffs(4)*q(2,:) ];

%
%------ adaptive timestep ---------
%

function h = dt(q,lam,nu,M,l)

h = 0.25 * M^(-l) * min( (q+1) ./ abs(nu*lam) );






