%% Import data from excel
clear;
[~,~,PAPavTrainingData] = xlsread('PAPavTrainingDataALL.xlsx');

%% Assign data to corresponding variables within a structure
VarNames = PAPavTrainingData(1,:);
Data = PAPavTrainingData(2:end,:);

PAPav = struct();

for i=1:42
    PAPav.(VarNames{i}) = Data(1:end,(i));
end

PAPav.CSPlusOnset = Data(1:end,43:103);
PAPav.CSMinOnset = Data(1:end,104:163);
PAPav.PETimestamps = Data(1:end,225:942);

%% Convert rat numbers from strings into numbers and append to a separate field
for i = 1:44
   PAPav(:).RatNum(i,1) = str2num(PAPav.Rat{i}(3:4));
end

for i = 45:length(PAPav.Rat)
   PAPav(:).RatNum(i,1) = str2num(PAPav.Rat{i}(3:5));
end

%% Sort matrix by rat
[sortedRat,index] = sort([PAPav(:).Rat]);
fields = fieldnames(PAPav);
for i = 1:length(fields)
    PAPav.(fields{i}) = PAPav.(fields{i})(index,:);
end

%% Number the start dates for each rat and display any errors
Rat = PAPav.Rat(1);
k = 1;
for i=1:length(PAPav.StartDate)
   if  strcmp(PAPav.Rat(i), Rat)
       PAPav.Day(i,:) = k;
       k=k+1;
       if i < length(PAPav.StartDate)
           if PAPav.StartDate{i} == PAPav.StartDate{i+1};
               fprintf('repeated date %d with rat %s\n', PAPav.StartDate{i}, PAPav.Rat{i});
           end
       end
   else
       k = 1;
       Rat = PAPav.Rat(i);
       PAPav.Day(i,:) = k;
       k=k+1;
   end
end

%% Determine stages of training
% all stages of training are 5

for i=1:length(PAPav.Duration)
    PAPav.Stage{i,:} = 5;
end

%% Calculate trial-by-trial PE latencies for CSPlus and CSMin

% Convert necessary cells to matrices, setting all strings to NaN
for y=1:length(PAPav.CSPlusOnset)
    for x=1:size(PAPav.CSPlusOnset,2)
        if ischar(PAPav.CSPlusOnset{y,x})
           PAPav.CSPlusOnset{y,x} = NaN;
        end
    end
end

for y=1:length(PAPav.CSMinOnset)
    for x=1:size(PAPav.CSMinOnset,2)
        if ischar(PAPav.CSMinOnset{y,x})
           PAPav.CSMinOnset{y,x} = NaN;
        end
    end
end

for y=1:size(PAPav.PETimestamps,1)
    for x=1:size(PAPav.PETimestamps,2)
        if ischar(PAPav.PETimestamps{y,x})
           PAPav.PETimestamps{y,x} = NaN;
        end
    end
end

PAPav.CSPlusOnset = cell2mat(PAPav.CSPlusOnset);
PAPav.CSMinOnset = cell2mat(PAPav.CSMinOnset);
PAPav.PETimestamps = cell2mat(PAPav.PETimestamps);

% Calculate CSPlus absolute and cue period latencies
%   Conditions:
%       if next csplus > 0 but less than pe time, lat = nan
%       if next csmin > 0 but less than pe time, lat =nan

for y=1:length(PAPav.CSPlusOnset)
    for x=1:size(PAPav.CSPlusOnset,2)
       nextPlus = find(PAPav.CSPlusOnset(y,:) > PAPav.CSPlusOnset(y,x), 1);
       nextMin = find(PAPav.CSMinOnset(y,:) > PAPav.CSPlusOnset(y,x), 1);
       firstPE = find(PAPav.PETimestamps(y,:) > PAPav.CSPlusOnset(y,x), 1);
       PAPav.CSPlusAbsLat{y,x} = PAPav.PETimestamps(y,firstPE) - PAPav.CSPlusOnset(y,x);
       if PAPav.CSPlusAbsLat{y,x} < PAPav.Duration{y,1}
           PAPav.CSPlusCPLat{y,x} = PAPav.CSPlusAbsLat{y,x};
       else
           PAPav.CSPlusCPLat{y,x} = NaN;
       end
       if PAPav.CSMinOnset(y,nextMin) > 0 & PAPav.CSMinOnset(y,nextMin) < PAPav.PETimestamps(y,firstPE)
           PAPav.CSPlusAbsLat{y,x} = NaN;
           PAPav.CSPlusCPLat{y,x} = NaN;
       elseif PAPav.CSPlusOnset(y,nextPlus) > 0 & PAPav.CSPlusOnset(y,nextPlus) < PAPav.PETimestamps(y,firstPE)
           PAPav.CSPlusAbsLat{y,x} = NaN;
           PAPav.CSPlusCPLat{y,x} = NaN;
       elseif isnan(PAPav.CSPlusOnset(y,x))
           PAPav.CSPlusAbsLat{y,x} = NaN;
           PAPav.CSPlusCPLat{y,x} = NaN;
       elseif PAPav.CSPlusOnset(y,x) == 0
           PAPav.CSPlusAbsLat{y,x} = NaN;
           PAPav.CSPlusCPLat{y,x} = NaN;
       elseif isempty(firstPE)
           PAPav.CSPlusAbsLat{y,x} = NaN;
           PAPav.CSPlusCPLat{y,x} = NaN;
       end
    end
