% Code to calculate the temperature and salinity dependence for
% the A. fundyense growth function
%
% Each data set was normalized by its maximum rate.  A cubic polynomial is 
% then fit to the data.  A linear extrapolation of the slope at 5 degrees
% is used for the region < 5 degrees to prevent abnormal curvature in the
% region poorly constrained by data.  The salinity dependence is also fit
% with a cubic polynomial.
%
% Charlie Stock
% cstock@alum.mit.edu
%

clear all

addpath Keafer;
% Keafer - GtCA29 (divisions/day)
load GTCA29.txt;
% convert from divisions per day
GTCA29(:,2) = log(2)*GTCA29(:,2);
% normalize
GTCA29(:,2) = GTCA29(:,2)./max(GTCA29(:,2));

% repeat steps for GTCA28 (also in divisions/day)
load GTCA28.txt;
% convert from divisions per day
GTCA28(:,2) = log(2)*GTCA28(:,2);
GTCA28(:,2) = GTCA28(:,2)./max(GTCA28(:,2));

% Watras et al. data (in days^-1).  There are two series of 
% temerature/growth parings. The first is at 32.5 ppt and the second is at
% 25.5 ppt.
addpath Watras;
load watras.txt;
% normalize 32.5 ppt time series
watras(:,2) = watras(:,2)./max(watras(:,2));
% normalize the 25.5 ppt time series
watras(:,3) = watras(:,3)./max(watras(:,3));

% Langdon data, divisions per day
addpath Langdon
load langdon_gt6.txt;
langdon = langdon_gt6;
langdon(:,2) = log(2)*langdon(:,2);
langdon(:,2) = langdon(:,2)./max(langdon(:,2));

% Gulf of Maine strains from Etheridge, BOF = Bay of Fundy
% Etheridge Data (BOF), data in day-1
addpath Etheridge
load BOF_temp.txt;
BOF = BOF_temp;
BOF(:,2) = BOF(:,2)./max(BOF(:,2));

% Etheridge Data (Monhegan Island)
load MI_temp.txt;
MI = MI_temp;
MI(:,2) = MI(:,2)./max(MI(:,2));

