clear all;
%To use this code, just enter your input file prefix here - 
filename = 'MG-GFP-1series-4';
%and the path to where your data should be saved.  If you don't specify an
%output location, it will  use the current directory.
% outputpath = '~/Desktop/Chassis-System/Data/PlateReader/DedicatedTranslation/';


[TimeData, LabelData] = platereaderread2MG(filename);

%%
%============
%set smoothing parameters
%============
method = 'fb';
if method == 'fb' | method == 'ff'
    cutofffreq = 0.15;
    neighbors = [];
elseif method == 'co'
    cutofffreq = 0.15;
    neighbors = 9;
elseif method == 'ma'
    cutofffreq = [];
    neighbors = 9;
end
plotunsmooth = 1;
plotsmooth = 1; 

Absbackgroundsamples = [];
GFPbackgroundsamples = [];

experimentsamples = [30:32];

names = ['RBS 2.8-1A3 Induced       ';...
         'RBS 2.8-1A3 Uninduced     ';...
         'RBS 0.3-1A3 Induced       ';...
         'RBS 0.3-1A3 Uninduced     ';...
         'RBS 8-1A3 Induced         ';...
         'RBS 8-1A3 Uninduced       ';...
         'RBS 4B-1A3 Induced        ';...
         'RBS 4B-1A3 Uninduced      ';...
         'RBS 2-1A3 Induced         ';...
         'RBS 2-1A3 Uninduced       ';...
         'Non-RBS HK Induced        ';...
         'Non-RBS HK Uninduced      ';...
         'Non-RBS MG Induced        ';...
         'Non-RBS MG Uninduced      ';...
         'RBS 2.8-4A3 Induced       ';...
         'RBS 2.8-4A3 Uninduced     ';...
         'RBS 0.3-4A3 Induced       ';...
         'RBS 0.3-4A3 Uninduced     ';...
         'RBS 8-4A3 Induced         ';...
         'RBS 8-4A3 Uninduced       ';...
         'RBS 4B-4A3 Induced        ';...
         'RBS 4B-4A3 Uninduced      ';...
         'RBS 2-4A3 Induced         ';...
         'RBS 2-4A3 Uninduced       ';...
         'Vector 1A3 Empty Induced  ';...
         'Vector 1A3 Empty Uninduced';...
         'Vector 4A3 Empty Induced  ';...
         'Vector 4A3 Empty Uninduced';...
         'MG1655                    ';...
         'Blank - EZ                ';...
         'Blank - EZ + Amp          ';...
         'Blank - EZ + Amp + IPTG   '];

wellnames = cellstr(names(experimentsamples,:));

% All RBS's 1A3: 1:14,25,26,29
% All RBS's 1A3 induced: 1,3,5,7,9,11,13,25,29
% All RBS's 1A3 uninduced: 2,4,6,8,10,12,14,26,29
% All RBS's 4A3: 15:24,27,28,29
% All RBS's 4A3 induced: 15,17,19,21,23,27,29
% All RBS's 4A3 uininduced: 16,18,20,22,24,28,29
% Empty Vectors: 25:28
% Media: 30:32

% Subtraction Scheme
% Note: The EZ and EZ + Amp wells (30,31) were contaminated, so I used EZ +
% Amp + IPTG (32). 
% Absbackgroundsamples = [32];
% GFPbackgroundsamples = [29];

% Time Interval Used
% timepoints = [15:85];

replicates = 3;
samples = [experimentsamples,GFPbackgroundsamples,Absbackgroundsamples];
% Extract the data that we are interested in into three new variables.  By
% default, consider all timepoints and wells
% timepoints = [1:size(LabelData,1)-1];
wells = [1:size(LabelData,2)];
timepoints = [15:85];
% wells = [1:6];

rawGFP = LabelData(timepoints,wells,2);
rawAbs = LabelData(timepoints,wells,1);
Time = TimeData(timepoints);

%If your data has a fixed number of replicates for each sample and control,
%and the replicates are in adjacent wells, use this code to reshape both 
%matrices so that the third dimension represents the different replicates
orderedAbs = reshape(rawAbs,[size(rawAbs,1),replicates,size(rawAbs,2)/replicates]);
orderedGFP = reshape(rawGFP,[size(rawGFP,1),replicates,size(rawGFP,2)/replicates]);
%reorder the dimensions of the array so that the third dimension is the
%number of replicates.
orderedAbs = permute(orderedAbs,[1,3,2]);
orderedGFP = permute(orderedGFP,[1,3,2]);
orderedAbs = orderedAbs(:,samples,:);
orderedGFP = orderedGFP(:,samples,:);

