clear; clc;

colors = colormap(jet);
colors_to_use = colors([1,14,26,42,50,62],:);
temps = [17; 20; 23; 25; 27; 31];

%% read in the control data
[num words] = xlsread('DidiniumPopCyclesTemps_final_14_02_2020','Population densities per mL');

pc_17_c = num(38:41,[4 6]); 
    pc_r_17 = log(pc_17_c(:,2)./pc_17_c(:,1))./2;
    pc_r_17_mean = mean(pc_r_17);
    pc_r_17_std = std(pc_r_17);
pc_20_c = num(42:45,[4 6]); 
    pc_r_20 = log(pc_20_c(:,2)./pc_20_c(:,1))./2;
    pc_r_20_mean = mean(pc_r_20);
    pc_r_20_std = std(pc_r_20);
pc_23_c = num(46:49,[4 6]); 
    pc_r_23 = log(pc_23_c(:,2)./pc_23_c(:,1))./2;
    pc_r_23_mean = mean(pc_r_23);
    pc_r_23_std = std(pc_r_23);
pc_25_c = num(50:53,[4 6]); 
    pc_r_25 = log(pc_25_c(:,2)./pc_25_c(:,1))./2;
    pc_r_25_mean = mean(pc_r_25);
    pc_r_25_std = std(pc_r_25);
pc_27_c = num(54:57,[4 6]); 
    pc_r_27 = log(pc_27_c(:,2)./pc_27_c(:,1))./2;
    pc_r_27_mean = mean(pc_r_27);
    pc_r_27_std = std(pc_r_27);
pc_31_c = num(58:59,[4 6]); 
    pc_r_31 = log(pc_31_c(:,2)./pc_31_c(:,1))./2;
    pc_r_31_mean = mean(pc_r_31);
    pc_r_31_std = std(pc_r_31);

pc_r_mean = [pc_r_17_mean; pc_r_20_mean; pc_r_23_mean; pc_r_25_mean; pc_r_27_mean; pc_r_31_mean];
pc_r_std = [pc_r_17_std; pc_r_20_std; pc_r_23_std; pc_r_25_std; pc_r_27_std; pc_r_31_std];

