clear all
%%
%FOLDERS SELECTION

%Selecting assay folders containing data and avi files
%(Warning: make sure hdf data file is in the same folder as all_trials.avi file!)
files=uipickfiles('FilterSpec','/Users/stevenknafo1/Documents/MATLAB/Aequorin_data/Charlie_data/Aequorin/Mnx','Prompt','Select folders with data files');

%Assay parameters
 durationTrials=2500; %duration of each trial in ms
 binningTime=10; %binning time for photons count in ms

%TIMELINE ANALYSIS

for numFile=1:length(files)
    
    FilePath = files{1,numFile};
    [PathName,FolderName]=fileparts(FilePath);
    DataFile=fullfile(PathName,FolderName,[FolderName '_0.hdf']);
    
    %Importing the data file
    data = dataset('File',fullfile(DataFile),'delimiter',';','HeaderLines',3,'ReadVarNames',false);
    %Naming columns and creating variables
    data.Properties.VarNames = {'index', 'image', 'timestamp', 'shutter', 'PMT', 'stim1', 'stim2'};

    %Determine assay parameters
    assayPhotons=data(:,:);
    numTrials=round(length(data)/durationTrials);
    durationAssay=numTrials*durationTrials;
    %Create stimuli timeline
    stim=zeros(durationTrials,1);
    stim=assayPhotons.stim1+assayPhotons.stim2;
    stim(length(stim):durationAssay)=stim(length(stim));
    %Create frames timeline
    image=zeros(durationTrials,1);
    image=assayPhotons.image;
    image(length(image):durationAssay)=image(length(image));

    %Creating cell arrays 
    trial_stimulus = cell (numTrials,1);
    trial_stim_start = cell(numTrials,1);
    trial_frame_start = cell(numTrials,1);
    trial_frame_end = cell(numTrials,1);
    trial_frame = cell(numTrials,1);
    assay_frame_stim_start= cell(numTrials,1);
    trial_frame_stim_start = cell(numTrials,1);

    %Filling the cell arrays with discrete values for each trial
    for i= 1:numTrials;
        Limit_inf = (i-1)*2500+1;
        Limit_sup = (Limit_inf+2498);
        trial_stimulus {i} = stim (Limit_inf : Limit_sup);
        trial_frame{i}=image(Limit_inf : Limit_sup);
        trial_frame_start{i}=image(Limit_inf);
        trial_frame_end{i}=image(Limit_sup);
        trial_stim_start{i} = find(trial_stimulus{i},1,'first');
        assay_frame_stim_start{i}=trial_frame{i}(trial_stim_start{i},1);
        trial_frame_stim_start{i}=assay_frame_stim_start{i}-trial_frame_start{i};
    end

    Frame_Stim = [trial_frame_start,trial_frame_end,trial_stim_start,assay_frame_stim_start,trial_frame_stim_start];
    Timeline = cell2dataset (Frame_Stim,'ReadVarNames',false, 'VarNames', {'trial_frame_start','trial_frame_end', 'trial_stim_start', 'assay_frame_stim_start', 'trial_frame_stim_start'});
    save(fullfile(strrep(DataFile,'.hdf','_Quantitative')),'Timeline');

end


clearvars -except files numTrials durationTrials durationAssay binningTime

%%
% BIOLUMINESCENCE ANALYSIS