%Background subtraction
%If you are subtracting one time course from another as a background
%subtraction, it would be better to wait until after the data has been
%filtered.  For now, we'll just subtract the mean of a media only well from
%the absorbance.
   
if ~isempty(Absbackgroundsamples)
    Absbackground = mean(orderedAbs(:,end,:),3);
    cellAbs = orderedAbs - repmat(Absbackground,[1,size(orderedAbs,2),size(orderedAbs,3)]);
   
else 
    cellAbs = orderedAbs;
    
end

cellOD = 3.113.*cellAbs-0.0158;

% cellOD = cellAbs;
%==============================
%Data Smoothing
%Put the data to be smoothed into a cell array so the function can handle
%arbitrary numbers of labels.
if ~isempty(method)
    timeseries{1} = cellOD;
    timeseries{2} = orderedGFP;
    [smoothtimeseries] = platereadersmooth(timeseries,method,cutofffreq,neighbors);
    %now get the smoothed data back out into individual arrays.
    fOD = smoothtimeseries{1};
    OD = cellOD;
    fGFP = smoothtimeseries{2};
    GFP = orderedGFP;
else
    fOD = cellOD;
    fGFP = orderedGFP;
end


%===========================
%Background Subtraction
%You should probably be subtracting the fluorescence of cells not
%expressing GFP so the GFP signal is truly proportional to GFP.  This
%implementation assumes the OD of your samples is similar to that of the
%control at each time point.
% backgroundsample = size(fGFP,2);
if ~isempty(GFPbackgroundsamples)
    meanbackground = mean(fGFP(:,end-1,:),3);
    %Need to make this array the same size as the full array so we can subtract
    %arrays rather than using slow looping.
    meanbackgroundarray = repmat(meanbackground,[1,size(fGFP,2),size(fGFP,3)]);
    correctedGFP = fGFP - meanbackgroundarray;
    
    
    umeanbackground = mean(GFP(:,end-1,:),3);
    %Need to make this array the same size as the full array so we can subtract
    %arrays rather than using slow looping.
    umeanbackgroundarray = repmat(umeanbackground,[1,size(GFP,2),size(GFP,3)]);
    ucorrectedGFP = GFP - umeanbackgroundarray;
    
else
    correctedGFP = fGFP;
    
    
    ucorrectedGFP = GFP;
    
end


%===========================
%GFP Synthesis Rate Calculations

%find the amount of GFP per OD
GFPperOD = correctedGFP./fOD;
uGFPperOD = ucorrectedGFP./OD;
%Calculate the rate of accumulation of GFP per OD per unit time
dGFPperOD = diff(GFPperOD);
udGFPperOD = diff(uGFPperOD);
dt = mean(diff(Time));
dGFPperODdt = dGFPperOD./dt;
udGFPperODdt = udGFPperOD./dt;
%Calculate the dilution rate, assuming exponential growth
dOD = diff(fOD);
udOD = diff(OD);

dODdt = dOD./dt;
udODdt = udOD./dt;

gammad = (1/log(2))*(dODdt./fOD(2:end,:,:));
ugammad = (1/log(2))*(udODdt./OD(2:end,:,:));
%calculate the degradation rate of GFP per OD
degradation = gammad.*GFPperOD(2:end,:,:);
udegradation = ugammad.*uGFPperOD(2:end,:,:);
%The synthesis rate is the sum of the accumulation rate and the degradation
%rate.
synthesis = dGFPperODdt + degradation;
usynthesis = udGFPperODdt + udegradation;

%At this point, we should average each of the replicates
meansynthesis = mean(synthesis,3);
umeansynthesis = mean(usynthesis,3);

unsmoothGFP = mean(orderedGFP,3);
unsmoothOD = mean(cellOD,3);
smoothGFP = mean(fGFP,3);
smoothOD = mean(fOD,3);

meanGFPperOD = mean(GFPperOD,3);
umeanGFPperOD = mean(uGFPperOD,3);