figure(1);clf(1); % not used
for i = 1:6
    box on; hold on;    
    plot([temps(i) temps(i)],[pc_r_mean(i)-pc_r_std(i) pc_r_mean(i)+pc_r_std(i)],'-','Color',colors_to_use(i,:));
    plot(temps(i),pc_r_mean(i),'ok','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
    ylabel('\itr','FontSize',12);
    xlim([16 32]);
end
%% read in averaged count (per mL) data

days = 0:18;
columns_to_read_para = [2 4 6 8 10 12];
columns_to_read_did = columns_to_read_para+1;
num_reps = [5 6 3 5 6 3];

[num words] = xlsread('DidiniumPopCyclesTemps_final_14_02_2020','Combined for plotting'); % read in appropriate data

pb_avg = num(1:19,columns_to_read_para); % pull out paramecium means
pb_std = num(23:41,columns_to_read_para); % pull out paramecium stds
did_avg = num(1:19,columns_to_read_did); % pull out didinium means
did_std = num(23:41,columns_to_read_did); % pull out didinium stds

for i = 1:6
    pb_ste(:,i) = pb_std(:,i)./num_reps(i); % calculate std across replicates
    did_ste(:,i) = did_std(:,i)./num_reps(i); % calculate std across replicates
end

%% plot raw densities through time and in phase space
% produces ms figure #1

figure(1);clf(1);
for i = 1:6
    subplot(221);
    box on; hold on;
        plot(days(~isnan(pb_avg(:,i))),pb_avg(~isnan(pb_avg(:,i)),i),'-o','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
        ylabel('{\itParamecium} mL^{-1}','FontSize',12);

    subplot(223);
    box on; hold on;
        plot(days(~isnan(did_avg(:,i))),did_avg(~isnan(did_avg(:,i)),i),'-o','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
        ylabel('{\itDidinium} mL^{-1}','FontSize',12);
        xlabel('Time (days)','FontSize',12);

    subplot(122);
    box on; hold on;
        plot(pb_avg(~isnan(pb_avg(:,i)),i),did_avg(~isnan(did_avg(:,i)),i),'-o','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
        ylabel('{\itDidinium} mL^{-1}','FontSize',12);
        xlabel('{\itParamecium} mL^{-1}','FontSize',12);
end
legend(num2str(temps));
shg;

gtext('A','FontSize',12);
gtext('B','FontSize',12);
gtext('C','FontSize',12);

%% read in fitted parameters
[num words] = xlsread('Fitted parameters_14_Feb_20','Compiled parameter sets'); 
rows_for_parameters = [1 11 21 31 41 51]; % r is on row 1, rslope on row 2, etc.

for i = 1:6
    r_avg(i) = num(rows_for_parameters(i),3); r_cis(i,1:2) = num(rows_for_parameters(i),4:5);    
    rslope_avg(i) = num(rows_for_parameters(i)+1,3); rslope_cis(i,1:2) = num(rows_for_parameters(i)+1,4:5);
    a_avg(i) = num(rows_for_parameters(i)+2,3); a_cis(i,1:2) = num(rows_for_parameters(i)+2,4:5);
    m_avg(i) = -num(rows_for_parameters(i)+3,3); m_cis(i,1:2) = -num(rows_for_parameters(i)+3,4:5);
    h_avg(i) = num(rows_for_parameters(i)+4,3); h_cis(i,1:2) = num(rows_for_parameters(i)+4,4:5);
    e_avg(i) = num(rows_for_parameters(i)+5,3); e_cis(i,1:2) = num(rows_for_parameters(i)+5,4:5);
    d_avg(i) = num(rows_for_parameters(i)+6,3); d_cis(i,1:2) = num(rows_for_parameters(i)+6,4:5);
    Cd_avg(i) = num(rows_for_parameters(i)+7,3); Cd_cis(i,1:2) = num(rows_for_parameters(i)+7,4:5);
end

%% starting values for ODE
y0 = [8.33 0.33]; % initial prey and predator densities
tspan = [0 18]; % start end times

% produces ms figure 2
figure(2);clf(2);

for i = 1:6
ode = @(t,y) MR_HVH_model(t,y,r_avg(i),rslope_avg(i),a_avg(i),m_avg(i),h_avg(i),e_avg(i),d_avg(i),1,Cd_avg(i)); % compile function and call
[t1,y1] = ode45(ode, [0:0.25:18], y0); % return time and population density vectors  

if i > 1
    y1(min(find(y1(:,1) < 0.6 & t1 > 3)):end,1) = 0;
    y1(min(find(y1(:,2) < 0.6 & t1 > 3)):end,2) = 0;
end

subplot(2,6,i);
box on; hold on;
    jbfill(days(~isnan(pb_avg(:,i))),pb_avg(~isnan(pb_avg(:,i)),i)'+pb_ste(~isnan(pb_avg(:,i)),i)',pb_avg(~isnan(pb_avg(:,i)),i)'-pb_ste(~isnan(pb_avg(:,i)),i)',colors_to_use(i,:).*0.8,'w',1,0.2); hold on;
    plot(days,pb_avg(:,i),'o','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
    plot(t1,y1(:,1),'-','Color',colors_to_use(i,:),'LineWidth',2);
    title(num2str(temps(i)));
    ylim([0 150]);

subplot(2,6,i+6);
box on; hold on; 
    jbfill(days(~isnan(did_avg(:,i))),did_avg(~isnan(did_avg(:,i)),i)'+did_ste(~isnan(did_avg(:,i)),i)',did_avg(~isnan(did_avg(:,i)),i)'-did_ste(~isnan(did_avg(:,i)),i)',colors_to_use(i,:).*0.8,'w',1,0.2); hold on;
    plot(days,did_avg(:,i),'o','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
    plot(t1,y1(:,2),'-','Color',colors_to_use(i,:),'LineWidth',2);
    ylim([0 30]);
end
subplot(2,6,1);
    ylabel('{\itParamecium} mL^{-1}','FontSize',12);
subplot(2,6,7);
    ylabel('{\itDidinium} mL^{-1}','FontSize',12);
subplot(2,6,10);
    xlabel('Time (days)','FontSize',12);
shg;

gtext('A','FontSize',12);
gtext('B','FontSize',12);
gtext('C','FontSize',12);
gtext('D','FontSize',12);
gtext('E','FontSize',12);
gtext('F','FontSize',12);
gtext('G','FontSize',12);
gtext('H','FontSize',12);
gtext('I','FontSize',12);
gtext('J','FontSize',12);
gtext('K','FontSize',12);
gtext('L','FontSize',12);

%% plot parameters

param_means = [r_avg' rslope_avg' a_avg' m_avg' h_avg' e_avg' d_avg' Cd_avg'];
param_cis = cat(3,r_cis,rslope_cis,a_cis,m_cis,h_cis,e_cis,d_cis,Cd_cis);
paramswithunits = {'{\itr} (day^{-1})' '{\itr}_{slope} (mL prey^{-1} day^{-1} )' '{\ita} (mL pred^{-1} day^{-1})',...
    '{\itm} (-)' '{\ith} (day)' '{\ite} (pred prey^{-1})' '{\itd} (day^{-1})' '{\itCd} (prey mL^{-1})'};

param_means_full = param_means;
param_means_full(1,:) = [1.03 0.0078 0.16 0.022 0.045 0.06 0.47 0.009];
param_means_full(2,:) = [1.60 0.01 0.37 0.12 0.10 0.07 0.55 0.00];
param_means_full(6,:) = [2.36 0.0012 1.22 0.0023 0.0073 0.048 3.84 2.67E-07];

% produces ms figure # 3
figure(3); clf(3);

for j = 1:8 % sequence of parameters
for i = 1:6 % sequence of temperatures
subplot(4,2,j);
    box on; hold on;
    plot([temps(i) temps(i)],[param_cis(i,1,j) param_cis(i,2,j)],'-','Color',colors_to_use(i,:));
    plot(temps(i),param_means(i,j),'ok','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
    ylabel(paramswithunits{j},'FontSize',12);
    xlim([16 32]);
    set(gca,'XTick',[17 20 23 25 27 31]);
end
end
subplot(4,2,7);
    xlabel('Temperature ({^o}C)','FontSize',12);
subplot(4,2,8);
    ylim([1e-7 1e-0]);
    set(gca,'YTick',[1e-6 1e-1]);

subplot(4,2,1);
hold on;
for i = 1:6
    plot([temps(i) temps(i)],[pc_r_mean(i)-pc_r_std(i) pc_r_mean(i)+pc_r_std(i)],'-','Color',colors_to_use(i,:));
    plot(temps(i),pc_r_mean(i),'o','Color',colors_to_use(i,:),'MarkerFaceColor','w');
end

shg;

gtext('A','FontSize',12);
gtext('B','FontSize',12);
gtext('C','FontSize',12);
gtext('D','FontSize',12);
gtext('E','FontSize',12);
gtext('F','FontSize',12);
gtext('G','FontSize',12);
gtext('H','FontSize',12);

ylim([0 0.3]);
    
%% assess individual parameter effects on dynamics

mean_paramter_set = mean(param_means); % take mean across all means
titles = {'\itr' '{\itr}_{slope}' '{\ita}' '{\itm}' '{\ith}' '{\ite}' '{\itd}' '{\itCd}'}; % parameter titles

% solve for dynamics of mean parameter set
ode = @(t,y) MR_HVH_model(t,y,mean(r_avg),mean(rslope_avg),mean(a_avg),mean(m_avg),mean(h_avg),mean(e_avg),mean(d_avg),1,mean(Cd_avg)); % compile function and call
[t1,y1] = ode45(ode, [0:0.25:18], y0); % return time and population density vectors  
    y1(min(find(y1(:,1) < 0.6 & t1 > 3)):end,1) = 0;
    y1(min(find(y1(:,2) < 0.6 & t1 > 3)):end,2) = 0;

% produces ms figure # 4    
figure(5);clf(5);
hold on; box on;
for i = 1:8 % run through each parameter
    mean_paramter_set = mean(param_means); % take mean across all means
    mean_paramter_set(i) = max(param_means(:,i)); % swap out max for mean of parameter i
    ode = @(t,y) MR_HVH_model(t,y,mean_paramter_set(1),mean_paramter_set(2),mean_paramter_set(3),...
        mean_paramter_set(4),mean_paramter_set(5),mean_paramter_set(6),mean_paramter_set(7),1,mean_paramter_set(8)); % compile function and call
        [t2,y2] = ode45(ode, [0:0.25:18], y0); % return time and population density vectors  
        y2(min(find(y2(:,1) < 0.6 & t1 > 3)):end,1) = 0;
        y2(min(find(y2(:,2) < 0.6 & t1 > 3)):end,2) = 0;
    
    mean_paramter_set(i) = min(param_means(:,i)); % swap out min for mean of parameter i
    ode = @(t,y) MR_HVH_model(t,y,mean_paramter_set(1),mean_paramter_set(2),mean_paramter_set(3),...
        mean_paramter_set(4),mean_paramter_set(5),mean_paramter_set(6),mean_paramter_set(7),1,mean_paramter_set(8)); % compile function and call
        [t3,y3] = ode45(ode, [0:0.25:18], y0); % return time and population density vectors  
        y3(min(find(y3(:,1) < 0.6 & t1 > 3)):end,1) = 0;
        y3(min(find(y3(:,2) < 0.6 & t1 > 3)):end,2) = 0;

    subplot(4,2,i);
    box on; hold on;
        h1 = plot(t1,y1(:,1),'-','Color',[0.5 0.5 0.5],'LineWidth',2);
        h2 =plot(t1,y1(:,2),'--','Color',[0.5 0.5 0.5],'LineWidth',2);
        h3 = plot(t2,y2(:,1),'-','Color',colors_to_use(4,:),'LineWidth',2);
        h4 = plot(t2,y2(:,2),'--','Color',colors_to_use(4,:),'LineWidth',2);
        h5 = plot(t3,y3(:,1),'-','Color',colors_to_use(6,:),'LineWidth',2);
        h6 = plot(t3,y3(:,2),'--','Color',colors_to_use(6,:),'LineWidth',2);
        title(titles(i));
            xlim([0 18]);
        if i > 6
            ylim([0 20]);
        else
            ylim([0 200]);
        end
end
legend([h1 h2 h3 h4 h5 h6],'Mean-prey','Mean-predator','Max-prey','Max-predator','Min-prey','Min-predator');
    xlabel('Time (days)','FontSize',12);
    ylabel('Population density (cells mL^{-1})','FontSize',12);

gtext('A','FontSize',12);
gtext('B','FontSize',12);
gtext('C','FontSize',12);
gtext('D','FontSize',12);
gtext('E','FontSize',12);
gtext('F','FontSize',12);
gtext('G','FontSize',12);
gtext('H','FontSize',12);

%% read in the cell size data
[num words] = xlsread('DidiniumPopCyclesTemps','Sheet3');

temp_sizes = num(10:202,1);
cell_volumes = num(10:202,7);
[mean_cell_volume std_cell_volume] = grpstats(cell_volumes,temp_sizes,{'mean','std'});

cell_volumes_start = num(1:19,7);
cell_volumes_start_mean = mean(cell_volumes_start);

% produces ms figure # 5
figure(6); clf(6);
box on; hold on;

for i = 1:6    
    plot(temps(i),mean_cell_volume(i),'ok','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
    plot([temps(i) temps(i)],...
        [mean_cell_volume(i)-std_cell_volume(i) mean_cell_volume(i)+std_cell_volume(i)],'-','Color',colors_to_use(i,:),'MarkerFaceColor',colors_to_use(i,:));
end

xlabel('Temperature ({^o}C)','FontSize',12);
ylabel('Cell volume at peak (\mu m^3)','FontSize',12);
shg;