for numFile=1:length(files);

    FilePath = files{1,numFile};
    [PathName,FolderName]=fileparts(FilePath);
    DataFile=fullfile(PathName,FolderName,[FolderName '_0.hdf']);
    load(fullfile(strrep(DataFile,'.hdf','_Quantitative')));
    
    %Importing the data file
    data = dataset('File',fullfile(DataFile),'delimiter',';','HeaderLines',3,'ReadVarNames',false);
    %Naming columns and creating variables
    data.Properties.VarNames = {'index', 'image', 'timestamp', 'shutter', 'PMT', 'stim1', 'stim2'};
    %Determine number of photons per ms
    data.photons_ms = diff(data.PMT);

    %Keep only assay lines
    assayPhotons=data.photons_ms(:,:);
    assayPhotons(length(assayPhotons):durationAssay)=assayPhotons(length(assayPhotons));
    assayPhotons=assayPhotons(1:durationAssay);
    %Sum the number of photons according to binning time
    A=reshape(assayPhotons,binningTime,length(assayPhotons)/binningTime);
    photons_binned = sum(reshape(assayPhotons,binningTime,length(assayPhotons)/binningTime));

    %Replacing PMT artefact outliers by median value + add first value
    outlier_threshold=100;
    photons_binned(photons_binned>outlier_threshold)=median(median(photons_binned));

    %Creating cell arrays 
    photons = cell(numTrials,1);

    %Filling the cell arrays with discrete values for each trial
    for i= 1:numTrials;
        Limit_inf = (i-1)*(durationTrials/binningTime)+1;
        Limit_sup = (i)*(durationTrials/binningTime);
        photons {i} = photons_binned(Limit_inf:Limit_sup);
    end

    %Grouping all trials within one dataset
    B_all_trials = photons;

    %Filtering signal
    filtered_signal=zeros((durationTrials/binningTime),numTrials);
    B_per_trial=cell(numTrials,durationTrials/binningTime);
    windowSize = 10;
    for i=1:numTrials;
         B_per_trial{i,1}=B_all_trials{i,1};
         filtered_signal(:,i)=filter(ones(1,windowSize)/windowSize,1,B_per_trial{i,1});
    end;

    %Linear fit of noise (based on first 500ms of each trial)
    coeff_noise=zeros(2,numTrials);
    cum_noise_only_trial=zeros((durationTrials/binningTime),numTrials);
    cum_signal_only_trial=zeros((durationTrials/binningTime));
    signal_only_trial=zeros((durationTrials/binningTime)-1,numTrials);
    time=1:binningTime:2500;
    for i=1:numTrials, 
        B_baseline=cumsum(filtered_signal(1:50,i));
        T_baseline=rot90(1:10:500,3);
        linear_noise=fit(T_baseline,B_baseline,'poly1');
        coeff_noise(:,i)=coeffvalues(linear_noise);
        cum_noise_only_trial(:,i)=coeff_noise(1,i)*time+coeff_noise(2,i);
        cum_signal_only_trial(:,i)=(cumsum(filtered_signal(:,i))-cum_noise_only_trial(:,i));
        signal_only_trial(:,i)=diff(cum_signal_only_trial(:,i));
    end

    %Quantitative analysis of bioluminescence signals
    signal_start_time = cell (numTrials,1);
    signal_end_time = cell (numTrials,1);
    variationSignal=cell(numTrials,1);
    signal = cell (numTrials,1);
    latency = cell (numTrials,1);
    amplitude = cell (numTrials,1);
    integrale = cell (numTrials,1);
    duration = cell (numTrials,1);
    time2peak = cell (numTrials,1);
    decay = cell (numTrials,1);

    threshold_start=0.4; N_start=3; 
    threshold_end=0.2; N_end=3;

    for i=1:numTrials;

        variationSignal{i,1}=diff(signal_only_trial(:,i)); variationSignal{i,1}(durationTrials/binningTime)=0;
        
        start_signal{i,1}=find(abs(variationSignal{i,1})>threshold_start); start_signal{i,1}=rot90(start_signal{i,1});
        x=diff(start_signal{i,1})==1;
        f=find([false,x]~=[x,false]);
        g=find(f(2:2:end)-f(1:2:end-1)>=N_start,1,'first');
        signal_start_time{i,1} = start_signal{i,1}(f(2*g-1)); % First value followed by >=N consecutive numbers
        
        variationSignal{i,1}=variationSignal{i,1}(signal_start_time{i,1}:end);
        end_signal{i,1}=find(abs(variationSignal{i,1})<threshold_end); end_signal{i,1}=rot90(end_signal{i,1});
        y=diff(end_signal{i,1})==1;
        u=find([false,y]~=[y,false]);
        h=find(u(2:2:end)-u(1:2:end-1)>=N_end,1,'first');
        
        empty_cells = cellfun('isempty',signal_start_time);
        signal_start_time(empty_cells) = {0};
        empty_cells = cellfun('isempty',signal_end_time);
        signal_end_time(empty_cells) = {0};
        signal_end_time{i,1} = end_signal{i,1}(u(2*h-1))+signal_start_time{i,1};

        signal{i}= signal_only_trial(signal_start_time{i}:signal_end_time{i},i);
        if signal_start_time{i}==0;
           latency{i}=0;
        else
        latency{i}=signal_start_time{i}*10-Timeline.trial_stim_start(i); %*10 pour avoir la valeur en ms
        end
        [amplitude{i},time2peak{i}] = max (signal{i});
        duration{i} = length (signal{i})*binningTime; 
        time2peak{i}=time2peak{i}*binningTime; 
        decay{i} = duration{i}-time2peak{i};
        integrale{i} = trapz (signal{i});  
    end

    %Display bioluminescence signals for each trial
        figure('units','normalized','outerposition',[0 0 1 1]);
        X_axis_top=max(max(signal_only_trial(:,:)));
    for i=1:numTrials, 
       subplot(1*numTrials/10,10,i),plot(signal_only_trial(:,i),'g', 'LineWidth',1), title(['trial ',num2str(i)])
       if signal_start_time{i,1} ~= 0 
       hold on, plot(signal_start_time{i,1},signal_only_trial(signal_start_time{i,1},i),'og'); end
       if signal_end_time{i,1} ~= 0
       hold on, plot(signal_end_time{i,1},signal_only_trial(signal_end_time{i,1},i),'og'); end
       hold on, line([Timeline.trial_stim_start(i)/10,Timeline.trial_stim_start(i)/10],[-2,X_axis_top],'color','r')
       axis([0 (durationTrials/binningTime) -2 X_axis_top]); set(gca,'FontSize',14)
    end
    %%
    %Manually determine signals that were not properly detected
    feedback=input('are there trials to modify?');
    feedback_discard=zeros(numTrials,1);
    for i=1:length(feedback);
        j=feedback(i);
        figure('units','normalized','outerposition',[0 0 1 1]),plot(signal_only_trial(:,j),'g', 'LineWidth',2)
        hold on, line([Timeline.trial_stim_start(j)/10,Timeline.trial_stim_start(j)/10],[-2,X_axis_top],'color','r')
        axis([0 (durationTrials/binningTime) -2 X_axis_top]);
        [x,y]=ginput(1); signal_start_time{j,1}=round(x); 
        [x,y]=ginput(1); signal_end_time{j,1}=round(x);
        signal{j}= signal_only_trial(signal_start_time{j}:signal_end_time{j},j);
        if isempty(x)==1, feedback_discard(j,1)=1; signal_start_time{j}=0; end
        close;
        if signal_start_time{j}==0; latency{j}=0;
        else latency{j}=signal_start_time{j}*10-Timeline.trial_stim_start(j); end
        [amplitude{j},time2peak{j}] = max (signal{j});
        duration{j} = length (signal{j})*binningTime; 
        time2peak{j}=time2peak{j}*binningTime; 
        decay{j} = duration{j}-time2peak{j};
        integrale{j} = trapz (signal{j}); 
    end 
    close;
    %%
    %Calculate rising and decay tau coefficients
    rising_coeff=cell (numTrials,1); decay_coeff=cell (numTrials,1);
    rising_signal=cell (numTrials,1); decay_signal=cell (numTrials,1);
    rising_time=cell (numTrials,1); decay_time=cell (numTrials,1);
    rising_fit=cell (numTrials,1); decay_fit=cell (numTrials,1);
    rising_slope=cell (numTrials,1); decay_tau=cell (numTrials,1);
    rising_gof=cell (numTrials,1); decay_gof=cell (numTrials,1);
    rising_rsquare=cell (numTrials,1); decay_rsquare=cell (numTrials,1);
    for i=1:numTrials,
        if signal_start_time{i}==0; rising_slope{i}=0; decay_tau{i}=0;
        else
        rising_signal{i}=signal_only_trial(signal_start_time{i}:(signal_start_time{i}+time2peak{i}/10),i);
        decay_signal{i}=signal_only_trial((signal_start_time{i}+time2peak{i}/10):signal_end_time{i},i);
        rising_time{i}=rot90(signal_start_time{i}:1:(signal_start_time{i}+time2peak{i}/10),3);
        decay_time{i}=rot90((signal_start_time{i}+time2peak{i}/10):1:signal_end_time{i},3);
        [rising_fit{i},rising_gof{i}]=fit(rising_time{i},rising_signal{i},'poly1'); rising_coeff{i}=coeffvalues(rising_fit{i}); 
        rising_slope{i}= rising_coeff{i,1}(1); rising_rsquare{i}=rising_gof{i}.rsquare;
        [decay_fit{i},decay_gof{i}]=fit(decay_time{i},decay_signal{i},'exp1'); decay_coeff{i}=coeffvalues(decay_fit{i}); 
        decay_tau{i}=-1/decay_coeff{i,1}(2)*10; decay_rsquare{i}=decay_gof{i}.rsquare;