meandGFPperODdt = mean(dGFPperODdt,3);
umeandGFPperODdt = mean(udGFPperODdt,3);

meangammad = mean(gammad,3);
umeangammad = mean(ugammad,3);

meandegradation = mean(degradation,3);
umeandegradation = mean(udegradation,3);

%=======================
%Data Plotting

%Define some figure handles
raw = 1;
processed = 2;
surface = 4;
ODwindow = smoothOD > 0.1 & smoothOD < 0.3;
figure1columns = size(experimentsamples+1,2);
figure2columns = size(experimentsamples,2);
% timewindow = 10:97;
set(0,'DefaultAxesColorOrder',jet(figure1columns))
figure(raw);
subplot(1,2,1)
    hold off;
    if plotunsmooth
        semilogy(Time,unsmoothOD(:,1:figure1columns),'o','MarkerSize',5)
        hold on;
     
    end
    
    if plotsmooth
      semilogy(Time,smoothOD(:,1:figure1columns),'LineWidth',1.5); 
    end
    
    title('Raw and smoothed OD')
    xlabel('Time (mins)')
    ylabel('OD600');
legend(wellnames,0);
   % ylim([0.1 1]);
subplot(1,2,2)
    hold off;
    if plotunsmooth
        plot(smoothOD(:,1:figure1columns),unsmoothGFP(:,1:figure1columns),'o','MarkerSize',5)
        hold on;
    end
    
    if plotsmooth
        plot(smoothOD(:,1:figure1columns),smoothGFP(:,1:figure1columns),'LineWidth',1.5)
    end
    title('Raw and smoothed GFP')
    xlabel('OD600')
    ylabel('GFP counts');
