% DobosiewiczCode -- Dobosiewicz M, Liu Q, and Bargmann CI, eLife
% Input a 301-frame calcium trace "FullMatrix1"


function Output = DobosiewiczCode(FullMatrix1)
Pulses1 = Cleanup(FullMatrix1); %removes traces that have NaN
tracenum1 = size(Pulses1,1);

% Normalize and smooth trace
Baseline1=NaN([tracenum1 1]); Norm1=NaN(size(Pulses1)); ToUse1=NaN(size(Pulses1));Tderiv1 = NaN(size(Pulses1));
for m=1:tracenum1
    F0=median(Pulses1(m,1:99)); 
    Baseline1(m,1)=F0;
    Norm1(m,:)=((Pulses1(m,:)-F0)/F0*100);
    ToUse1(m,:)=smooth(Norm1(m,:),5);
    Tderiv1(m,:) = tderiv(ToUse1(m,:));
end

% EVENT DETECTION
eventtimes1=NaN([size(ToUse1,1) 1]);
start=98;% the frame of stimulus onset, adjusted because of STD calculation, but gets adjusted back
for trace=1:tracenum1
    STD1=std(ToUse1(trace,1:97));
    STDtderiv1=std(Tderiv1(trace,3:97));
    eventtimes1(trace)=findevent(ToUse1(trace,start:end),STD1,Tderiv1(trace,start:end),STDtderiv1);
end
[~,I]=sort(eventtimes1);Order1=I;
ReorderedEventTimes1=sort(eventtimes1);
ReorderedSmoothNorm1=ToUse1(Order1,:);

% COUNT NUMBER OF RESPONSES, REORDER
Numberofresponses1=nnz(~isnan(ReorderedEventTimes1));
ind = ~isnan(ReorderedEventTimes1);ReorderedEventTimesResponses1=ReorderedEventTimes1(ind);
Responses1=NaN([Numberofresponses1 301]);
for trace=1:Numberofresponses1
    if ~isnan(ReorderedEventTimes1(trace))
        Responses1(trace,:)=ReorderedSmoothNorm1(trace,:);
    end
end

% Calculate mean, ste of responses(R) vs. all(A)
FmeanR1=mean(Responses1);FsteR1=ste(Responses1);
FmeanA1=mean(ReorderedSmoothNorm1);FsteA1=ste(ReorderedSmoothNorm1);

% CREATE SHIFTED RESPONSES OF EARLY RESPONSES ONLY, IE ALIGNED TO RESPONSE START
count = 0;
for trace = 1:size(Responses1)
    if ReorderedEventTimesResponses1(trace) < 50 %ie only if within first 5 seconds of stimulus
        count = count + 1;
        Shifted1(count,:) = aligntostart(Responses1(trace,:),ReorderedEventTimesResponses1(trace));
        TderivShifted1(count,:) = tderiv(Shifted1(count,:));
    end
end

%FIND RISE TIMES OF ALIGNED
RiseTime1 = [];
for trace = 1:size(Shifted1)
    P66= prctile(Shifted1(trace,47:97),66);
    P33= prctile(Shifted1(trace,47:97),33);
    for i=47:100
        if Shifted1(trace,i) > P66
            FrameP66 = i;
            break
        end
    end
    for i=47:100
        if Shifted1(trace,i) > P33
            FrameP33 = i;
            break
        end
    end
    RiseTime1(trace) = FrameP66-FrameP33;
end

% Cumulative Response Time Profile plots
EVENTTIMES1 = eventtimes1;EVENTTIMES1(isnan(EVENTTIMES1)) = inf; %change nonevents' event times to infinity
[cdf1,x1]=ecdf(EVENTTIMES1);
figure();plot(x1,cdf1,'b');xlim([0 50]);hold on;ylim([0 1]);hold on;


end

%---------------%
% SUB-FUNCTIONS %
%---------------%
function event = findevent(Trace,STD,tderiv,STDtderiv)
event = NaN;
for t=3:103 %Note: all t's are shifted by 2 as input
    if mean(Trace(t-2:t+10))>(2*STD) && mean(tderiv(t-2:t-1))>=(1*STDtderiv)
        event=t-2;
        break
    end
end
end

function derivative = tderiv(vector,window,fps)
%TDERIV Takes the time derivative of VECTOR using dt = WINDOW. WINDOW
%should be an odd number.  If not, the program will automatically add one
%to it.If no WINDOW is defined, it is assumed WINDOW = 1. If frames per
%second (FPS) is defined, then the output is in dX/dt (AU/sec).  If not,
%output is dX/dframe.
if nargin < 2
    window = 5;
end
if nargin < 3
    fps = 10;
end
derivative = NaN(size(vector));
if mod(window,2)==0
    window = window+1;
end
stepsize = floor(window/2);
for m=stepsize+1:length(vector)-stepsize
    derivative(m) = (vector(m+stepsize) - vector(m-stepsize))/window;
end
derivative = derivative*fps;
end

function std_error = ste(X,dim)
% Takes standard error of sample mean
if nargin < 2
    dim = 1;
end
if min(size(X)) == 1
    std_error = nanstd(X)./sqrt(length(X));
else    
    std_error = nanstd(X,0,dim)./sqrt(size(X,dim));
end
end

function cleaned = Cleanup(RawMatrix)
%Takes raw matrix and first transposes to proper direction
%then removes NaN traces
PulsesB = RawMatrix;
[t1, t2]=size(PulsesB);
PulsesB=PulsesB';
cleaned=[];
for m=1:size(PulsesB)
    if ~isnan(max(PulsesB(m,:),[],'includenan'))
        cleaned=[cleaned;PulsesB(m,:)];
    end
end
end

function shifted = aligntostart(Response,Time)
%takes the responses and shifts the trace to align to event start
%to look at shape of responses, regardless of when it begins
RealTime = Time + 100;%because eventtimes gives event time within stimulus pulse, not within entire video
shifted = Response(RealTime-50:RealTime+99);
end