% endo_plot_all.m
% This program reads in the endogenous clock data and plots up the first
% cycle from each data set on a common plot for comparison.  Detailed
% descriptions of the rationale for the various steps can be found in the
% posted detailed description of the germination function.
%
% Charlie Stock
% cstock@alum.mit.edu
%

clf
clear all

% This routine uses several functions within the "timeplt" toolbox from the
% USGS (Rich Signell, Chuck Denham, 
addpath timeplt

% read Matrai Data, note that textread is a discontinued matlab
% function (replaced by textscanf).  I have included a copy in this
% directory.  "mat_date" is short for Matrai date.
[mat_date,st1,st2,st3] = textread('Matrai_NaN.txt','%s %d %d %d');
mat_date = char(mat_date);
mat_date = mat_date;

% loop through each date and isolate the day and the month based on the
% position of the /, the julian function is from Rich Signell's toolbox and
% is eventually used for plotting.
for i = 1:size(mat_date,1);
    aa = find(mat_date(i,:) == '/');
    mon_mat(i) = str2num(mat_date(i,1:(aa(1)-1)));
    day_mat(i) = str2num(mat_date(i,(aa(1)+1):(aa(2)-1)));
    year_mat(i) = str2num(mat_date(i,(aa(2)+1:end)));
    jday_mat(i) = julian([year_mat(i),mon_mat(i),day_mat(i),12,0,0]);
end

% read 84-88 data from Anderson lab, each of three isolates, 
% repeat extraction of date information as with Matrai
[old_date,iso1,iso2,iso3] = textread('clock84_88.txt','%s %d %d %d');
old_date = char(old_date);
for i = 1:size(old_date,1);
    aa = find(old_date(i,:) == '/');
    mon_old(i) = str2num(old_date(i,1:(aa(1)-1)));
    day_old(i) = str2num(old_date(i,(aa(1)+1):(aa(2)-1)));
    year_old(i) = str2num(old_date(i,(aa(2)+1:end)));
    jday_old(i) = julian([year_old(i),mon_old(i),day_old(i),12,0,0]);
end

% read in 1998 data from Amy Bronzino, 1 isolate
[amy_date,amy,dum] = textread('clock_1998.txt','%s %f %f');
amy_date = char(amy_date);

for i = 1:size(amy_date,1);
    aa = find(amy_date(i,:) == '/');
    mon_amy(i) = str2num(amy_date(i,1:(aa(1)-1)));
    day_amy(i) = str2num(amy_date(i,(aa(1)+1):(aa(2)-1)));
    year_amy(i) = str2num(amy_date(i,(aa(2)+1:end)));
    jday_amy(i) = julian([year_amy(i)+1900,mon_amy(i),day_amy(i),12,0,0]);
end

% plot each isolate separately and put on a julian day axis.  Then color in
% dots from the first year of isolation and before second rise in
% germination potential.  There are a few other details regarding early
% points (see detailed germ. function description, particularly the 
% captions of Figs. 9-11)

title('First Cycle Endogenous Clock Data')
hold on

% Figure 10 of appendix A

% Matrai, station 1
figure(1)
subplot(3,1,1)
% trick the legend
plot(0,0,'Marker','none','LineStyle','none');
hold on
aa = find(isfinite(st1));
% plot all stations w/finite values
plot(jday_mat(aa),st1(aa),'Color','k','Marker','o','MarkerSize',6)
% first 22 finite stations are in the first year
plot(jday_mat(aa(1:22)),st1(aa(1:22)),'Color','k','Marker','o', ...
     'MarkerSize',6,'MarkerFaceColor',[0 0 0]);
jday_min = julian([1999 11 1 0 0 0]);
jday_max = julian([2004 11 1 0 0 0]);
jd = jday_min:1:jday_max;
% slightly modified version of Signell's gregaxm routine to get rid of some
% pesky labeling issues.
gregaxm_mod(jd,3);
axis([jday_min jday_max 0 100]);
ylabel('germination potential');
legend('station 61')
legend boxoff

% Repeat general steps from above for all stations

% Matrai, station 2
subplot(3,1,2)
aa = find(isfinite(st2));
% trick the legend
plot(0,0,'Marker','none','LineStyle','none');
hold on
plot(jday_mat(aa),st2(aa),'Color','k','Marker','o','MarkerSize',6)
% First 22 points are within the first year and before the next rise
% in germination.
plot(jday_mat(aa(1:22)),st2(aa(1:22)),'Color','k','Marker','o', ...
     'MarkerSize',6,'MarkerFaceColor',[0 0 0]);
jday_min = julian([1999 11 1 0 0 0]);
jday_max = julian([2004 11 1 0 0 0]);
jd = jday_min:1:jday_max;
gregaxm_mod(jd,3);
axis([jday_min jday_max 0 100]);
ylabel('germination potential');
legend('station 94')
legend boxoff

% Matrai, station 3
subplot(3,1,3)
aa = find(isfinite(st3));
% trick the legend
plot(0,0,'Marker','none','LineStyle','none');
hold on
plot(jday_mat(aa),st3(aa),'Color','k','Marker','o','MarkerSize',6)
% First 19 points are within the first year of germination
plot(jday_mat(aa(1:19)),st3(aa(1:19)),'Color','k','Marker','o', ...
     'MarkerSize',6,'MarkerFaceColor',[0 0 0]);
jday_min = julian([1999 11 1 0 0 0]);
jday_max = julian([2004 11 1 0 0 0]);
jd = jday_min:1:jday_max;
gregaxm_mod(jd,3);
axis([jday_min jday_max 0 100]);
legend('station 4');
legend boxoff
ylabel('germination potential','VerticalAlignment','top');

% Figure 9 of appendix A
figure(2)
% Old data, isolation 1
subplot(3,1,1)
% trick the legend
plot(0,0,'Marker','none','LineStyle','none');
hold on
aa = find(isfinite(iso1));
plot(jday_old(aa),iso1(aa),'Color','k','Marker','o','MarkerSize',6)
hold on
plot(jday_old(aa(2:7)),iso1(aa(2:7)),'Color','k','Marker','o', ...
     'MarkerSize',6,'MarkerFaceColor',[0 0 0]);
jday_min = julian([1983 11 1 0 0 0]);
jday_max = julian([1988 11 1 0 0 0]);
jd = jday_min:1:jday_max;
gregaxm_mod(jd,3);
axis([jday_min jday_max 0 100]);
ylabel('germination potential')
legend('station 29 (1984)');
legend boxoff

% Old data, isolation 2
subplot(3,1,2)
% trick the legend
plot(0,0,'Marker','none','LineStyle','none');
hold on
aa = find(isfinite(iso2));
plot(jday_old(aa),iso2(aa),'Color','k','Marker','o','MarkerSize',6)
hold on
plot(jday_old(aa(1:11)),iso2(aa(1:11)),'Color','k','Marker','o', ...
     'MarkerSize',6,'MarkerFaceColor',[0 0 0]);
jday_min = julian([1984 11 1 0 0 0]);
jday_max = julian([1989 11 1 0 0 0]);
jd = jday_min:1:jday_max;
gregaxm_mod(jd,3);
axis([jday_min jday_max 0 100]);
ylabel('germination potential')
legend('station 29 (1985)')
legend boxoff

% Old data, isolation 3: Note that this sample was isolated in June even
% though experiments weren't started until later.  Thus, it may seem like
% less then a year was taken based on the first data point, but data was
% actually taken up to a year after the cyst isolation.
subplot(3,1,3)
% trick the legend
plot(0,0,'Marker','none','LineStyle','none');
hold on
aa = find(isfinite(iso3));
plot(jday_old(aa),iso3(aa),'Color','k','Marker','o','MarkerSize',6)
plot(jday_old(aa(1:3)),iso3(aa(1:3)),'Color','k','Marker','o', ...
     'MarkerSize',6,'MarkerFaceColor',[0 0 0]);
jday_min = julian([1986 11 1 0 0 0]);
jday_max = julian([1991 11 1 0 0 0]);
jd = jday_min:1:jday_max;
gregaxm_mod(jd,3);
axis([jday_min jday_max 0 100]);
legend('station 29 (1986)');
legend boxoff
ylabel('germination potential');

figure(3)
% 1998 amy data
subplot(3,1,1)
% trick the legend
plot(0,0,'Marker','none','LineStyle','none');
hold on
aa = find(isfinite(amy));
plot(jday_amy(aa),amy(aa),'Color','k','Marker','o','MarkerSize',6)
plot(jday_amy(aa(2:27)),amy(aa(2:27)),'Color','k','Marker','o', ...
     'MarkerSize',6,'MarkerFaceColor',[0 0 0]);
jday_min = julian([1997 11 1 0 0 0]);
jday_max = julian([2002 11 1 0 0 0]);
jd = jday_min:1:jday_max;
gregaxm_mod(jd,3);
axis([jday_min jday_max 0 100])
legend('station 38');
legend boxoff
ylabel('germination potential')

% Isolate the data 1 year after isolation and before the second measured
% cycle starts.

% Experiment initiated for the first point on 12/27/1999, last point on
% 10/30/2000, clock function starts to begin a second cycle after the 
% 10/30/2000 point
aa = find(isfinite(st1));
mod61(:,1) = st1(aa(1:22));
mod61(:,2) = mon_mat(aa(1:22))';
mod61(:,3) = day_mat(aa(1:22))';

% Experiment initiated for the first point on 1/5/2000, last point on
% 11/21/2000, clock function starts to begin a second cycle after the 
% last point
aa = find(isfinite(st2));
mod94(:,1) = st2(aa(1:22));
mod94(:,2) = mon_mat(aa(1:22))';
mod94(:,3) = day_mat(aa(1:22))';

% Experiment initiated for the first point on 2/23/2000, last point on
% 11/14/2000, clock function starts to begin a second cycle after the 
% last point
aa = find(isfinite(st3));
mod4(:,1) = st3(aa(1:19));
mod4(:,2) = mon_mat(aa(1:19))';
mod4(:,3) = day_mat(aa(1:19))';

% The first point in this experiment was on 5/18/1984 and from an earlier
% isolation.  The first point is therefore considered to be on 11/23/1984.
% The last point is on 8/26/1985, clock function starts to begin a second 
% cycle after the last point.
aa = find(isfinite(iso1));
mod29_84(:,1) = iso1(aa(2:7));
mod29_84(:,2) = mon_old(aa(2:7))';
mod29_84(:,3) = day_old(aa(2:7))';

% The first point in this experiment was on 8/26/1985, the last point is
% from 7/21/1986 - this one has a somewhat odd pattern.
aa = find(isfinite(iso2));
mod29_85(:,1) = iso2(aa(1:11));
mod29_85(:,2) = mon_old(aa(1:11))';
mod29_85(:,3) = day_old(aa(1:11))';

% This isolate was isolated in June but experiments were not started until
% later.  The first point is on 1/22/1987, the last is on 6/3/1987
aa = find(isfinite(iso3));
mod29_87(:,1) = iso3(aa(1:3));
mod29_87(:,2) = mon_old(aa(1:3))';
mod29_87(:,3) = day_old(aa(1:3))';

% The first point in this series is from a previous isolation.  The first
% point considered is thus 2/12/1998, the last 1/28/1999
aa = find(isfinite(amy));
mod38(:,1) = amy(aa(2:27));
mod38(:,2) = mon_amy(aa(2:27))';
mod38(:,3) = day_amy(aa(2:27))';

% plot on uniform - do conversion to year day using function
% "yearday" which considers january 1st yearday 1.
% neglect the influence of the leap-year - we are trying to map 
% onto the same 365 day axis - we need to either treat them all
% as leap-years or non - leap-years...either way there is only
% a slight inaccuracy.

% Also, plot against mid-point of experiment, some experiments were 6 weeks
% long (add 21 days), others were 4 (add 14 days).  
bb = size(mod61,1);
for i = 1:bb
  temp = yearday([1995 mod61(i,2:3)]) + [0 21];
  yd61(i) = temp(1,2);
end
aa = find(yd61 > 365);
yd61(aa) = yd61(aa) - 365;

% sort in increasing order
clear temp;
temp = mod61;
clear mod61;
[yd61 sort_ind] = sort(yd61);
mod61 = temp(sort_ind);
clear temp;

% Repeat process for station 94
bb = size(mod94,1);
for i = 1:bb
  temp = yearday([1995 mod94(i,2:3)]) + [0 21];
  yd94(i) = temp(1,2);
end
aa = find(yd94 > 365);
yd94(aa) = yd94(aa) - 365;

clear temp;
temp = mod94;
clear mod94;
[yd94 sort_ind] = sort(yd94);
mod94 = temp(sort_ind);
clear temp;

% Repeat for station 4
bb = size(mod4,1);
for i = 1:bb
  temp = yearday([1995 mod4(i,2:3)]) + [0 21];
  yd4(i) = temp(1,2);
end
aa = find(yd4 > 365);
yd4(aa) = yd4(aa) - 365;

clear temp;
temp = mod4;
clear mod4;
[yd4 sort_ind] = sort(yd4);
mod4 = temp(sort_ind);
clear temp;

% 1984 isolate
bb = size(mod29_84,1);
for i = 1:bb
  temp = yearday([1995 mod29_84(i,2:3)]) + [0 14];
  yd29_84(i) = temp(1,2);
end
aa = find(yd29_84 > 365);
yd29_84(aa) = yd29_84(aa) - 365;

clear temp;
temp = mod29_84;
clear mod29_84;
[yd29_84 sort_ind] = sort(yd29_84);
mod29_84 = temp(sort_ind);
clear temp;

% 1985 isolate
bb = size(mod29_85,1);
for i = 1:bb
  temp = yearday([1995 mod29_85(i,2:3)]) + [0 14];
  yd29_85(i) = temp(1,2);
end
aa = find(yd29_85 > 365);
yd29_85(aa) = yd29_85(aa) - 365;

clear temp;
temp = mod29_85;
clear mod29_85;
[yd29_85 sort_ind] = sort(yd29_85);
mod29_85 = temp(sort_ind);
clear temp;

% 1986, 87 isolate
bb = size(mod29_87,1);
for i = 1:bb
  temp = yearday([1995 mod29_87(i,2:3)]) + [0 14];
  yd29_87(i) = temp(1,2);
end
aa = find(yd29_87 > 365);
yd29_87(aa) = yd29_87(aa) - 365;

clear temp;
temp = mod29_87;
clear mod29_87;
[yd29_87 sort_ind] = sort(yd29_87);
mod29_87 = temp(sort_ind);
clear temp;

bb = size(mod38,1);
for i = 1:bb
  temp = yearday([1995 mod38(i,2:3)]) + [0 14];
  yd38(i) = temp(1,2);
end
aa = find(yd38 > 365);
yd38(aa) = yd38(aa) - 365;

clear temp;
temp = mod38;
clear mod38;
[yd38 sort_ind] = sort(yd38);
mod38 = temp(sort_ind);
clear temp;

% plot all together with different color codings

figure(4)

hold on
plot(yd61,mod61,'k-','Marker','.','MarkerSize',20);
plot(yd94,mod94,'r-','Marker','.','MarkerSize',20);
plot(yd4,mod4,'b-','Marker','.','MarkerSize',20);
plot(yd29_84,mod29_84,'g-','Marker','.','MarkerSize',20);
plot(yd29_85,mod29_85,'y-','Marker','.','MarkerSize',20);
plot(yd29_87,mod29_87,'c-','Marker','.','MarkerSize',20);
plot(yd38,mod38,'m-','Marker','.','MarkerSize',20);
ylabel('germination potential')
xlabel('year-day, January 1 = day 1');
legend('st. 61','st. 94','st 4','st. 29, 84','st 29, 85','st 29, 87','st 38');
set(gca,'Xtick',[0:30:360]);

% lump all data together and calculate means and medians at monthly and
% bimonthly intervals.  Ended up using the monthly median (Table 4,
% Appendix A of Stock thesis

col1 = [yd61 yd94 yd4 yd29_84 yd29_85 yd29_87 yd38]';
col2 = [mod61 mod94 mod4 mod29_84 mod29_85 mod29_87 mod38]';
mod_data(:,1) = col1;
mod_data(:,2) = col2;

% monthly periods
mon_div = [0 31 59 90 120 151 181 212 243 273 304 334 365];

for i = 1:12
  i
  aa = find((mod_data(:,1) > mon_div(i)) & (mod_data(:,1) <= mon_div(i+1)));
  mod_data(aa,2)
  monthly_mean(i) = mean(mod_data(aa,2));
  monthly_median(i) = median(mod_data(aa,2));
  monthly_size(i) = size(aa,1);
  mon_xplot(i) = mon_div(i) + (mon_div(i+1) - mon_div(i))/2; 
end

% bi-monthly periods

bmon_div = [0 59 120 181 243 304 365];

for i = 1:6
  aa = find((mod_data(:,1) > bmon_div(i)) & (mod_data(:,1) <= bmon_div(i+1)));
  bmonthly_mean(i) = mean(mod_data(aa,2));
  bmonthly_median(i) = median(mod_data(aa,2));
  bmonthly_size(i) = size(aa,1);
  bmon_xplot(i) = bmon_div(i) + (bmon_div(i+1) - bmon_div(i))/2;
end

figure(5)
subplot(2,1,1);
plot(mon_xplot,monthly_median,'k-','Marker','o','MarkerSize',8,'MarkerFaceColor',[0 0 0])
hold on
plot(mon_xplot,monthly_mean,'--','Marker','o','MarkerSize',8,'MarkerFaceColor', ...
     [0.7 0.7 0.7],'Color',[0.7 0.7 0.7]);
legend('median','mean');
title('Monthly Bins');
ylabel('germination potential')
xlabel('year-day, January 1 = day 1');
set(gca,'Xtick',[0:30:360]);
axis([0 370 0 100]);

subplot(2,1,2)
plot(bmon_xplot,bmonthly_median,'k-','Marker','o','MarkerSize',8,'MarkerFaceColor',[0 0 0])
hold on
a = plot(bmon_xplot,bmonthly_mean,'--','Marker','o','MarkerSize',8,'MarkerFaceColor', ...
     [0.7 0.7 0.7],'MarkerEdgeColor',[0.7 0.7 0.7],'Color',[0.7 0.7 0.7]);
legend('median','mean');
title('Bi-monthly Bins')
ylabel('germination potential')
xlabel('year-day, January 1 = day 1');
set(gca,'Xtick',[0:30:360]);
axis([0 370 0 100]);