%         subplot(2*numTrials/10,10,i),plot(rising_fit{i},rising_time{i},rising_signal{i}),legend('off');
%         subplot(2*numTrials/10,10,10+i),plot(decay_fit{i},decay_time{i},decay_signal{i}),legend('off');
        end
    end
    %%
    %Create quantitative results dataset
    parameters = [signal_start_time,signal_end_time,latency,amplitude, duration, time2peak, decay, integrale, rising_slope,rising_rsquare,decay_tau,decay_rsquare];
    empty_cells = cellfun('isempty',parameters); parameters(empty_cells) = {0};
    %%
    %Replace values with NaN in discarded trials
    feedback_discard=input('are there trials to discard?');
    for i=1:length(feedback_discard);
        for j=1:length(parameters(i,:)); parameters{feedback_discard(i),j}=NaN; end
    end
    %%
    Bioluminescence_quantitative = cell2dataset (parameters, 'ReadVarNames',false, 'VarNames', {'signal_start','signal_end','latency','amplitude','duration', 'time2peak', 'decay', 'integrale','rising_slope','rising_rsquare','decay_tau','decay_rsquare'}); 
     
   %Calculate normalized B parameters with both z-scores and Min-Max
   %Exclude discarded trials
   Bioluminescence_quantitative_nonNaN=dataset2cell(Bioluminescence_quantitative);
   Bioluminescence_quantitative_nonNaN(1,:)=[];
   Bioluminescence_quantitative_nonNaN=cell2mat(Bioluminescence_quantitative_nonNaN);
   Bioluminescence_quantitative_nonNaN(any(isnan(Bioluminescence_quantitative_nonNaN(:,3)),2),:)=[];
   Bioluminescence_quantitative_nonNaN=mat2dataset(Bioluminescence_quantitative_nonNaN, 'VarNames', {'signal_start','signal_end','latency','amplitude','duration', 'time2peak', 'decay', 'integrale','rising_slope','rising_rsquare','decay_tau','decay_rsquare'});

   for i=1:numTrials
        %Z-score and Min-Max normalization for B_amplitude
        Bioluminescence_quantitative.B_amp_Z(i)=(Bioluminescence_quantitative.amplitude(i)-mean(Bioluminescence_quantitative_nonNaN.amplitude))/std(Bioluminescence_quantitative_nonNaN.amplitude);
        Bioluminescence_quantitative.B_amp_N(i)=(Bioluminescence_quantitative.amplitude(i)-min(Bioluminescence_quantitative.amplitude))/(max(Bioluminescence_quantitative.amplitude)-min((Bioluminescence_quantitative.amplitude))); 
        
        %Z-score and Min-Max normalization for B_rising_slope
        Bioluminescence_quantitative.B_rise_Z(i)=(Bioluminescence_quantitative.rising_slope(i)-mean(Bioluminescence_quantitative_nonNaN.rising_slope))/std(Bioluminescence_quantitative_nonNaN.rising_slope);
        Bioluminescence_quantitative.B_rise_N(i)=(Bioluminescence_quantitative.rising_slope(i)-min(Bioluminescence_quantitative.rising_slope))/(max(Bioluminescence_quantitative.rising_slope)-min((Bioluminescence_quantitative.rising_slope)));  
        
        %Z-score and Min-Max normalization for B_decay_tau
        Bioluminescence_quantitative.B_decay_Z(i)=(Bioluminescence_quantitative.decay_tau(i)-mean(Bioluminescence_quantitative_nonNaN.decay_tau))/std(Bioluminescence_quantitative_nonNaN.decay_tau);
        Bioluminescence_quantitative.B_decay_N(i)=(Bioluminescence_quantitative.decay_tau(i)-min(Bioluminescence_quantitative.decay_tau))/(max(Bioluminescence_quantitative.decay_tau)-min((Bioluminescence_quantitative.decay_tau)));
    end
    
    %Group signal data per trial as: raw ; filtered ; start ; end ; max ; locs
    signal_data=cell(numTrials,6);
    for i=1:numTrials;
        signal_data{i,1}=B_per_trial{i,1};
        signal_data{i,2}=signal_only_trial(:,i);
        signal_data{i,3}=signal_start_time{i,1};
        signal_data{i,4}=signal_end_time{i,1};
        signal_data{i,5}=amplitude{i};
        signal_data{i,6}=time2peak{i}+signal_start_time{i,1};
    end
    
    save(fullfile(strrep(DataFile,'.hdf','_Quantitative')),'Bioluminescence_quantitative','signal_data','-append');

