%% multiplot_dwelltimes.m
% Written for Berro Lab at Yale University
% Last edited by Mike Lacy, January 2019

% Script to take [sets] and get all track lifetimes, excluding tracks that
% started before time cutoff. Plots PDFs, CDFs, and more.

% If you want to use a different startcutoff for each set, first
% process the sets() struct separately to remove individual tracks below  
% each chosen cutoff and then call plot_dwelltimes(sets, 0).

function [lts, setmeans, setns] = multiplot_dwelltimes(sets, startcutoff)
%% Inputs:
% sets = data struct array that is the output from process_tracks. You can 
% pass the whole struct or one or a few members, e.g. mysets(1)

minlt = 0.3;      % time, in sec, of minimum track length
doPlots = 1;
finalsetgray = 0; % for plots. set to 1 if the last member of sets is a control

lts = {};
for i=1:numel(sets)
lts{i} = [];
for j=1:numel(sets(i).datasets)
    if(~isempty(sets(i).datasets(j).tracks))
        datasetstart = min([sets(i).datasets(j).tracks.start]);
        allstart = [sets(i).datasets(j).tracks.start] - datasetstart;
        l = [sets(i).datasets(j).tracks.lifetime];
        lts{i} = [lts{i} l(allstart>startcutoff & l>=minlt)];
    end
end
end

if(doPlots)
%make CDF plot containing all the histograms in one:
figure, hold on;
leg = {};
for i=1:numel(sets)
    col = [i/numel(sets) 0 1-(i/numel(sets))];
    if(i==numel(sets) && finalsetgray)
        col = [0.2 0.2 0.2];
    end
    histogram(lts{i},(minlt-0.001):0.1:25,'Normalization','cdf','DisplayStyle','stairs','LineWidth',2,'EdgeColor',col);
        % (minlt-0.001) is to avoid plotting artifacts of time binning,
        % since 1/10 ~= 0.1000 due to floating point precision
    leg{i} = sprintf('%s, n=%d, mean = %.2g +/- %.2g',sets(i).dirname, numel(lts{i}),mean(lts{i}),std(lts{i}));
end
legend(leg,'Location','southeast');
ylim([0 1]);
xlabel('lifetime, sec'); ylabel('cumulative fraction of events');
title(['histograms of all sets, tracks with start > ', num2str(startcutoff)]);
a=gca;
a.FontSize = 12;

% make a series of strips for each pdf plot
figure
for i=1:numel(lts)
subplot(numel(lts),1,i);
hold on
col = [i/numel(sets) 0 1-(i/numel(sets))];
if(i==numel(sets) && finalsetgray)
    col = [0.2 0.2 0.2];
end
histogram(lts{i},[(minlt-0.001):0.1:6, 30],'Normalization','probability','DisplayStyle','stairs','LineWidth',1.5,'EdgeColor',col);
a=gca;
a.FontSize = 10;
title(sprintf('%sp, n=%d, mean = %.2g s, 95%%ile = %.2g s',sets(i).dirname, numel(lts{i}),mean(lts{i}),prctile(lts{i},95)));
ylim([0 0.1]); xlim([0 6.1]);
if(i<numel(lts))
    a.XTickLabels = {};
end
end
xlabel('lifetime, sec')

[setmeans, setns] = datasets_error(sets);

end %end plotting
end %end function

%% Reproducibility error %%
% Compare lifetimes data across datsets within each set
%
% Run after sets() has been generated from plotmetrics

function [dataset_means, dataset_ns] = datasets_error(sets)
minlifetime = 0.3;
minN = 40;  %minimum n spots in a dataset to include it in the overall mean and SD.

datasetSizes = [];
m = [];
stdevs = [];
fullmeans = [];
dataset_means = {};
dataset_ns = {};
fullmedians = [];
pct95s = [];
sprintf('set\t| name\t| median, s\t| mean, s\t| stdev, s\t| 95%%ile \t| total n tracks \t| n datasets of >40 tracks')

for i=1:numel(sets)
    l = {};
    %get all the lifetimes in this dataset that meet criteria for start
    %cutoff and min lifetime. basically the same as lts{} array
    for j=1:numel(sets(i).datasets)
          if(~isempty(sets(i).datasets(j).tracks))
            datasetstart = min([sets(i).datasets(j).tracks.start]);
            allstart = [sets(i).datasets(j).tracks.start] - datasetstart;
            temp = [sets(i).datasets(j).tracks.lifetime];
            l{j} = temp(allstart > startcutoff & temp>minlifetime);
          end
    end

    meanlts = [];
    for j=1:numel(l)
        datasetSizes = [datasetSizes; numel(l{j})];
        if numel(l{j})>minN
            meanlts(j) = mean(l{j});
        else
            meanlts(j) = NaN;
        end
    end
    nans = isnan(meanlts);
    meanlts(nans) = [];
    nbig = numel(meanlts);
    m(i) = mean(meanlts);
    stdevs(i) = std(meanlts);
    fullmeans(i) = mean([l{:}]);
    dataset_means{i} = [];
    dataset_ns{i} = [];
    for j=1:numel(l)
        dataset_means{i}(j) = mean(l{j});
        dataset_ns{i}(j) = numel(l{j});
    end
    fullmedians(i) = median([l{:}]);
    pct95s(i) = prctile([l{:}],95);
    % print data table to the console:
    sprintf('%d \t%s \t%.3g \t%.3g \t%.3g \t%.3g \t%d \t%d',...
        i, sets(i).dirname, fullmedians(i),fullmeans(i),stdevs(i),pct95s(i),numel([l{:}]),nbig)
end

figure;
hold on;
bar(1:numel(sets), fullmeans);
errorbar(1:numel(sets),fullmeans,stdevs, '.k');
for i=1:numel(dataset_means)
    jitterx = rand(size(dataset_means{i}(dataset_ns{i}>minN)))*0.4 - 0.2; %random number from -0.2 to 0.2
    plot(jitterx + i, dataset_means{i}(dataset_ns{i}>minN), 'ok', 'MarkerSize', 4, 'MarkerFaceColor', [0.75 0.75 0.75]);
    
end

ylabel('mean lifetimes, sec')
title(sprintf('sm tracking lifetimes and stdev across datasets, cutoff starttime > %g s', startcutoff));
legend({'full mean of all tracks in set',...
    ['stdev across datasets in set' newline ' only taking datasets of n > ', num2str(minN)],...
    'individual dataset means'},'Location', 'southoutside','FontSize',12);
a=gca;
a.XTick = 1:numel(sets);
a.XTickLabels = {sets.dirname};
a.XTickLabelRotation = 30;
a.FontSize = 11;
end