%     xlim([0.1 0.4]);
figure(processed(1));
subplot(3,2,1)
    hold off;
    if plotunsmooth
        plot(smoothOD(:,1:figure2columns),umeanGFPperOD(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
    end
    
    if plotsmooth
        plot(smoothOD(:,1:figure2columns),meanGFPperOD(:,1:figure2columns),'LineWidth',1.5)
    end
    title('GFP per OD');
% ylim([-2000, 150000]);
    subplot(3,2,2)
    hold off; 
    if plotunsmooth
        plot(smoothOD(2:end,1:figure2columns),umeandGFPperODdt(:,1:figure2columns),'o','MarkerSize',5)
        hold on;

    end
    if plotsmooth
            plot(smoothOD(2:end,1:figure2columns),meandGFPperODdt(:,1:figure2columns),'LineWidth',1.5)
    end
    title('Accumulation rate of GFP per OD');
    %         ylim([-200, 500]);
subplot(3,2,3)
    hold off;
    if plotunsmooth
        plot(smoothOD(2:end,1:figure2columns),umeangammad(:,1:figure2columns).*60,'o','MarkerSize',5)
        hold on;
    end
    if plotsmooth
        plot(smoothOD(2:end,1:figure2columns),meangammad(:,1:figure2columns).*60,'LineWidth',1.5)
    end
    title('Dilution rate (doublings/hour)');
    ylim([-1,2]);
subplot(3,2,4)
    hold off;
    if plotunsmooth
        plot(smoothOD(2:end,1:figure2columns),umeandegradation(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
    end
    if plotsmooth
        plot(smoothOD(2:end,1:figure2columns),meandegradation(:,1:figure2columns),'LineWidth',1.5)
    end
   title('Degradation rate');
% ylim([-200,500])
    subplot(3,2,5)
    hold off;
    if plotunsmooth
        plot(smoothOD(2:end,1:figure2columns),umeansynthesis(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
    end
    if plotsmooth
        plot(smoothOD(2:end,1:figure2columns),meansynthesis(:,1:figure2columns),'LineWidth',1.5)
    end
    title('Synthesis rate');
%     ylim([-200, 700]);
legend(wellnames,0);

% figure(surface);
% timeaxis = Time(2:end);
% samplestoplot = [1:4:29];
% [xx,yy] = meshgrid(timeaxis,1:length(samplestoplot));
% mesh(xx,yy,meansynthesis(:,samplestoplot)');
% xlabel('Time (mins)')
% ylabel('Plotted Sample #');
%%

indices = fOD>0.1 & fOD<0.3;

numberofsamplesinrange = sum(indices);
correspondingsynthesisrates = synthesis.*indices(2:end,:,:);
meansynthesisrates = sum(correspondingsynthesisrates)./numberofsamplesinrange;

synthesisratesminusmean = synthesis-repmat(meansynthesisrates,[size(synthesis,1),1,1]);
synthesisratesminusmean = synthesisratesminusmean.*indices(2:end,:,:);
synthesisratesminusmeansquared = synthesisratesminusmean.^2;
varsynthesisrates = sum(synthesisratesminusmeansquared)./numberofsamplesinrange;
sdsynthesisrates = varsynthesisrates.^0.5;
sesynthesisrates = sdsynthesisrates./sqrt(numberofsamplesinrange);
samplemeansynthesisrates = squeeze(mean(meansynthesisrates,3));
samplestsynthesisrates = squeeze(mean(sdsynthesisrates,3));
samplesesynthesisrates = squeeze(mean(sesynthesisrates,3));

% figure


% [b,a] = butter(2,0.125);
% label = 1;
% d = diff(LabelData(timepoints,wells,label));
% fd = filter(b,a,d);
% f = filter(b,a,LabelData(timepoints,wells,label) - mean(LabelData(1,wells,label)));
% f = f + mean(LabelData(1,wells,label));
% s = smooth(LabelData(timepoints,wells(1),label),13,'moving');
% noise = LabelData(timepoints,wells(1),label) - s;
% 
% 
% df = diff(f);
% % wells = [48];
% fftplot = 3;
% figure(fftplot)
% hold off;
% subplot(1,3,1)
% hold off;
% PlotFFT(noise(:,1),1/mean(diff(TimeData)));
% % PlotFFT(diff(LabelData(timepoints,wells(1),label)),1/mean(diff(TimeData)));
% hold on;
% % PlotFFTred(LabelData(timepoints,wells(1),label)-mean(LabelData(timepoints
% % ,wells(2),label)),1/mean(diff(TimeData)));
% 
% subplot(1,3,2)
% hold off;
% plot(TimeData(timepoints(1:end-1)),d,'.-b');
% hold on;
% plot(TimeData(timepoints(1:end-1)),df,'.-r');
% plot(TimeData(timepoints(1:end-1)),fd,'.-g');
% subplot(1,3,3)
% plot(TimeData(timepoints),LabelData(timepoints,wells,label),'.-b');
% hold on;
% plot(TimeData(timepoints),f,'.-r');
%%
% plot(fdifferentials2,'-g');
% fdifferentials(50:end) = 0;
% differentials2 = ifft(fdifferentials,512);
% differentials2 = abs(differentials2(1:143));


% sabsorbance = smooth(LabelData(:,wells,1),17,'rloess');
% sfluorescence = smooth(LabelData(:,wells,2),17,'rloess');
% snormalized = sfluorescence./sabsorbance;
% sdifferentials = diff(snormalized)./mean(diff(TimeData));


% s2absfft = fft(LabelData(:,wells,1),512);
% s2absfft(255:end) = 0;
% s2absorbance = ifft(s2absfft,512);
% s2absorbance = s2absorbance(1:length(TimeData));
% s2flufft = fft(LabelData(:,wells,2),512);
% s2flufft(255:end) = 0;
% s2fluorescence = ifft(s2flufft,512);
% s2fluorescence = s2fluorescence(1:length(TimeData));
% 
% s2normalized = s2fluorescence./s2absorbance;
% s2differentials = diff(s2normalized)./mean(diff(TimeData));

% platereaderplotting2(TimeData,LabelData(:,wells,1),LabelData(:,wells,2),normalized,differentials)
% % figure
% % platereaderplotting2(TimeData,sabsorbance,sfluorescence,snormalized,sdifferentials)
% figure
% platereaderplotting2(TimeData,LabelData(:,wells,1),LabelData(:,wells,2),normalized,fdifferentials)

% platereaderplotting2(TimeData,s2absorbance(1:length(TimeData)),s2fluorescence,s2normalized,s2differentials)

%legend(LabelNames(wells,:),'Location','EastOutside');


%[dadt] = calculategrowthrate(time, absorbance)
%%
% figure
% transformed = fft(LabelData(:,wells,1),512);
% P = transformed.*conj(transformed)/512;
% f = 0.2*(0:256)/512;
% plot(f,P(1:257))
% %%
% diffs = fft(differentials,512);
% P3 = diffs.*conj(diffs)/512;
% %figure
% plot(f,P3(1:257))
% [b,a] = butter(2,0.15)
% filtered = filter(b,a,differentials);
% 
% range = 132:138;
% diffs(range) = 0;
% diffs(512-range) = 0;
% rev = ifft(diffs,512);
% figure
% hold on;
% plot(TimeData,rev(1:length(TimeData)),'-g')
% %plot(TimeData(1:end-1),differentials,'-r')
% plot(TimeData(1:end-1),filtered,'-b')

% windowSize = 5;
% manAbs = filter(ones(1,windowSize)/windowSize,1,nAbs(:,wells)-repmat(nAbs(1,wells),length(timepoints),1));
% manAbs = manAbs+repmat(nAbs(1,wells),length(timepoints),1);
% 
% manGFP = filter(ones(1,windowSize)/windowSize,1,nGFP(:,wells)-repmat(nGFP(1,wells),length(timepoints),1));
% manGFP = manGFP+repmat(nGFP(1,wells),length(timepoints),1);
% 
% sAbs(:,1) = smooth(nAbs(:,wells(1)),5,'moving');
% sAbs(:,2) = smooth(nAbs(:,wells(2)),5,'moving');
% 
% sGFP(:,1) = smooth(nGFP(:,wells(1)),5,'moving');
% sGFP(:,2) = smooth(nGFP(:,wells(2)),5,'moving');
% 
% figure(5);
% 
% subplot(1,2,1)
% hold off;
% plot(Time,nGFP(:,wells),'-b');
% hold on;
% plot(Time,fGFP,'-r');
% plot(Time,manGFP,'-g');
% plot(Time,sGFP,'-k');
% 
% subplot(1,2,2)
% hold off;
% plot(Time,nAbs(:,wells),'-b');
% hold on;
% plot(Time,fAbs,'-r');
% plot(Time,manAbs,'-g');
% plot(Time,sAbs,'-k');

% nGFP = LabelData(timepoints,:,2) - repmat(LabelData(1,:,2),length(timepoints),1);


% [b,a] = butter(2,0.15);
% fGFP = filter(b,a,nGFP(:,wells)-repmat(nGFP(1,wells),length(timepoints),1));
% fGFP = fGFP+repmat(nGFP(1,wells),length(timepoints),1);
% 
% [b2,a2] = butter(2,0.3);
% fAbs = filter(b2,a2,nAbs(:,wells)-repmat(nAbs(1,wells),length(timepoints),1));
% fAbs = fAbs+repmat(nAbs(1,wells),length(timepoints),1);

% noise = diff(nAbs(:,wells));
% figure(10)
% plot(Time(1:end-1),noise);
% figure(11)
% PlotFFT(noise(:,1),1/mean(diff(TimeData)));
% hold on;
% PlotFFTred(noise(:,2),1/mean(diff(TimeData)));
% PlotFFT(noise(:,3),1/mean(diff(TimeData)));
% PlotFFTred(noise(:,4),1/mean(diff(TimeData)));
% PlotFFT(noise(:,5),1/mean(diff(TimeData)));
% PlotFFTred(noise(:,6),1/mean(diff(TimeData)));

%asmoothOD = mean(smoothOD,1)
%bsmoothOD = mean(asmoothOD,2)
%deviation = std(asmoothOD)

%differences = [asmoothOD(:,1);asmoothOD(:,2);asmoothOD(:,3);asmoothOD(:,4);asmoothOD(:,5);asmoothOD(:,6)]-asmoothOD(:,7)
%deviations = differences/.0037

% smoothOD = mean(fOD,3)
% newOD = smoothOD(15,:)
% 
% differences = [newOD(1)-newOD(2),newOD(1)-newOD(3),newOD(1)-newOD(4),newOD(1)-newOD(5);
%               newOD(2)-newOD(3),newOD(2)-newOD(4),newOD(2)-newOD(5),0;
%               newOD(3)-newOD(4),newOD(3)-newOD(5),0,0;
%               newOD(4)-newOD(5),0,0,0]
%           
% positivedifferences = abs(differences)
% 
% standarddivs = positivedifferences/.034131