end

clearvars -except files durationTrials binningTime

%%
%KINEMATICS ANALYSIS

%Read movies and place tail base points
numFile=length(files);
base_point=cell(numFile,4);

for numFile=1:length(files);
    FilePath = files{1,numFile};
    movie_path=fullfile(FilePath,'All_trials.avi');
    all_trials_movie(numFile)=VideoReader(movie_path);
    a1=figure;
    S_image=read(all_trials_movie(numFile),1);
    imshow(S_image);
    %Get tail base coordinates
    tail_base=impoint(gca,[]);
    tail_base=wait(tail_base);
    hold on;
    scatter(tail_base(1),tail_base(2),3,'filled','g');
    base_point{numFile,1}=tail_base(1);
    base_point{numFile,2}=tail_base(2);
    %Get tail tip coordinates
    tail_tip=impoint(gca,[]);
    tail_tip=wait(tail_tip);
    hold on;
    scatter(tail_base(1),tail_base(2),3,'filled','g');
    base_point{numFile,3}=tail_tip(1);
    base_point{numFile,4}=tail_tip(2);
    
    close(a1)
end

%Tail tracking and kinematics quantitative analysis
for numFile=1:length(files);

    %Setting files
    FilePath = files{1,numFile};
    [PathName,FolderName]=fileparts(FilePath);
    DataFile=fullfile(PathName,FolderName,[FolderName '_0.hdf']);
    
    load(fullfile(strrep(DataFile,'.hdf','_Quantitative')));
    
    [PathName,FolderName]=fileparts(FilePath);
    DataFile=fullfile(PathName,FolderName,[FolderName '_0.hdf']);
    
    %Setting tracking parameters
    number_segments=10;
    segment_length=(base_point{numFile,1}-base_point{numFile,3})/number_segments;
    
    %Setting assay parameters
    numTrials=length(Timeline);
    durationAssay=numTrials*durationTrials;
    numFrames=Timeline.trial_frame_end(numTrials);
    
    %Lance le tail tracker pour récupérer les coordonnées des segments
    [TAIL]=Tail_tracker(all_trials_movie,segment_length,number_segments,numFrames,numFile,base_point);
    
    %Replace Nan by preceding value in TAIL
        for k=1:2,
            for i=1:size(TAIL,1),
                for j=2:size(TAIL,2),
                    if (isnan(TAIL(i,j,k)))
                        TAIL(i,j,k)=TAIL(i,j-1,k);
                    end
                end
            end
        end
        
    %Calculate alpha angle between tail base and tip
     alpha_rad=zeros(length(TAIL),1);
     alpha=zeros(length(TAIL),1);
     rapport=zeros(length(TAIL),1);
     for frame=1:length(TAIL)
         
         Ytip=TAIL(frame,11,1);
         Ybase=TAIL(frame,1,1);
         Xtip=TAIL(frame,11,2);
         Xbase=TAIL(frame,1,2);
         
         if Ytip-Ybase<=0 %flexion de la queue sur la gauche: alpha>0
             if Xtip-Xbase<=0 %flexion < 90°: Xa-Xo<0
                rapport(frame,1)=(Ybase-Ytip)/(Xbase-Xtip);
                alpha_rad(frame,1)=atan(rapport(frame,1));
                alpha(frame,1)=alpha_rad(frame)*57.29578; 
             elseif Xtip-Xbase>0 %flexion > 90°: Xa-Xo>0
                rapport(frame,1)=(Ybase-Ytip)/(Xtip-Xbase);
                alpha_rad(frame,1)=atan(rapport(frame,1));
                alpha(frame,1)=180-(alpha_rad(frame)*57.29578); 
             end
         elseif Ytip-Ybase>0% flexion de la queue sur la droite: alpha<0
            if Xtip-Xbase<0 %flexion < 90°: Xa-Xo<0
                rapport(frame,1)=(Ybase-Ytip)/(Xtip-Xbase);
                alpha_rad(frame,1)=-(atan(rapport(frame,1)));
                alpha(frame,1)=alpha_rad(frame)*57.29578; 
            elseif Xtip-Xbase>=0 %flexion > 90°: Xa-Xo>0
                rapport(frame,1)=(Ytip-Ybase)/(Xtip-Xbase);
                alpha_rad(frame,1)=-(atan(rapport(frame,1)));
                alpha(frame,1)=-(180+alpha_rad(frame)*57.29578); 
            end
         end
     end
    
    %Apply median filtering to alpha angle
    alpha_filtered=medfilt2(alpha,[10 1]);
 
     %Split angle data for each trial according to timeline limits
    trial_alpha=cell(numTrials,2);
     for i=1:numTrials;
        trial_alpha{i,1}=alpha(Timeline.trial_frame_start(i):Timeline.trial_frame_end(i));
        trial_alpha{i,1}(length(trial_alpha{i,1}):durationTrials)=trial_alpha{i,1}(length(trial_alpha{i,1}));
        trial_alpha{i,2}=alpha_filtered(Timeline.trial_frame_start(i):Timeline.trial_frame_end(i));
        trial_alpha{i,2}(length(trial_alpha{i,2}):durationTrials)=trial_alpha{i,2}(length(trial_alpha{i,2}));
     end
     
     %Replace Nan by interpolated value in trial_alpha
    for i = 1:numTrials
        times = 1:length(trial_alpha{i,2});
        mask =  ~isnan(trial_alpha{i,2});
        trial_alpha{i,3} = trial_alpha{i,2};
        trial_alpha{i,3}(~mask) = interp1(times(mask), trial_alpha{i,3}(mask), times(~mask));
    end

    %substract alpha baseline
    for i = 1:numTrials
      alpha_baseline=mean(trial_alpha{i,3}(1:400));
      trial_alpha{i,3}=trial_alpha{i,3}-alpha_baseline;
    end

    %Find start and end of the escape response for each trial
    variationAlpha=cell(numTrials,1);
    movement=cell(numTrials,1);
    start_frame = cell (numTrials,1);
    end_frame = cell (numTrials,1);
    threshold_start=0.08;%Variation threshold of alpha in degree between two consecutive frames
    N_start=3; %Required number of consecutive frames above threshold to detect movement  
    threshold_end=0.1; N_end=20;
    
     for i=1:numTrials;
        variationAlpha{i,1}=diff(trial_alpha{i,3}); variationAlpha{i,1}(2500)=0;
        movement_start{i,1}=find(abs(variationAlpha{i,1})>threshold_start); movement_start{i,1}=rot90(movement_start{i,1});
        x=diff(movement_start{i,1})==1;
        f=find([false,x]~=[x,false]);
        g=find(f(2:2:end)-f(1:2:end-1)>=N_start,1,'first');
        % First value followed by >=N consecutive numbers
        start_frame{i,1} = movement_start{i,1}(f(2*g-1)); 
        %only consider start_frame that are after frame_stim
        shift=0;
        while start_frame{i,1} < Timeline.trial_frame_stim_start(i);
              start_frame{i,1}=start_frame{i,1}+shift;
              shift=shift+1;
            continue
        end
         

        %%
        variationAlpha{i,1}=variationAlpha{i,1}(start_frame{i,1}:end);
        movement_end{i,1}=find(abs(variationAlpha{i,1})<threshold_end); movement_end{i,1}=rot90(movement_end{i,1});
        y=diff(movement_end{i,1})==1;
        u=find([false,y]~=[y,false]);
        h=find(u(2:2:end)-u(1:2:end-1)>=N_end,1,'first');
        empty_cells = cellfun(@isempty,start_frame); start_frame(empty_cells) = {0};
        end_frame{i,1} = movement_end{i,1}(u(2*h-1))+start_frame{i,1};
        empty_cells = cellfun(@isempty,end_frame); end_frame(empty_cells) = {0};
     end
   
    
    %Compute quantitative kinematics analysis for each escape
    ER = cell (numTrials,1);
    latency = cell (numTrials,1);
    amplitude = cell (numTrials,1);
    duration = cell (numTrials,1);
    N_cycles = cell (numTrials,1);
    frequency = cell (numTrials,1);
    pks= cell (numTrials,1);
    locs=cell (numTrials,1);


    for i=1:numTrials;
        if end_frame{i,1}-start_frame{i,1} > 3 ;
            ER{i}=trial_alpha{i,3}(start_frame{i,1}:end_frame{i,1});
            latency{i}=start_frame{i}-Timeline.trial_frame_stim_start(i);
            amplitude{i}=max(abs(ER{i}));
            duration{i}=length(ER{i});
            
            %Find all extrema, both max and min
            [pks_max{i},locs_max{i}] = findpeaks(ER{i});
            DataInv = 1.01*max(ER{i})-ER{i};
            [pks_min{i},locs_min{i}] = findpeaks(DataInv);
            pks_min{i} = ER{i}(locs_min{i});
            %Group max and min together and sort them by time
            locs_max{i}=locs_max{i}+start_frame{i};
            locs_min{i}=locs_min{i}+start_frame{i};
            locs{i}=vertcat(locs_min{i},locs_max{i});locs{i}=sort(locs{i});
            %Select only extrema that are located > 5 ms from each other
            for k=1:length(locs{i,1})-1
                if locs{i,1}(k+1)-locs{i,1}(k)<2
                locs{i,1}(k)=0;
                end
            end
            locs{i,1}(all(~locs{i,1},2),:)=[];
            %Among those, select only extrema with signal > 0.5 from n-10 signal
            for k=1:length(locs{i,1})
                if abs(trial_alpha{i,3}(locs{i}(k))-trial_alpha{i,3}(locs{i}(k)-5))<1
                locs{i,1}(k)=0;
                end
            end
            locs{i,1}(all(~locs{i,1},2),:)=[];
            %Determine signal values corresponding to selected locations
            pks{i}=ER{i}(locs{i}-start_frame{i});
            
            N_cycles{i}=length(locs{i})/2;
            frequency{i}=N_cycles{i}/duration{i}*1000;
        else
            ER{i}=0; duration{i}=0; latency{i}=0;
        end
    end
   
       %Display alpha angle for each trial
        figure('units','normalized','outerposition',[0 0 1 1]); title('Kinematics_analysis')
    for i=1:numTrials;
        subplot(1*numTrials/10,10,i),plot(ER{i}), title(['trial ',num2str(i)])
        if start_frame{i,1}~=0 && end_frame{i,1}~=0;
        hold on, subplot(1*numTrials/10,10,i),plot(locs{i}-start_frame{i},ER{i}(locs{i}-start_frame{i}),'^k')
        end
        AX_top=round(max(Kinematics_quantitative.amplitude));
        if AX_top > 0; axis([0 length(ER{i}) -AX_top AX_top]);
        else axis([0 length(ER{i}) -180 180]);
        end
         set(gca,'FontSize',14)
        if numTrials==30; set(gcf,'PaperUnits','centimeters','PaperPosition',[0,0,35,35])
            else if numTrials==10; set(gcf,'PaperUnits','centimeters','PaperPosition',[0,0,35,5])
            end
        end
    end
   
    %Manually determine ER that were not properly detected
    feedback=input('are there trials to modify?');
    for i=1:length(feedback);
        j=feedback(i);
        figure('units','normalized','outerposition',[0 0 1 1]), plot(trial_alpha{j,3},'b', 'LineWidth',1)
        hold on, line([Timeline.trial_frame_stim_start(j),Timeline.trial_frame_stim_start(j)],[-180;180],'color','r')
        axis([0 durationTrials -180 180]);
        [x,y]=ginput(1); start_frame{j,1}=round(x);
        [x,y]=ginput(1); end_frame{j,1}=round(x);
        [x,y]=ginput; pks{j}=round(y); locs{j}=round(x);
        if isempty(x)==1, feedback_discard(j,1)=1;end
        close;

        ER{j}=trial_alpha{j,3}(start_frame{j,1}:end_frame{j,1});
        latency{j}=start_frame{j}-Timeline.trial_frame_stim_start(j);
        amplitude{j}=max(abs(ER{j}));
        duration{j}=length(ER{j});
        empty_cells = cellfun('isempty',ER);
        ER(empty_cells)={0};
        N_cycles{j}=length(locs{j})/2;
        frequency{j}=N_cycles{j}/duration{j}*1000;
    end
    
 
    %Classify ER as: 0 = no movement; 1 = escape only, 2 = slow swim only, 3 = escape+slow swim
    classification=cell(numTrials,1);
    for i=1:numTrials;
        if N_cycles{i}==(0.5||1);
            classification{i}=1;
        elseif isempty(N_cycles{i})==1;
            classification{i}=0;
            elseif N_cycles{i} > 1 && amplitude{i} < 45;
                classification{i}=2; 
                else N_cycles{i} > 1 && amplitude{i} > 45;
                classification{i}=3;
        end
    end

    
    %Group alpha data per trial as: raw ; filtered ; start ; end ; max ; locs
    alpha_data=cell(numTrials,6);
    for i=1:numTrials;
        alpha_data{i,1}=trial_alpha{i,1};
        alpha_data{i,2}=trial_alpha{i,3};
        alpha_data{i,3}=start_frame{i};
        alpha_data{i,4}=end_frame{i};
        alpha_data{i,5}=pks{i};
        alpha_data{i,6}=locs{i};
    end
    
    parameters = [start_frame, end_frame, latency, amplitude, duration, N_cycles, frequency, classification];
    empty_cells = cellfun('isempty',parameters); parameters(empty_cells) = {0};
    
    %Replace values with NaN in discarded trials (ex: struggles, non-canonical escapes)
    feedback_discard=input('are there trials to discard?');
    for i=1:length(feedback_discard);
        for j=1:length(parameters(i,:)); parameters{feedback_discard(i),j}=NaN; end
    end
   
    Kinematics_quantitative = cell2dataset (parameters, 'ReadVarNames',false, 'VarNames', {'start_frame','end_frame','latency','amplitude','duration','N_cycles','frequency','classification'});
   %% 
    %Calculate instantaneous frequency and angular velocity for movement
    for i=1:numTrials 
        if length(alpha_data{i,6})> 2
            %instantaneous frequency and velocity for first bend
            alpha_data{i,7}(1)=1000*1/(2*(alpha_data{i,6}(1)-alpha_data{i,3}));
            alpha_data{i,8}(1)=abs(alpha_data{i,5}(1))+(abs(alpha_data{i,2}(alpha_data{i,3})))/(alpha_data{i,6}(1)-alpha_data{i,3});
            for k=2:length(alpha_data{i,6})-1
                %instantaneous frequency for following bends
                alpha_data{i,7}(k)=1000*1/(2*(alpha_data{i,6}(k+1)-alpha_data{i,6}(k)));
                %Adding instantaneous frequency mean and max to dataset
                Kinematics_quantitative.K_instfqcy_mean(i)=mean(alpha_data{i,7});
                Kinematics_quantitative.K_instfqcy_max(i)=max(alpha_data{i,7});
                %instantaneous angular velocity for following bends
                alpha_data{i,8}(k)=abs(alpha_data{i,5}(k+1))+(abs(alpha_data{i,5}(k)))/(alpha_data{i,6}(k+1)-alpha_data{i,6}(k)); 
            end
        end
            
        if length(alpha_data{i,6})==2
            %instantaneous frequency for first and second bend
            alpha_data{i,7}(1)=1000*1/(2*(alpha_data{i,6}(1)-alpha_data{i,3}));
            alpha_data{i,7}(2)=1000*1/(2*(alpha_data{i,6}(2)-alpha_data{i,6}(1)));
            %instantaneous angular velocity for first bend
            alpha_data{i,8}(1)=abs(alpha_data{i,5}(1))+(abs(alpha_data{i,2}(alpha_data{i,3})))/(alpha_data{i,6}(1)-alpha_data{i,3});
            alpha_data{i,8}(2)=abs(alpha_data{i,5}(2))+(abs(alpha_data{i,5}(1)))/(alpha_data{i,6}(2)-alpha_data{i,6}(1));
        end
        if length(alpha_data{i,6})==1 
                alpha_data{i,7}(1)=1000*1/(2*(alpha_data{i,6}(1)-alpha_data{i,3}));
                alpha_data{i,8}(1)=abs(alpha_data{i,5}(1))+(abs(alpha_data{i,2}(alpha_data{i,3})))/(alpha_data{i,6}(1)-alpha_data{i,3});
        end
        
        if isempty(alpha_data{i,6})==1
                Kinematics_quantitative.K_instfqcy_mean(i)=0;
                Kinematics_quantitative.K_instfqcy_max(i)=0;
                Kinematics_quantitative.K_instvlcty_mean(i)=0;
                Kinematics_quantitative.K_instvlcty_max(i)=0;
        else
        	Kinematics_quantitative.K_instfqcy_mean(i)=mean(alpha_data{i,7});
            Kinematics_quantitative.K_instfqcy_max(i)=max(alpha_data{i,7});
            Kinematics_quantitative.K_instvlcty_mean(i)=mean(alpha_data{i,8});
            Kinematics_quantitative.K_instvlcty_max(i)=max(alpha_data{i,8});
        end
        
    end
    
   
   %% 
   %Calculate kinematics normalized parameters per fish 
   
   %Exclude discarded trials
   Kinematics_quantitative_nonNaN=dataset2cell(Kinematics_quantitative);
   Kinematics_quantitative_nonNaN(1,:)=[];
   Kinematics_quantitative_nonNaN=cell2mat(Kinematics_quantitative_nonNaN);
   Kinematics_quantitative_nonNaN(any(isnan(Kinematics_quantitative_nonNaN(:,7)),2),:)=[];
   Kinematics_quantitative_nonNaN=mat2dataset(Kinematics_quantitative_nonNaN, 'VarNames', {'start_frame','end_frame','latency','amplitude','duration','N_cycles','frequency','classification','K_instfqcy_mean','K_instfqcy_max','K_instvlcty_mean','K_instvlcty_max'});

   
   for i=1:numTrials 
        %Normalized amplitude
        Kinematics_quantitative.K_amp_Z(i)=(Kinematics_quantitative.amplitude(i)-mean(Kinematics_quantitative_nonNaN.amplitude))/std(Kinematics_quantitative_nonNaN.amplitude);
        Kinematics_quantitative.K_amp_N(i)=(Kinematics_quantitative.amplitude(i)-min(Kinematics_quantitative.amplitude))/(max(Kinematics_quantitative.amplitude)-min(Kinematics_quantitative.amplitude));
        %Normalized duration
        Kinematics_quantitative.K_duration_Z(i)=(Kinematics_quantitative.duration(i)-mean(Kinematics_quantitative_nonNaN.duration))/std(Kinematics_quantitative_nonNaN.duration);
        Kinematics_quantitative.K_duration_N(i)=(Kinematics_quantitative.duration(i)-min(Kinematics_quantitative.duration))/(max(Kinematics_quantitative.duration)-min(Kinematics_quantitative.duration));
        %Normalized instantaneous mean frequency
        Kinematics_quantitative.K_instfqcy_mean_Z(i)=(Kinematics_quantitative.K_instfqcy_mean(i)-mean(Kinematics_quantitative_nonNaN.K_instfqcy_mean))/std(Kinematics_quantitative_nonNaN.K_instfqcy_mean);
        Kinematics_quantitative.K_instfqcy_mean_N(i)=(Kinematics_quantitative.K_instfqcy_mean(i)-min(Kinematics_quantitative.K_instfqcy_mean))/(max(Kinematics_quantitative.K_instfqcy_mean)-min(Kinematics_quantitative.K_instfqcy_mean));
        %Normalized instantaneous max frequency
        Kinematics_quantitative.K_instfqcy_max_Z(i)=(Kinematics_quantitative.K_instfqcy_max(i)-mean(Kinematics_quantitative_nonNaN.K_instfqcy_max))/std(Kinematics_quantitative_nonNaN.K_instfqcy_max);
        Kinematics_quantitative.K_instfqcy_max_N(i)=(Kinematics_quantitative.K_instfqcy_max(i)-min(Kinematics_quantitative.K_instfqcy_max))/(max(Kinematics_quantitative.K_instfqcy_max)-min(Kinematics_quantitative.K_instfqcy_max));
        %Normalized instantaneous mean angular velocity
        Kinematics_quantitative.K_instvlcty_mean_Z(i)=(Kinematics_quantitative.K_instvlcty_mean(i)-mean(Kinematics_quantitative_nonNaN.K_instvlcty_mean))/std(Kinematics_quantitative_nonNaN.K_instvlcty_mean);
        Kinematics_quantitative.K_instvlcty_mean_N(i)=(Kinematics_quantitative.K_instvlcty_mean(i)-min(Kinematics_quantitative.K_instvlcty_mean))/(max(Kinematics_quantitative.K_instvlcty_mean)-min(Kinematics_quantitative.K_instvlcty_mean));
        %Normalized instantaneous max angular velocity
        Kinematics_quantitative.K_instvlcty_max_Z(i)=(Kinematics_quantitative.K_instvlcty_max(i)-mean(Kinematics_quantitative_nonNaN.K_instvlcty_max))/std(Kinematics_quantitative_nonNaN.K_instvlcty_max);
        Kinematics_quantitative.K_instvlcty_max_N(i)=(Kinematics_quantitative.K_instvlcty_max(i)-min(Kinematics_quantitative.K_instvlcty_max))/(max(Kinematics_quantitative.K_instvlcty_max)-min(Kinematics_quantitative.K_instvlcty_max));

   end
    
    save(fullfile(strrep(DataFile,'.hdf','_Quantitative')),'Kinematics_quantitative','DataFile','TAIL','alpha_data','-append');  