end

% Calculate CSPlus absolute and cue period latency means

PAPav.CSPlusAbsLat = cell2mat(PAPav.CSPlusAbsLat);
PAPav.CSPlusCPLat = cell2mat(PAPav.CSPlusCPLat);
PAPav.CSPlusAbsLatMean = nanmean(PAPav.CSPlusAbsLat,2);
PAPav.CSPlusCPLatMean = nanmean(PAPav.CSPlusCPLat,2);

% Calculate CSMin absolute and cue period latencies
%   Conditions:
%       if next csplus > 0 but less than pe time, lat = nan
%       if next csmin > 0 but less than pe time, lat =nan

for y=1:length(PAPav.CSMinOnset)
    for x=1:size(PAPav.CSMinOnset,2)
       nextPlus = find(PAPav.CSPlusOnset(y,:) > PAPav.CSMinOnset(y,x), 1);
       nextMin = find(PAPav.CSMinOnset(y,:) > PAPav.CSMinOnset(y,x), 1);
       firstPE = find(PAPav.PETimestamps(y,:) > PAPav.CSMinOnset(y,x), 1);
       PAPav.CSMinAbsLat{y,x} = PAPav.PETimestamps(y,firstPE) - PAPav.CSMinOnset(y,x);
       if PAPav.CSMinAbsLat{y,x} < PAPav.Duration{y,1}
           PAPav.CSMinCPLat{y,x} = PAPav.CSMinAbsLat{y,x};
       else
           PAPav.CSMinCPLat{y,x} = NaN;
       end
       if PAPav.CSMinOnset(y,nextMin) > 0 & PAPav.CSMinOnset(y,nextMin) < PAPav.PETimestamps(y,firstPE)
           PAPav.CSMinAbsLat{y,x} = NaN;
           PAPav.CSMinCPLat{y,x} = NaN;
       elseif PAPav.CSPlusOnset(y,nextPlus) > 0 & PAPav.CSPlusOnset(y,nextPlus) < PAPav.PETimestamps(y,firstPE)
           PAPav.CSMinAbsLat{y,x} = NaN;
           PAPav.CSMinCPLat{y,x} = NaN;
       elseif isnan(PAPav.CSMinOnset(y,x))
           PAPav.CSMinAbsLat{y,x} = NaN;
           PAPav.CSMinCPLat{y,x} = NaN;
       elseif PAPav.CSMinOnset(y,x) == 0
           PAPav.CSMinAbsLat{y,x} = NaN;
           PAPav.CSMinCPLat{y,x} = NaN;
       elseif isempty(firstPE)
           PAPav.CSMinAbsLat{y,x} = NaN;
           PAPav.CSMinCPLat{y,x} = NaN;
       end
    end
end

% Calculate CSMin absolute and cue period latency means

PAPav.CSMinAbsLat = cell2mat(PAPav.CSMinAbsLat);
PAPav.CSMinCPLat = cell2mat(PAPav.CSMinCPLat);
PAPav.CSMinAbsLatMean = nanmean(PAPav.CSMinAbsLat,2);
PAPav.CSMinCPLatMean = nanmean(PAPav.CSMinCPLat,2);

%% Determine when rats reach criteria
% FOR EACH RAT:
% criteria = 0 if cs+ratio < .7 or cs-ratio > .3
% if cs+ratio > .7 and cs-ratio < .3, then
%   for the first time this happens, criteria = 0
%   for every other time this happens and even if they don't meet criteria in later dates (don't need to check cs ratios), criteria = 1
% if new rat, start over

