function [ latTimes, lsqCurveTimes, lsqFittedCurve, lsqFittedParam, lsqCurveFitErr] = computeLatencies( MuaeData, MuaeTimes, lsqTimeRes )
% COMPUTE LATENCY function written by Demetrio Ferro demetrio.ferro@iit.it
% Based on (Roelfsema, Tolboom, Khayat, Neuron 2007).
% developed for (Ferro, van Kempen, Boyd, Panzeri, Thiele, PNAS 2021)
%
% The function computes Latency indices base don MUAes fitted by cumulative 
% gaussian function with 5 (degrees of freedom) parameters (mu, sigma, alpha, c, d).
% The best fit is assigned based on Least Squares (LS) of iterative estimation.
%
% Inputs
% MuaeData          [NxT]  Multi-channel (N channels) MUAes in time (T points) 
% MuaeTimes         [1xT]  Timestamps of MUAe recordings (T points)
% lsqTimeRes        [1x1]  Resolution of time series for Latency estimation
%
% Outputs
% latTimes          [1xN]  Latency indices at t*:={f(t*)=1/3*max(f(t)),t in MuaeTimes}
% lsqCurveTimes     [1x(T*lsqTimeRes)] Timestamps of fitted curves
% lsqFittedCurve    [Nx(T*lsqTimeRes)] Fitted curves at LS parameters
% lsqFittedParam    [Nx5]  LS parameters
% lsqFitErr         [1xN]  LS error of the fit

[N,T]=size(MuaeData);
latTimes=zeros(1,N);

lsqCurveTimes=linspace(MuaeTimes(1),MuaeTimes(end),ceil(T*lsqTimeRes));
lsqFittedCurve=zeros(N,ceil(T*lsqTimeRes));
lsqFittedParam=zeros(N,5);
lsqCurveFitErr=zeros(1,N);

% LsqParam0 = [mu, sigma, alpha, c, d];
% LsqParam0 is the starting point for the fitting routine;
lsqParam0=[0 0 0 0 0]; % starts with naive values for V1;

% Ex-Gaussian Function to fit
exGfunToFit = @(p,tdata)(p(5)*exp(p(1)*p(3)+0.5*p(2)^2*p(3)^2-p(3)*tdata)...
    .*normcdf(tdata,p(1)+p(2)^2*p(3),p(2))...
    +p(4)*normcdf(tdata,p(1),p(2)));

movMeanResolution=5;
mu0ValuesOrigVec=MuaeTimes;

options = optimoptions('lsqcurvefit','Display','off');

% LOOP OVER CHANNELS
for ch=1:N

    MuaeDataMovMean=movmean(MuaeData(ch,:),movMeanResolution);
    IntMuaeDataMM=interp1(MuaeTimes,MuaeDataMovMean,mu0ValuesOrigVec);
    DiffMuaeDataMM=diff([IntMuaeDataMM IntMuaeDataMM(end)]);

    % EMPIRICAL SETTING: TAKE ONLY POSITIVE SLOPES TO SPEED UP
    mu0ValuesVec=mu0ValuesOrigVec(DiffMuaeDataMM>=0); %.2*min(DiffMuaeDataMM));    
    
    lsqMuErrMovMean=-1*ones(1,length(mu0ValuesVec));
    lsqFittedParamMovMean=zeros(length(mu0ValuesVec),5);
    
    % ATTENTION: FIT THE MOVMEAN THEN USE THE MIN MSE PARAMETERS
    for muIdx=1:length(mu0ValuesVec)
        lsqParam0(1)=mu0ValuesVec(muIdx);
        try
        lsqFittedParamMovMean(muIdx,:)=lsqcurvefit(exGfunToFit,lsqParam0,MuaeTimes,MuaeDataMovMean,[],[],options);
        lsqMuErrMovMean(muIdx)=norm(feval(exGfunToFit,lsqFittedParamMovMean(muIdx,:),MuaeTimes)-MuaeData(ch,:));
        catch
            warning('lsqcurvefit does not converge for the parameters chosen!');
        end
    end
    % Remove non-tested mu parameters
    lsqMuErrMovMean(lsqMuErrMovMean==-1)=[];
    lsqFittedParamMovMean(lsqMuErrMovMean==-1)=[];
    
    [~,lsqMinMuErrIdxMovMean]=min(lsqMuErrMovMean);
    
    lsqFittedCurve(ch,:)=feval(exGfunToFit,lsqFittedParamMovMean(lsqMinMuErrIdxMovMean,:),lsqCurveTimes);
    lsqCurveFitErr(ch)=norm(interp1(lsqCurveTimes,lsqFittedCurve(ch,:),MuaeTimes)-MuaeData(ch,:));
    lsqFittedParam(ch,:)=lsqFittedParamMovMean(lsqMinMuErrIdxMovMean,:);
    
    latIndex=find(lsqFittedCurve(ch,:)>=1/3*max(lsqFittedCurve(ch,:)),1,'first');
    if ~isempty(latIndex)
        latTimes(ch)=lsqCurveTimes(latIndex);
    else
        latTimes(ch)=NaN;
    end
    
end
end