%%
 % QUALITATIVE RESULTS
 
numTrials=length(Timeline);

    for i=1:numTrials;
        trial_alpha{i,1}=alpha_data{i,1};
        trial_alpha{i,3}=alpha_data{i,2};
        start_frame{i}=alpha_data{i,3};
        end_frame{i}=alpha_data{i,4};
        pks{i}=alpha_data{i,5};
        locs{i}=alpha_data{i,6};
    end

        
 %Write combined B+K plots for qualitative results
 AX1_top=round(max(Bioluminescence_quantitative.amplitude));
 AX2_top=round(max(Kinematics_quantitative.amplitude));

 figure(1); 
   for i=1:numTrials
        
        B=signal_data{i,2};
        K=alpha_data{i,2};
        
        %Calculate frame shift due to Hiris delayed recording trigger
        if Timeline.trial_stim_start(i)-Timeline.trial_frame_stim_start(i)>0;
        K_shift=Timeline.trial_stim_start(i)-Timeline.trial_frame_stim_start(i);
        K(K_shift:2500-1)=K(1:2500-K_shift);
        else K_shift=0;
        end
        
        %Subplot bioluminescence and alpha data on combined plots
        B_last=find(B,1,'last'); B(B_last:250)=B_last;
        K2=K(1:10:2500);
        T=1:1:250;
        if numTrials==30;
            figure(1),subplot_tight(1*numTrials/10,10,i);[AX,H1,H2]=plotyy(T(1:250),B(1:250),T(1:250),K2(1:250),'plot');title(['trial ',num2str(i)]);
        else if numTrials==10;
            figure(1),subplot_tight(1*numTrials/10*2,5,i);[AX,H1,H2]=plotyy(T(1:250),B(1:250),T(1:250),K2(1:250),'plot');title(['trial ',num2str(i)]);
            end
        end
        line([Timeline.trial_stim_start(i)/10,Timeline.trial_stim_start(i)/10],[-2,AX1_top],'color','r')
        set(H1,'color','g', 'LineWidth',1)
        set(H2,'color','b','LineWidth',1)

        %Set bioluminescence start/end and axis parameters
        hold(AX(1),'on')
        axes(AX(1))
        if Bioluminescence_quantitative.signal_start(i)~=0 && isnan(Bioluminescence_quantitative.signal_start(i))==0;
            plot(AX(1),Bioluminescence_quantitative.signal_start(i),signal_data{i,2}(Bioluminescence_quantitative.signal_start(i)),'og')
            plot(AX(1),Bioluminescence_quantitative.signal_end(i),signal_data{i,2}(Bioluminescence_quantitative.signal_end(i)),'og')
        end  
        set(AX(1),'XLim', [0 250])
        set(AX(1),'XTick',0:50:250)
        if AX1_top > 0; set(AX(1),'YLim', [-2 AX1_top]); set(AX(1),'YTick',-2:4:AX1_top);
        else set(AX(1),'YLim', [-2 10]); set(AX(2),'YTick',-2:4:10)
        end
        set(AX(1),'YColor','g')
        
        %Set alpha angle start/end/max and axis parameters
        hold(AX(2),'on')
        axes(AX(2))
        if start_frame{i,1}~=0 %%&& end_frame{i,1}~=0;
        plot(AX(2),(start_frame{i,1}+K_shift)/10,alpha_data{i,2}(start_frame{i,1}),'ob')
        plot(AX(2),(end_frame{i,1}+K_shift)/10,alpha_data{i,2}(end_frame{i,1}),'ob')
        plot(AX(2),(locs{i}+K_shift)/10,alpha_data{i,2}(locs{i}),'^k')
        end
        set(AX(2),'XLim', [0 250])
        if AX2_top > 0; set(AX(2),'YLim', [-AX2_top AX2_top]); set(AX(2),'YTick',-AX2_top:40:AX2_top)
        else set(AX(2),'YLim', [-180 180]); set(AX(2),'YTick',-180:40:180)
        end
        set(AX(2),'YColor','b')
        hold(AX(2),'off')
   end

    %Save results as matlab figure
    delete(fullfile(strrep(DataFile,'.hdf','_Qualitative.fig')));
    imageFileName=fullfile(strrep(DataFile,'.hdf','_Qualitative.fig'));
    saveas(figure(1),imageFileName,'fig');
    %%
    
end

clear all