for i=1:length(PAPav.CSPlusRatio)
    if ischar(PAPav.CSPlusRatio{i})
        fprintf('no cs+ ratio on %d with rat %s\n', PAPav.StartDate{i}, PAPav.Rat{i})
        PAPav.CSPlusRatio{i} = NaN;
    end
end

for i=1:length(PAPav.CSMinRatio)
    if ischar(PAPav.CSMinRatio{i})
        fprintf('no cs+ ratio on %d with rat %s\n', PAPav.StartDate{i}, PAPav.Rat{i})
        PAPav.CSMinRatio{i} = NaN;
    end
end

criteria_met = false;
Rat = PAPav.Rat(1);
k = 1;

for i=1:length(PAPav.Stage)
    if strcmp(PAPav.Rat(i), Rat)
        if criteria_met == true % this condition is set from criteria being met
            PAPav.Criteria{i,1} = 1;
        elseif PAPav.CSPlusRatio{i} >= 0.7 & PAPav.CSMinRatio{i} <= 0.3 % correct criteria
            PAPav.Criteria{i,1} = 0;
            PAPav.CriteriaDay(k,1) = PAPav.Day(i,1); 
            PAPav.CriteriaDayRat(k,1) = PAPav.RatNum(i,1);
            PAPav.CriteriaDayCSPlusRatio(k,1) = PAPav.CSPlusRatio(i,1);
            PAPav.CriteriaDayCSMinRatio(k,1) = PAPav.CSMinRatio(i,1);
            PAPav.CriteriaDayCSPlusLat(k,1) = PAPav.CSPlusCPLatMean(i,1);
            PAPav.CriteriaDayCSMinLat(k,1) = PAPav.CSMinCPLatMean(i,1);
            k = k+1;
            criteria_met = true; % false switched to true
        else
            PAPav.Criteria{i,1} = 0;
        end
    else
       Rat = PAPav.Rat(i);
       PAPav.Criteria{i,1} = 0;
       criteria_met = false;
    end
end

for i=1:length(PAPav.Criteria)
    if isnan(PAPav.CSPlusRatio{i,1})
        PAPav.Criteria{i,1} = NaN;
    end
end

for i=1:length(PAPav.Criteria)
    if isnan(PAPav.CSMinRatio{i,1})
        PAPav.Criteria{i,1} = NaN;
    end
end

PAPav.Criteria = cell2mat(PAPav.Criteria);
PAPav.CSPlusRatio = cell2mat(PAPav.CSPlusRatio);
PAPav.CSMinRatio = cell2mat(PAPav.CSMinRatio);

%% Calculate averages of latency means on days where criteria = 0 during cue period (group data)

% To find averages of CS+ and CS- Cue Period Latency Means by day

RatSelection=ismember(PAPav.RatNum,[36 122 137 140 154]);


for i=1:max(PAPav.Day(PAPav.Criteria == 0))
    PAPav.AverageCSPlusCPLatMean(i) = nanmean(PAPav.CSPlusCPLatMean(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection));
    PAPav.Error_AverageCSPlusCPLatMean(i) = nanste(PAPav.CSPlusCPLatMean(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection), 1);
end

for i=1:max(PAPav.Day(PAPav.Criteria == 0))
    PAPav.AverageCSMinCPLatMean(i) = nanmean(PAPav.CSMinCPLatMean(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection));
    PAPav.Error_AverageCSMinCPLatMean(i) = nanste(PAPav.CSMinCPLatMean(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection), 1);
end


% To find averages of CS+ and CS- Ratios by day
for i=1:max(PAPav.Day(PAPav.Criteria == 0))
    PAPav.AverageCSPlusRatio(i) = nanmean(PAPav.CSPlusRatio(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection));
    PAPav.Error_AverageCSPlusRatio(i) = nanste(PAPav.CSPlusRatio(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection), 1);
end

for i=1:max(PAPav.Day(PAPav.Criteria == 0))
    PAPav.AverageCSMinRatio(i) = nanmean(PAPav.CSMinRatio(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection));
    PAPav.Error_AverageCSMinRatio(i) = nanste(PAPav.CSMinRatio(PAPav.Criteria == 0 & PAPav.Day == i & RatSelection), 1);
end



%% Save file
filename = 'PavTrainingAnalysis.mat';
save(filename);