% % 
% % %%
% % set (parcluster, 'NumWorkers', 10)
% % matlabpool open 10
% %%
%%

parms.mup = (log(2)/22); 
parms.n = 5; 
tau = 29; 
parms.p0=385;  
parms.x1=log(2)/20; 
parms.x2=log(2)/35;
parms.a1 = parms.x1;
parms.a2 = parms.x2;
parms.p1=80; 
parms.murr=log(2)/1000;
parms.alphar=1;
parms.alphap=2;
parms.alpham=1;
parms.p2=100; 
parms.n2=5;
parms.n1=5;
parms.p3=280; 
parms.n3=5;

% simulate for 8 hours. Use as initial conds
parms.N = 1;
v=zeros(6,1);

v(1)=20;
v(2)=400;
v(3) = 0.00001;

v(4) = parms.N*v(1);
v(5) = parms.N*v(2);
v(6) = parms.N*v(3);

Tend=720;

timestep = 0.1;
tau = 29;

Nvec = [0.2, 1, 50];
for i = 1:length(Nvec);
    
    parms.N = Nvec(i);
    v(4) = parms.N*v(1);
    v(5) = parms.N*v(2);
    v(6) = parms.N*v(3);

    [y,dy]  = ChemicalLangevinEquationFinal( tau,v,Tend,timestep,parms );

    h = tau/timestep;
    y(:,[1:h])=[];
    dy(:,[1:h])=[];
        
end

m1=y(1, :);
p1=y(2, :);
r1=y(3, :);

t=timestep:timestep:Tend;

% use end points as initial condition starting points...
mIND = m1(end);
pIND = p1(end);
rIND = r1(end);

plot(t,m1);%,t,p1,t,r1)

%%

% tic
totalreps = 3000;
Nvec = [0.2, 1, 50];

v=zeros(6,1);
v(1)=mIND;
v(2)=pIND;
v(3) = rIND;

Tend=10000;
timestep = 0.1;

%%
% 50/50 divisions with stochastic dynamics
for i = 1:length(Nvec);
    i
    parms.N = Nvec(i);
    v(4) = parms.N*0.5*v(1);
    v(5) = parms.N*0.5*v(2);
    v(6) = parms.N*0.5*v(3);

    [ Dist] = CLEoptTIMINGwhile( tau,v,Tend,timestep,parms,totalreps );
    distTOT(i,:) = Dist;
    
end

distCLEunpert = distTOT - tau;
disp('done CLEunpert')

%%

for i = 1:length(Nvec);
    i
    parms.N = Nvec(i);
    v(4) = parms.N*v(1);
    v(5) = parms.N*v(2);
    v(6) = parms.N*v(3);

    [ Dist1, Dist2] = ChemicalLangevinDistBINOMIALdet( tau,v,Tend,timestep,parms,totalreps );
    distTOT1(i,:) = Dist1;
    distTOT2(i,:) = Dist2;
    
end

distCLEbinomial = distTOT1 - tau;
distCLEdet = distTOT2 - tau;
disp('done CLEuniform')
% toc

%%
% add pure deterministic
Tend = 5000;
parms.N = 1;
vdet = 0.5*v;
v(4) = parms.N*vdet(1);
v(5) = parms.N*vdet(2);
v(6) = parms.N*vdet(3);

[y,dy]  = ChemicalLangevinEquationFinal( tau,vdet,Tend,timestep,parms );

h = tau/timestep;
y(:,[1:h])=[];
dy(:,[1:h])=[];

mDet(1,:) = y(1,:);
pDet(1,:) = y(2,:);
rDet(1,:) = y(3,:);

t=timestep:timestep:Tend;

plot(t,pDet);


ix = find(pDet>20,1,'last');
distDet = t(ix);
distdet = distDet*ones(1,3000);

%%
% save distCLEbinomialDETlowsystem
% toc

%%

% calculate Kolmogorov distance from histogram

[h1,p1,ks1stat1] = kstest2(distCLEunpert(1,:),distCLEbinomial(1,:));
[h2,p2,ks1stat2] = kstest2(distCLEunpert(2,:),distCLEbinomial(2,:));
[h3,p3,ks1stat3] = kstest2(distCLEunpert(3,:),distCLEbinomial(3,:));


%%
% need 9 standard deviations
std1 = std(distCLEunpert(1,:));
std2 = std(distCLEunpert(2,:));
std3 = std(distCLEunpert(3,:));
std4 = std(distCLEbinomial(1,:));
std5 = std(distCLEbinomial(2,:));
std6 = std(distCLEbinomial(3,:));
std7 = std(distCLEdet(1,:));
std8 = std(distCLEdet(2,:));
std9 = std(distCLEdet(3,:));
distCLEuniform = distCLEbinomial;

