%========================
%Data Input
%========================
clear all;
%Enter your input filename here (including file extension) - 
filename = 'MG-GFP-1series-1.txt';
%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, LabelNames] = platereaderread3(filename);

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

%========================
%set data parameters
%========================
timepoints = [1:98];
Absbackgroundsamples = [];
GFPbackgroundsamples = [];
experimentsamples = [1:3];
names = ['1';...
         '2';...
         '3'];
wellnames = cellstr(names);
replicates = 3;
samples = [experimentsamples,GFPbackgroundsamples,Absbackgroundsamples];
Absconversion = 'OD';
%========================
%set display parameters
%========================
plotunsmooth = 0;
plotsmooth = 1; 
%Define some figure handles
raw = 1;
processed = 2;
if Absconversion = 'Ab'
    cellunits = 'Absorbance @ 600nm';
elseif Absconversion = 'OD'
    cellunits = 'Optical Density (600nm)';
elseif Absconversion = 'ce'
    cellunits = 'Cells/ml';
end

%========================
%Reshape data
%========================

Time = TimeData(timepoints);
rawAbs = LabelData(timepoints,:,1);
rawGFP = LabelData(timepoints,:,2);
%reshape time series 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]);
uAbs = orderedAbs(:,samples,:);
uGFP = orderedGFP(:,samples,:);
   

%========================
%Time Series 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} = uAbs;
    timeseries{2} = uGFP;
    [smoothtimeseries] = platereadersmooth(timeseries,method,cutofffreq,neighbors);
    %now get the smoothed data back out into individual arrays.
    fAbs = smoothtimeseries{1};
    fGFP = smoothtimeseries{2};
else
    fAbs = uAbs;
    fGFP = uGFP;
end


%========================
%Background subtraction
%========================
%This implementation assumes the OD of your samples is similar to that of the
%control at each time point.

if ~isempty(Absbackgroundsamples)
    correcteduAbs = Absbackgroundsubtract(uAbs);
    correctedfAbs = Absbackgroundsubtract(fAbs);
else
    correcteduAbs = uAbs;
    correctedfAbs = fAbs;
end

if ~isempty(GFPbackgroundsamples)
    [correctedGFP]=FPbackgroundsubtract(fGFP);
    [ucorrectedGFP]=FPbackgroundsubtract(uGFP)
else
    correctedGFP = fGFP;
    ucorrectedGFP = uGFP;    
end

%========================
%Unit conversion
%========================
if Absconversion = 'Ab'
    ucells = uAbs;
    fcells = fAbs;
elseif Absconversion = 'OD'
    ucells = 3.113.*uAbs-0.0158;
    fcells = 3.113.*fAbs-0.0158;
elseif Absconversion = 'ce'

end
    
%===========================

%Calculate the dilution rate, assuming exponential growth
[Mfgammad, CIfgammad] = AbsProcess(Time,fcells);
[Mugammad, CIugammad] = AbsProcess(Time,ucells);

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

[MfGFPperOD, CIfGFPperOD, MdfGFPperODdt, CIdfGFPperODdt, Mfdegradation, CIfdegradation, Mfsynthesis, CIfsynthesis] = FPProcess(Time,fcells,fgammad,correctedGFP);
[MuGFPperOD, CIuGFPperOD, MduGFPperODdt, CIduGFPperODdt, Mudegradation, CIudegradation, Musynthesis, CIusynthesis] = FPProcess(Time,ucells,ugammad,ucorrectedGFP);

%We should average each of the raw/smooth time series replicates
[Mucells,CIucells] = platereaderstats(ucells);
[Mfcells,CIfcells] = platereaderstats(fcells);

[MuGFP,CIuGFP] = platereaderstats(uGFP);
[MfGFP,CIfGFP] = platesdreaderstats(fGFP);

%========================
%Data plotting
%========================

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,Mucells(:,1:figure1columns),'s','MarkerSize',5)
        hold on;
     
    end
    
    if plotsmooth
      semilogy(Time,Mfcells(:,1:figure1columns),'LineWidth',1.5); 
    end
    
    title('Raw and smoothed OD')
    xlabel('Time (mins)')
    ylabel(cellunits);
legend(num2str(experimentsamples'),'Location','best');
subplot(1,2,2)
    hold off;
    if plotunsmooth
        errorbar(Mfcells(:,1:figure1columns),MuGFP(:,1:figure1columns),CIuGFP(:,1:figure1columns),'o','MarkerSize',5)
        hold on;
    end
    
    if plotsmooth
        errorbar(Mfcells(:,1:figure1columns),MfGFP(:,1:figure1columns),CIfGFP(:,1:figure1columns),'LineWidth',0.5)
    end
    
    title('Raw and smoothed GFP')
    xlabel(cellunits);
    ylabel('GFP counts');
%     xlim([0.1 0.4]);
figure(processed(1));
subplot(3,2,1)
    hold off;
    if plotunsmooth
        errorbar(Mfcells(:,1:figure2columns),MuGFPperOD(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
    end
    
    if plotsmooth
        errorbar(Mfcells(:,1:figure2columns),MfGFPperOD(:,1:figure2columns),'LineWidth',0.5)
    end
    title('GFP per OD');
ylim([-2000, 150000]);
    subplot(3,2,2)
    hold off; 
    if plotunsmooth
        errorbar(Mfcells(2:end,1:figure2columns),MduGFPperODdt(:,1:figure2columns),CIduGFPperODdt(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
        ylim([-200, 500]);
    end
    if plotsmooth
            errorbar(Mfcells(2:end,1:figure2columns),MdfGFPperODdt(:,1:figure2columns),CIdfGFPperODdt(:,1:figure2columns),'LineWidth',0.5)
    end
    title('Accumulation rate of GFP per OD');
    
subplot(3,2,3)
    hold off;
    if plotunsmooth
        errorbar(Mfcells(2:end,1:figure2columns),Mugammad(:,1:figure2columns).*60,CIugammad(:,1:figure2columns).*60,'o','MarkerSize',5)
        hold on;
    end
    if plotsmooth
        errorbar(Mfcells(2:end,1:figure2columns),Mfgammad(:,1:figure2columns).*60,CIfgammad(:,1:figure2columns).*60,'LineWidth',0.5)
    end
    title('Dilution rate (doublings/hour)');
    ylim([-1,2]);
    xlabel(cellunits);

subplot(3,2,4)
    hold off;
    if plotunsmooth
        errorbar(Mfcells(2:end,1:figure2columns),Mudegradation(:,1:figure2columns),CIudegradation(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
    end
    if plotsmooth
        errorbar(Mfcells(2:end,1:figure2columns),Mfdegradation(:,1:figure2columns),CIfdegradation(:,1:figure2columns),'LineWidth',0.5)
    end
   title('Degradation rate');
    ylim([-200,500])
    xlabel(cellunits);

subplot(3,2,5)
    hold off;
    if plotunsmooth
        errorbar(Mfcells(2:end,1:figure2columns),Musynthesis(:,1:figure2columns),CIusynthesis(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
    end
    if plotsmooth
        errorbar(Mfcells(2:end,1:figure2columns),Mfsynthesis(:,1:figure2columns),CIfsynthesis(:,1:figure2columns),'LineWidth',0.5)
    end
    title('Synthesis rate');
    ylim([-200, 700]);
    xlabel(cellunits);
    legend(wellnames,'Location','Best');

% 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));