% put all growth data together
growth = [GTCA29(:,2)' GTCA28(:,2)' watras(:,2)' watras(:,3)' ...
          langdon(:,2)' BOF(:,2)' MI(:,2)'];

% put all temperature data together
temp = [GTCA29(:,1)' GTCA28(:,1)' watras(:,1)' watras(:,1)' ...
        langdon(:,1)' BOF(:,1)' MI(:,1)'];

% fit the data
[at, R] = polyfit(temp,growth,3);
x_gt5 = 5:0.1:27;
x_lt5 = [0:0.1:4.9];
% yt is the modeled fit to the temperature data before normalization
yt_gt5 = polyval(at,x_gt5);
% Rely on extrapolation of slope at 5 degrees to constrain cold region.
% This region is not well constrained by data, making polynomials behave
% unpredictably.
at_prime = 3*at(1)*5^2 + 2*at(2)*5 + at(3);
yt_lt5 = yt_gt5(1) - at_prime*[5-x_lt5];
yt = [yt_lt5 yt_gt5];
xt = [x_lt5 x_gt5];
% Calculate R2 for the curve before normalization - 2 points below 5 deg. C
yt_R2(1) = yt_gt5(1) - at_prime*(5-temp(1));
yt_R2(2) = yt_gt5(1) - at_prime*(5-temp(2));
yt_R2(3:size(temp,2)) = polyval(at,temp(3:end));
R2_t = 1-(sum((growth-yt_R2).^2))/(sum((growth-mean(growth)).^2));

% Normalize the polynomial fit to yield a maximum value of 1
ymax = max(yt);
atn = at/ymax;
% evaluate the normalized polynomial
ytn_gt5 = polyval(atn,x_gt5);
% calculate the first derivative of the relationship at 5 degrees
atn_prime = 3*atn(1)*5^2 + 2*atn(2)*5 + atn(3);
ytn_lt5 = ytn_gt5(1) - atn_prime*[5-x_lt5];
ytn = [ytn_lt5 ytn_gt5];
% Calculate R2 for the normalized curve
ytn_R2(1) = ytn_gt5(1) - atn_prime*(5-temp(1));
ytn_R2(2) = ytn_gt5(1) - atn_prime*(5-temp(2));
ytn_R2(3:size(temp,2)) = polyval(atn,temp(3:end));
R2_tn = 1-(sum((growth-ytn_R2).^2))/(sum((growth-mean(growth)).^2));

% Plot Temperature Relationship
figure(1)
clf
hold on
plot(GTCA29(:,1),GTCA29(:,2),'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(GTCA28(:,1),GTCA28(:,2),'ko','MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
plot(watras(:,1),watras(:,2),'ks','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(watras(:,1),watras(:,3),'ks','MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
plot(langdon(:,1),langdon(:,2),'kv','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(BOF(:,1),BOF(:,2),'kv','MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
plot(MI(:,1),MI(:,2),'kd','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
% Before Normalization
plot(xt,yt,'-','LineWidth',2,'Color',[0.5 0.5 0.5]);
% After Normalization
plot(xt,ytn,'k-','LineWidth',2);
ylabel('growth');
xlabel('temperature');
legend('GTCA29','GTCA28','GTMP(1)','GTMP(2)','GT6','BOF','MI','least-squares fit', ...
       'normalized')
title('Temperature Function')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constructing the salinity relationship
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Prakash
addpath Prakash
load prakash.txt;
prakash(:,2) = prakash(:,2)/max(prakash(:,2));
load BOF_salt.txt;
BOF_salt(:,2) = BOF_salt(:,2)/max(BOF_salt(:,2));
load MI_salt.txt;
MI_salt(:,2) = MI_salt(:,2)/max(MI_salt(:,2));
salt = [prakash(:,1)' BOF_salt(:,1)' MI_salt(:,1)'];
grow_salt = [prakash(:,2)' BOF_salt(:,2)' MI_salt(:,2)'];

[as] = polyfit(salt,grow_salt,3);
xs = [1:1:40];
% un-normalized relationship
ys = polyval(as,xs);
% calculate R2 for the unnormalized relationship
ys_R2 = polyval(as,salt);
R2_s = 1-sum((grow_salt-ys_R2).^2)/sum((grow_salt-mean(grow_salt)).^2);
% normalize to get a maximum of 1
ysmax = max(ys);
asn = as/ysmax;
ysn = polyval(asn,xs);
% do R2 calculation for the normalized curve
ysn_R2 = polyval(asn,salt);
R2_sn = 1-sum((grow_salt-ysn_R2).^2)/sum((grow_salt-mean(grow_salt)).^2);

figure(2)
clf
plot(prakash(:,1),prakash(:,2),'ko');
hold on
plot(BOF_salt(:,1),BOF_salt(:,2),'s','MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
plot(MI_salt(:,1),MI_salt(:,2),'k.','MarkerSize',14);
plot(xs,ys,'-','LineWidth',2,'Color',[0.5 0.5 0.5]);
plot(xs,ysn,'k-','LineWidth',2)
legend('Prakash','Etheridge, BOF','Etheridge, MI','least-squares fit', ...
       'normalized');
ylabel('growth');
xlabel('salt');

% Plot growth as a function of temperature and salinity
figure(3)
clf
[T S] = meshgrid([0:0.2:27],[5:0.25:40]);
T_FAC = atn(1)*T.^3 + atn(2)*T.^2 + atn(3)*T + atn(4);
T5 = atn(1)*5.^3 + atn(2)*5.^2 + atn(3)*5 + atn(4);
ee = find(T < 5);
T_FAC(ee) = T5 - atn_prime*(5-T(ee));
S_FAC = asn(1)*S.^3 + asn(2)*S.^2 + asn(3)*S + asn(4);
G = 0.58.*T_FAC.*S_FAC;

% sizes used to plot data points in T,S space
a = size(GTCA29,1);
b = size(GTCA28,1);
c = size(watras,1);
d = size(langdon,1);
e = size(BOF,1);
f = size(MI,1);
aa= size(prakash,1);
bb = size(BOF_salt,1);
cc = size(MI_salt,1);

% Marker size
ms = 4;

% Growth as a function of temperature and salinity
figure(3)
clf
[cs,h] = contour(T,S,G,[-0.2:0.1:0.60],'k-')
hold on
clabel(cs,h,[-0.2:0.1:0.6],'fontsize',7,'LabelSpacing',200);
xlabel('Temp., deg. C');
ylabel('Salinity (PSU)');

% plot data points (temp)
plot3(GTCA29(:,1),ones(a,1)*31,ones(a,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(GTCA28(:,1),ones(b,1)*31,ones(b,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(watras(:,1),ones(c,1)*25.5,ones(c,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(watras(:,1),ones(c,1)*32.5,ones(c,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(langdon(:,1),ones(d,1)*29.1,ones(d,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(BOF(:,1),ones(e,1)*30,ones(e,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(MI(:,1),ones(f,1)*30,ones(f,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
% salt data
plot3(ones(aa,1)*10,prakash(:,1),ones(aa,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(ones(bb,1)*20,BOF_salt(:,1),ones(bb,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
plot3(ones(cc,1)*20,MI_salt(:,1),ones(cc,1)*100,'ko','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',ms);