%%
% plot standard deviations for diff system sizes

upper = 6000;
edges = linspace(0,upper,40);
figure()

subplot(3,4,3)
hist(distCLEunpert(3,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 2000])
a = xlim();
b = ylim();
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'D'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'c'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

str1 = '\sigma';
str2 = sprintf(' = %.0f',std3);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str3 = '50/50 division + stoch dynamics';
text((a(2)-a(1))/2,1.25*b(2),{str3},...
    'FontSize',9,'HorizontalAlignment','center', 'VerticalAlignment','bottom','color','k')

subplot(3,4,4)
hist(distCLEuniform(3,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 2000])
a = xlim();
b = ylim();
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'E'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'d'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str1 = '\sigma';
str2 = sprintf(' = %.0f',std6);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str3 = 'Binomial division + stoch dynamics';
text((a(2)-a(1))/2,1.25*b(2),{str3},...
    'FontSize',9,'HorizontalAlignment','center', 'VerticalAlignment','bottom','color','k')

subplot(3,4,2)
hist(distCLEdet(3,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 3200])
a = xlim();
b = ylim();
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'C'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'b'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str1 = '\sigma';
str2 = sprintf(' = %.2f',std9);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str3 = 'Binomial division + det dynamics';
text((a(2)-a(1))/2,1.25*b(2),{str3},...
    'FontSize',9,'HorizontalAlignment','center', 'VerticalAlignment','bottom','color','k')

subplot(3,4,7)
hist(distCLEunpert(2,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 1000])
a = xlim();
b = ylim();
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'g'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

str1 = '\sigma';
str2 = sprintf(' = %.0f',std2);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

subplot(3,4,8)
hist(distCLEuniform(2,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 1000])
a = xlim();
b = ylim();
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'h'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str1 = '\sigma';
str2 = sprintf(' = %.0f',std5);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

subplot(3,4,6)
hist(distCLEdet(2,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 3000])
a = xlim();
b = ylim();
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'f'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str1 = '\sigma';
str2 = sprintf(' = %.1f',std8);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')



subplot(3,4,11)
hist(distCLEunpert(1,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 1000])
a = xlim();
b = ylim();
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'k'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

str1 = '\sigma';
str2 = sprintf(' = %.0f',std1);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

subplot(3,4,12)
hist(distCLEuniform(1,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 1000])
a = xlim();
b = ylim();
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'l'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str = sprintf('K-S = %.2f',ks1stat1);
str1 = '\sigma';
str2 = sprintf(' = %.0f',std4);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

subplot(3,4,10)
hist(distCLEdet(1,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 1500])
a = xlim();
b = ylim();
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'j'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% str = sprintf('K-S = %.2f',ks2stat1);
% title('Random initial conditions');
str1 = '\sigma';
str2 = sprintf(' = %.0f',std7);
str = [str1, str2];
text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

subplot(3,4,1)
hist(distdet(1,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 3500])
a = xlim();
b = ylim();
text(a(1)-0.45*(a(2)-a(1)),b(2)/2,{'\Omega = 50'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'B'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'a'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% str = sprintf('K-S = %.2f',ks2stat1);
% title('Random initial conditions');
% str1 = '\sigma';
% str2 = sprintf(' = %.0f',std7);
% str = [str1, str2];
% text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
str3 = '50/50 division + det dynamics';
text((a(2)-a(1))/2,1.25*b(2),{str3},...
    'FontSize',9,'HorizontalAlignment','center', 'VerticalAlignment','bottom','color','k')

subplot(3,4,5)
hist(distdet(1,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 3500])
a = xlim();
b = ylim();
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'e'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
text(a(1)-0.45*(a(2)-a(1)),b(2)/2,{'\Omega = 1'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% str = sprintf('K-S = %.2f',ks2stat1);
% title('Random initial conditions');
% str1 = '\sigma';
% str2 = sprintf(' = %.0f',std7);
% str = [str1, str2];
% text(a(1)+0.4*(a(2)-a(1)),b(2)-0.2*(b(2)-b(1)),{str},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')


subplot(3,4,9)
hist(distdet(1,:),edges)
xlim([0 upper])
xlabel('Time (mins)')
ylabel('Frequency')
ylim([0 3500])
a = xlim();
b = ylim();
text(a(1)-0.45*(a(2)-a(1)),b(2)/2,{'\Omega = 0.2'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')
% text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'i'},...
%     'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

