% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code has been written by Adelie Antoine, using functions from
% animaltags (M. Johnson, Animal Tag Tools Wiki (2017).
% URL https:// animaltags.org/), and following P. Goulet's work.
%
% It is designed to extract behavioral characteristics within dive clusters
%
% Last modification : February 2023
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all
close all

%% Load the dataset containing the results of the echograms visual analysis

% The structure must contain (with array names in quotes):
% - individual identification, "indiv"
% - yes/no answers for the presence of a signal corresponding to the bottom
% of the ocean on the echogram, "presence_bottom"
% - yes/no answers for the presence of a signal corresponding to a prey on 
% the echogram, "presence_prey"
% - yes/no answers for the presence of individual prey (as opposed to
% several preys), "indiv_prey"
% - yes/no answers for the presence of a signal corresponding to a reacting
% prey, "prey_reaction"
% - type of reaction, "reaction_type"
% - time at which the reaction is first spotted, "reaction_time"
% - distance at which the reaction is first spotted, "reaction_distance"
% - prey acoustic size, "prey_size"
% - index of the PrCA corresponding to the echogram that was analyzed,
% "index_jerk"

% The structure is hereafter called "sonar_game"

% Prey capture attemptes are hereafter noted (PCA).

% We suppose that the studied individual is identified as "314a"


%% load all necessary data
classif = readmatrix('.../lw19_314a_cluster.csv'); % load the clustering results 
load('.../lw19_314a_dives.mat') % load file containing dives for the studied individual

load('.../lw19_314a_jerks_fullbandwidth.mat') % load file containing jerk data at full bandwidth

addpath(genpath('.../dtag_and_sonar_tools/'));
load_nc('.../lw19_314a_sensor.nc'); % load file containing tag data at 5Hz


%% Extract PCA index + corresponding dive
% Echograms were analyzed in random order,
% get the PCA index in the random list of analyzed echograms for the individual to study
idx_pca_14a_random_file = find(sonar_game.indiv=="lw19_314a");

% get the PCA index in the original file (fullBW file)
idx_pca_14a = sonar_game.index_jerk(idx_pca_14a_random_file)';
clear idx_pca_14a_random_file

mean(jerks.depth(idx_pca_14a))
std (jerks.depth(idx_pca_14a))

% from now on, use idx_pca_14a in jerk file to index PCAs only among jerks

% get corresponding dive index
% get start and end time of jerks 
jerk_start = jerks.start_time(idx_pca_14a); % extract pca start time
jerk_end = jerks.end_time(idx_pca_14a); % extract pca end time
idx_dives = []; % initate dive index double
for ii = 1 : length(jerk_start) % for all PCA
    idx = find(dives.start <= jerk_start(ii,1) & dives.end >= jerk_end(ii,1)); % get corresponding dive index 
    idx_dives =cat(1,idx_dives,idx); % save index of all dives with PCA in a double, /!\ including repeated dives if several PCAs
end 

[GC,GR] = groupcounts(idx_dives); % GC: counts, GR: element that were counted (equivalent to unique(idx_dives))
mean(GC)
std(GC)


idx_dives_pca_14a = unique(idx_dives); % here are corresponding dive index
clear idx


% TO SELECT DIVES THAT ARE BOTH CLASSIFIED AND INCLUDE PCA

% extract corresponding index of classified dives with PCA 
a1 = intersect(classif(:,2), idx_dives_pca_14a(:,1)); % index of dives with PCA and that are classified, for all clusters 
classif_314a_pca_idx = [];
for ik = 1 : length (a1)
    idx = a1(ik,1); % dive index among all dives
    idx_clust = find(classif(:,2) == idx);  % index of classified dives with PCA in classif_314a file
    classif_314a_pca_idx = cat(1, classif_314a_pca_idx, idx_clust); % index of classified dives with PCA in classif_314a file
end 

% create file to replace classif_314a, containing cluster number for dives
% with PCA
classif_314a_pca = [];
for im = 1 : length (classif_314a_pca_idx)
    idx = classif_314a_pca_idx(im, 1); % index of classified dives with PCA in classif_314a file
    classif_314a_pca(im,1) = classif(idx,1);
    classif_314a_pca(im,2) = classif(idx,2); % new double that replaces classif_314a, containing cluster number only for dives with PCA
end

% idx of classified dives with prca
c1_pca = classif_314a_pca((find(classif_314a_pca(:,1)==1)),2);
c2_pca = classif_314a_pca((find(classif_314a_pca(:,1)==2)),2);
c3_pca = classif_314a_pca((find(classif_314a_pca(:,1)==3)),2);
c4_pca = classif_314a_pca((find(classif_314a_pca(:,1)==4)),2);
c5_pca = classif_314a_pca((find(classif_314a_pca(:,1)==5)),2);

%idx of non classified dives with prca
% there are doublons in idx_dives file, remove them
unique_idx_dives = unique(idx_dives);


% loop to save boxplot and violin plot 
for ij = 1 : 5

%% %%%%%%%%%%%%%%%%%%%%% PART 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Orientation: euler angles (pitch, roll, heading) %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 1 - Select jerks for classified dives 
    
classif_clust1 = find(classif_314a_pca(:,1)==ij); % extract row number of cluster ij dives with PCAs in file
clust1 = classif_314a_pca(classif_clust1,2); % index of dives in cluster ij that also have PCAS

dive_start_c1 = dives.start(clust1); % extract classified dive with PCA start time
dive_end_c1 = dives.end(clust1); % extract classified dive with PCA end time

classified_jerks_idx_c1 = []; % initate jerks index double
for ii = 1 : length(dive_start_c1) % for all classified dives in cluster 1
    classifjerk_c1 = find(jerks.start_time(idx_pca_14a)>=dive_start_c1(ii,1) & jerks.end_time(idx_pca_14a)<=dive_end_c1(ii,1)); % get corresponding jerk index (ie jerks for classified dives)
    classified_jerks_idx_c1 =cat(1,classified_jerks_idx_c1,classifjerk_c1); % save all classified jerks index in a double 
end 

classified_jerks_c1 = ([]); % create empty structure 
classified_jerks_c1.start_time = jerks.start_time(classified_jerks_idx_c1); % only need jerk start time for rest of script 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 2 - Compute orientation

% Apply a correction factor to the accelerometer and magnetometer axes so that when the seal is horizontal, its pitch is zero
v = depth_rate(P); % compute the vertical speed
vmax = quantile(v,0.9);
[pe,s] = infer_pitch_error(P,A,vmax);
corr_pitch = 5*ceil((pe*180/pi)/5); %round to the closest multiple of 5, pitch offset in degree

pitch_err = corr_pitch*pi/180 ;	% pitch offset of tag below horizontal in radians

A.data = A.data*euler2rotmat([pitch_err,0,0]) ;	% correct A for tag orientation
M.data = M.data*euler2rotmat([pitch_err,0,0]) ;	% correct M for tag orientation
A.rotation_matrix = [pitch_err,0,0];
M.rotation_matrix = [pitch_err,0,0];

% Compute orientation at 1Hz
[pitch,roll] = a2pr(decdc(A,A.samplingrate)) ;	% pitch and roll angles, directly with a 1Hz sampling frequency
%[pitch,roll] = a2pr(A) ; % here at 5Hz
heading =  m2h(decdc(M,M.samplingrate),decdc(A,A.samplingrate));	% heading angle, directly with a 1Hz sampling frequency	
%heading =  m2h(M,A); % here at 5Hz

% Collect Euler angles for each PrCA
fs = 1; % sampling frequency (1 Hz because we computed pitch and roll at 1 Hz)
lag = 10; % time before the PrCA

pitch_capt = NaN(length(classified_jerks_c1.start_time),11); % 1 line per PrCA, 1 column per second preceeding the PrCA
roll_capt = NaN(length(classified_jerks_c1.start_time),11);
heading_capt = NaN(length(classified_jerks_c1.start_time),11);
for i = 1:length(classified_jerks_c1.start_time)
    start_time = round(classified_jerks_c1.start_time(i))*fs-(lag*fs); % start_time/end time is in seconds, pitch (or A, P, J...) is in sample. 1 seconds = 1 * fs samples. 
    end_time = round(classified_jerks_c1.start_time(i))*fs;   
    pitch_capt(i,:) = pitch(start_time:end_time) *180/pi; % *180/pi: to have the angle in degree
    roll_capt(i,:) = roll(start_time:end_time)' *180/pi;
    heading_capt(i,:) = heading(start_time:end_time)' *180/pi;
end

% % if you want to visualise the pitch/roll during the approach phase:
figure
subplot(3,1,1)
boxplot(pitch_capt)
yline(0);
title ('cluster nb', ij)
ylim([-180 180]); ylabel('Pitch angle (°)')
subplot(3,1,2)
boxplot(roll_capt); ylabel('Roll angle (°)')
yline(0);
ylim([-180 180])
subplot(3, 1, 3)
boxplot(heading_capt); ylabel('Heading (°)')
yline(0);
filename = sprintf('.../314a_boxplot_cl_%d_10s', ij);
set(gcf, 'PaperUnits', 'inches', 'PaperPosition', [0 0 20 15]) % image size 2000 * 1500 pixels 
print('-dpng', filename);   % for PNG file (lossless)

% % violin plot
figure
subplot(3,1,1)
vs = violinplot(pitch_capt);
ylabel('Pitch angle (°)');
title ('cluster nb', ij)
subplot(3,1,2)
vs = violinplot(roll_capt);
ylabel('Roll (°)');
subplot(3,1,3)
vs = violinplot(heading_capt);
ylabel('Heading angle (°)');
filename = sprintf('.../314a_violin_cl_%d_10s', ij);
set(gcf, 'PaperUnits', 'inches', 'PaperPosition', [0 0 20 15]) % image size 2000 * 1500 pixels 
print('-dpng', filename);   % for PNG file (lossless)

close all
end 


%% %%%%%%%%%%%%%%%%%%%%% PART 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute number of flipper strokes during the approach phase for each PrCA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Apply a high pass filter to obtain high frequency acceleration data
fc = 0.6; % cutting-off frequency: used to separate high frequency movements due to flipper strokes from other movements due to orientation of the animal
Af = comp_filt(A,fc) ;
Ahigh = Af{2};	% high frequency A data
fs = A.samplingrate; % sampling frequency 
lag = 10; %time before the PrCA

% Compute nb of flipper strokes during the approach phase
nb_strokes = NaN(length(classified_jerks_c1.start_time),1); % 1 value per PrCA
for i = 1:length(classified_jerks_c1.start_time)
    start_time = round(classified_jerks_c1.start_time(i))*fs-(lag*fs);
    end_time = round(classified_jerks_c1.start_time(i))*fs;
    acc_capture = Ahigh(start_time:end_time,2);
    [K,s,KK] = zero_crossings(acc_capture,0.5);
    nb_strokes(i) = round(length(K)/2); 
end


number_strokes(ij,1) = mean(nb_strokes)


%% %%%%%%%%%%%%%%%%%%%%% PART 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% look at flipper strokes and lat acc for each PrCA %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v = depth_rate(P); % compute the vertical speed
vmax = quantile(v,0.9);
[pe,s] = infer_pitch_error(P,A,vmax);
corr_pitch = 5*ceil((pe*180/pi)/5); %round to the closest multiple of 5, pitch offset in degree
pitch_err = corr_pitch*pi/180 ;	% pitch offset of tag below horizontal in radians
A.data = A.data*euler2rotmat([pitch_err,0,0]) ;	% correct A for tag orientation
M.data = M.data*euler2rotmat([pitch_err,0,0]) ;	% correct M for tag orientation
A.rotation_matrix = [pitch_err,0,0];
M.rotation_matrix = [pitch_err,0,0];
% Compute orientation at 1Hz
[pitch,roll] = a2pr(decdc(A,A.samplingrate)) ;			% pitch and roll angles, directly with a 1Hz sampling frequency
heading =  m2h(decdc(M,M.samplingrate),decdc(A,A.samplingrate));	% heading angle, directly with a 1Hz sampling frequency	
prh  = [pitch roll heading];
% get lateral acceleration to visualise the flipper strokes
fc = 0.6;
Af = comp_filt(A,fc) ;
Ahigh = Af{2};	


% select PCA for classified dives 
classified_jerks_idx = []; % create empty structure 
for ij = 1:5 % for each cluster, to save all PCAs needed 
classif_clust1 = find(classif_314a_pca(:,1)==ij); % extract row number of cluster ij dives with PCAs in file
clust1 = classif_314a_pca(classif_clust1,2); % index of dives in cluster ij that also have PCAS

dive_start_c1 = dives.start(clust1); % extract classified dive with PCA start time
dive_end_c1 = dives.end(clust1); % extract classified dive with PCA end time

classified_jerks_idx_c1 = []; % initate jerks index double
for ii = 1 : length(dive_start_c1) % for all classified dives in cluster 1
    classifjerk_c1 = find(jerks.start_time(idx_pca_14a)>=dive_start_c1(ii,1) & jerks.end_time(idx_pca_14a)<=dive_end_c1(ii,1)); % get corresponding jerk index (ie jerks for classified dives)
    classified_jerks_idx_c1 =cat(1,classified_jerks_idx_c1,classifjerk_c1); % save all classified jerks index in a double 
end 

classified_jerks_idx = cat(1,classified_jerks_idx,classified_jerks_idx_c1); % create empty structure 
classified_jerks.start_time = jerks.start_time(classified_jerks_idx); % only need jerk start time for rest of script 
end 


% Plotting:
for i = 1:length(classified_jerks_idx) % for each pca 
for i = 182 % for a single pca

% crop the orientation + acceleration data at the time of the PrCA
    idx = classified_jerks_idx(i,1); % put the jerk index you are interested in
    start_time = round(jerks.start_time(idx)); 
    end_time = round(jerks.start_time(idx)); 
    lag = 10; % number of seconds before the PrCA
    fs_prh = 1; % sampling rate for the pitch and the roll
    fs_swim = 5; % sampling rate for acceleration data (used to 
    prh_prca = prh(start_time-lag*fs_prh:end_time,:); % crop body angles at the time of the PrCA
    lat_acc_prca = Ahigh(start_time*fs_swim-(fs_swim*lag):end_time*fs_swim,2); % crop the lateral acceleration of the filtered data at the time of the PrCA

% for acceleration plots
    start_time_acc = classified_jerks.start_time(i,1); % start time of pca, accounts for A.data index too
    acc_three_axes(:,1) = A.data (start_time_acc * fs_swim-(fs_swim*lag):start_time_acc* fs_swim,1); % crop acceleration data in 1st axis 10 sec before start 
    acc_three_axes(:,2) = A.data (start_time_acc * fs_swim-(fs_swim*lag):start_time_acc* fs_swim,2); % crop acceleration data in 2nd axis 10 sec before start 
    acc_three_axes(:,3) = A.data(start_time_acc * fs_swim-(fs_swim*lag):start_time_acc* fs_swim,3); % crop acceleration data in 3rd axis 10 sec before start 

    figure; 
ax(1) = subplot(5,1,1);
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,1)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Pitch')
   title('pca nb ',idx)

ax(2) = subplot(5,1,2);
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,2)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Roll')

ax(3) = subplot(5,1,3); % ajout heading ADELIE 11/02/22
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,3)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Heading')

ax(4) = subplot(5,1,4); 
plot( (1:length(lat_acc_prca))/fs_swim,    lat_acc_prca, 'k');
   grid;
%    linkaxes(ax, 'x');
   xlim([0 length(lat_acc_prca)/fs_swim]);
   ylabel('Lateral acceleration')

ax(5) = subplot(5, 1, 5);
plot ((1:length(acc_three_axes))/fs_swim, acc_three_axes); % plot acceleration in 3 axes
    grid;
   %xlim([0 length(prh_prca)/fs_prh]);
   xlim([0 length(acc_three_axes)/fs_swim]);
   ylabel('Acceleration in three axes')


   filename = sprintf('.../approach_pca_nb_%d', idx);
set(gcf, 'PaperUnits', 'inches', 'PaperPosition', [0 0 20 15]) % image size 2000 * 1500 pixels 
print('-dpng', filename);   % for PNG file (lossless)

close all 
end 




%% %%%%%%%%%%%%%%%%%%%%% PART 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% look at lateral acceleration for each PrCA with reacting prey %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find dives corresponding to PCA with reacting prey  

idx_reaction_all_types_yes = sonar_game.index_jerk(sonar_game.indiv=="lw19_314a" & sonar_game.prey_reaction == 1);
idx_reaction_all_types_dont_know = sonar_game.index_jerk(sonar_game.indiv=="lw19_314a" & sonar_game.prey_reaction == -1);
idx_reaction_all_types = cat(1,idx_reaction_all_types_yes,idx_reaction_all_types_dont_know);

% list all indices of PCA per prey reaction
idx_pca_escaping = sonar_game.index_jerk(sonar_game.indiv=="lw19_314a" & sonar_game.reaction_type == 2);
idx_pca_leaving = sonar_game.index_jerk(sonar_game.indiv=="lw19_314a" & sonar_game.reaction_type == 1);


% get dive index corresponding to each reaction type
idx_dives_escaping = []; 
for ii = 1 : length(idx_pca_escaping) 
    k = idx_pca_escaping(ii,1); % pca index
    jerks_start = jerks.start_time(k); % corresponding pca start time
    jerks_end = jerks.end_time(k); % corresponding pca end time
    idx = find (dives.start<=jerks_start & dives.end>=jerks_end); % find corresponding dive
    idx_dives_escaping = cat(1,idx_dives_escaping,idx); % store results, in same order as pca index
end

idx_dives_leaving = []; 
for ii = 1 : length(idx_pca_leaving) 
    k = idx_pca_leaving(ii,1); % pca index
    jerks_start = jerks.start_time(k); % corresponding pca start time
    jerks_end = jerks.end_time(k); % corresponding pca end time
    idx = find (dives.start<=jerks_start & dives.end>=jerks_end); % find corresponding dive
    idx_dives_leaving = cat(1,idx_dives_leaving,idx); % store results, in same order as pca index
end

% plot acceleration for each pca
v = depth_rate(P); % compute the vertical speed
vmax = quantile(v,0.9);
[pe,s] = infer_pitch_error(P,A,vmax);
corr_pitch = 5*ceil((pe*180/pi)/5); %round to the closest multiple of 5, pitch offset in degree
pitch_err = corr_pitch*pi/180 ;	% pitch offset of tag below horizontal in radians
A.data = A.data*euler2rotmat([pitch_err,0,0]) ;	% correct acceleration for tag orientation
M.data = M.data*euler2rotmat([pitch_err,0,0]) ;	% correct magnetometry for tag orientation
A.rotation_matrix = [pitch_err,0,0];
M.rotation_matrix = [pitch_err,0,0];
% Compute orientation at 1Hz
[pitch,roll] = a2pr(decdc(A,A.samplingrate)) ;	% pitch and roll angles, directly with a 1Hz sampling frequency
heading =  m2h(decdc(M,M.samplingrate),decdc(A,A.samplingrate)); % heading angle, directly with a 1Hz sampling frequency	
prh  = [pitch roll heading];
% get lateral acceleration to visualise the flipper strokes
fc = 0.6;
Af = comp_filt(A,fc) ;
Ahigh = Af{2};	

pca_escaping.start_time = jerks.start_time(idx_pca_escaping); 
pca_escaping.end_time = jerks.end_time(idx_pca_escaping); 
pca_leaving.start_time = jerks.start_time(idx_pca_leaving); 
pca_leaving.end_time = jerks.end_time(idx_pca_leaving); 

% plot acceleration and orientation for the whole pca 
%% escaping preys
for i = 1:length(idx_pca_escaping) % for each pca 

% crop the orientation + acceleration data at the time of the PrCA
    idx = idx_pca_escaping(i,1); % put the pca index you are interested in
    start_time = round(jerks.start_time(idx)); 
    end_time = round(jerks.end_time(idx));  % /!\ here we look at the WHOLE pca !!!!!!!
    lag = 10; % number of seconds before the PrCA
    lag_end = 5; % number of seconds after the end of PrCA
    fs_prh = 1; % sampling rate for the pitch and the roll
    fs_swim = 5; % sampling rate for acceleration data 
    prh_prca = prh(start_time-lag*fs_prh:end_time+lag_end*fs_prh,:); % crop body angles at the time of the PrCA
    lat_acc_prca = Ahigh(start_time*fs_swim-(fs_swim*lag):end_time*fs_swim+(fs_swim*lag_end),2); % crop the lateral acceleration of the filtered data at the time of the PrCA

% for acceleration plots
    acc_three_axes = [];
    start_time_acc = pca_escaping.start_time(i,1); % start time of pca, accounts for A.data index too
    end_time_acc = pca_escaping.end_time(i,1); % end time of pca, accounts for A.data index too
    acc_three_axes(:,1) = A.data (round(start_time_acc * fs_swim-(fs_swim*lag):end_time_acc* fs_swim+(fs_swim*lag_end)),1); % crop acceleration data in 1st axis 10 sec before start 
    acc_three_axes(:,2) = A.data (round(start_time_acc * fs_swim-(fs_swim*lag):end_time_acc* fs_swim+(fs_swim*lag_end)),2); % crop acceleration data in 2nd axis 10 sec before start 
    acc_three_axes(:,3) = A.data (round(start_time_acc * fs_swim-(fs_swim*lag):end_time_acc* fs_swim+(fs_swim*lag_end)),3); % crop acceleration data in 3rd axis 10 sec before start 

% to get time when prey escapes
    escape_time = sonar_game.reaction_time(sonar_game.index_jerk == idx & sonar_game.indiv=="lw19_314a") + (lag-5);
% lag-5 because lag = 10s but sonar game displayed 5s before pca 
% therefore recorded reaction time is 5s later than actual start of pca

    figure; 
ax(1) = subplot(5,1,1);
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,1)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Pitch')
   title('pca nb ',idx)
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r'); % line at escaping time
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); %line at start of pca
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(2) = subplot(5,1,2);
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,2)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Roll')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r');
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(3) = subplot(5,1,3); % ajout heading ADELIE 11/02/22
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,3)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Heading')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r');
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(4) = subplot(5,1,4); 
plot( (1:length(lat_acc_prca))/fs_swim,    lat_acc_prca, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(lat_acc_prca)/fs_swim]);
   ylabel('Lateral acceleration')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r');
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(5) = subplot(5, 1, 5);
plot ((1:length(acc_three_axes))/fs_swim, acc_three_axes); % plot acceleration in 3 axes
    grid;
   linkaxes(ax, 'x');
   xlim([0 length(acc_three_axes)/fs_swim]);
   ylabel('Acceleration in three axes')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r'); 
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca
   legend({'start/end of pca','prey escapes'},'Location','northwest')

   filename = sprintf('.../pca_nb_%d', idx);
set(gcf, 'PaperUnits', 'inches', 'PaperPosition', [0 0 20 15]) % image size 2000 * 1500 pixels 
print('-dpng', filename);   % for PNG file (lossless)

close all 
end 

%% leaving preys
for i = 1:length(idx_pca_leaving) % for each pca 
for i = 182
% crop the orientation + acceleration data at the time of the PrCA
    idx = idx_pca_leaving(i,1); % put the pca index you are interested in
    start_time = round(jerks.start_time(idx)); 
    end_time = round(jerks.end_time(idx));  % /!\ here we look at the WHOLE pca !!!!!!!
    lag = 10; % number of seconds before the PrCA
    lag_end = 5; % number of seconds after the end of PrCA
    fs_prh = 1; % sampling rate for the pitch and the roll
    fs_swim = 5; % sampling rate for acceleration data 
    prh_prca = prh(start_time-lag*fs_prh:end_time+lag_end*fs_prh,:); % crop body angles at the time of the PrCA
    lat_acc_prca = Ahigh(start_time*fs_swim-(fs_swim*lag):end_time*fs_swim+(fs_swim*lag_end),2); % crop the lateral acceleration of the filtered data at the time of the PrCA

% for acceleration plots
    acc_three_axes = [];
    start_time_acc = pca_leaving.start_time(i,1); % start time of pca, accounts for A.data index too
    end_time_acc = pca_leaving.end_time(i,1); % end time of pca, accounts for A.data index too
    acc_three_axes(:,1) = A.data (round(start_time_acc * fs_swim-(fs_swim*lag):end_time_acc* fs_swim+(fs_swim*lag_end)),1); % crop acceleration data in 1st axis 10 sec before start 
    acc_three_axes(:,2) = A.data (round(start_time_acc * fs_swim-(fs_swim*lag):end_time_acc* fs_swim+(fs_swim*lag_end)),2); % crop acceleration data in 2nd axis 10 sec before start 
    acc_three_axes(:,3) = A.data (round(start_time_acc * fs_swim-(fs_swim*lag):end_time_acc* fs_swim+(fs_swim*lag_end)),3); % crop acceleration data in 3rd axis 10 sec before start 

% to get time when prey escapes
    escape_time = sonar_game.reaction_time(sonar_game.index_jerk == idx & sonar_game.indiv=="lw19_314a") + (lag-5);


    figure; 
ax(1) = subplot(5,1,1);
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,1)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Pitch')
   title('pca nb ',idx)
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r'); % line at escaping time
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); %line at start of pca
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(2) = subplot(5,1,2);
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,2)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Roll')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r');
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(3) = subplot(5,1,3); % ajout heading ADELIE 11/02/22
plot( (1:length(prh_prca))/fs_prh,    prh_prca(:,3)*180/pi, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(prh_prca)/fs_prh]);
   ylim([-180 180])
   ylabel('Heading')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r');
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(4) = subplot(5,1,4); 
plot( (1:length(lat_acc_prca))/fs_swim,    lat_acc_prca, 'k');
   grid;
   linkaxes(ax, 'x');
   xlim([0 length(lat_acc_prca)/fs_swim]);
   ylabel('Lateral acceleration')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r');
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca

ax(5) = subplot(5, 1, 5);
plot ((1:length(acc_three_axes))/fs_swim, acc_three_axes); % plot acceleration in 3 axes
    grid;
   linkaxes(ax, 'x');
   xlim([0 length(acc_three_axes)/fs_swim]);
   ylabel('Acceleration in three axes')
   line([escape_time, escape_time], ylim, 'LineWidth', 2, 'Color', 'r'); 
   line([10, 10], ylim, 'LineWidth', 2, 'Color', 'b'); 
   line ([ (length(prh_prca))-5, (length(prh_prca))-5], ylim, 'LineWidth', 2, 'Color', 'b'); % line at end of pca
   legend({'start/end of pca','prey leaves'},'Location','northwest')

   filename = sprintf('.../pca_nb_%d', idx);
set(gcf, 'PaperUnits', 'inches', 'PaperPosition', [0 0 20 15]) % image size 2000 * 1500 pixels 
print('-dpng', filename);   % for PNG file (lossless)

%close all 
end 


% check which dive index are the PCA in 
classif = readmatrix('.../lw19_314a_cluster.csv');

% escaping preys
classif_dives_escaping = [];
for i = 1 : length (idx_dives_escaping)
    k = idx_dives_escaping(i,1);
    idx = find(classif(:,2)==k);
    classif = classif(idx,1);
    classif_dives_escaping = cat(1, classif_dives_escaping, classif);
end 

escape_cl_1 = sum (classif_dives_escaping(:,1) == 1)/1;
escape_cl_2 = sum (classif_dives_escaping(:,1) == 2)/2;


% leaving preys
classif_dives_leaving = [];
for i = 1 : length (idx_dives_leaving)
    k = idx_dives_leaving(i,1);
    idx = find(classif(:,2)==k);
    classif = classif(idx,1);
    classif_dives_leaving = cat(1, classif_dives_leaving, classif); 
end 

% to compare reacting prey size / seal's response and reacting prey
% distance / seal's response,
% look at mean + std of prey characteristic and compare in excel file with seal's reaction

% prey size 
idx_prey_size_escaping = find(sonar_game.indiv=="lw19_314a" & sonar_game.reaction_type ==2);
prey_size_escaping = sonar_game.prey_size(idx_prey_size_escaping);
nb_nan_size_escaping = sum(isnan(prey_size_escaping));
mean_size_escaping = mean(prey_size_escaping, 'omitnan');
std_size_escaping = std(prey_size_escaping, 'omitnan');
max_size_escaping = max(prey_size_escaping);
min_size_escaping = min(prey_size_escaping);

idx_prey_size_leaving = find(sonar_game.indiv=="lw19_314a" & sonar_game.reaction_type ==1);
prey_size_leaving = sonar_game.prey_size(idx_prey_size_leaving);
nb_nan_size_leaving = sum(isnan(prey_size_leaving));
mean_size_leaving = mean(prey_size_leaving, 'omitnan');
std_size_leaving = std(prey_size_leaving, 'omitnan');
max_size_leaving = max(prey_size_leaving);
min_size_leaving = min(prey_size_leaving);

prey_size_escaping(:,2) = sonar_game.index_jerk(idx_prey_size_escaping);
prey_size_leaving(:,2) = sonar_game.index_jerk(idx_prey_size_leaving);

% prey distance
idx_prey_distance_escaping = find(sonar_game.indiv=="lw19_314a" & sonar_game.reaction_type ==2);
prey_distance_escaping = sonar_game.reaction_distance(idx_prey_distance_escaping);
nb_nan_distance_escaping = sum(isnan(prey_distance_escaping));
mean_distance_escaping = mean(prey_distance_escaping, 'omitnan');
std_distance_escaping = std(prey_distance_escaping, 'omitnan');
max_distance_escaping = max(prey_distance_escaping);
min_distance_escaping = min(prey_distance_escaping);

idx_prey_distance_leaving = find(sonar_game.indiv=="lw19_314a" & sonar_game.reaction_type ==1);
prey_distance_leaving = sonar_game.reaction_distance(idx_prey_distance_leaving);
nb_nan_distance_leaving = sum(isnan(prey_distance_leaving));
mean_distance_leaving = mean(prey_distance_leaving, 'omitnan');
std_distance_leaving = std(prey_distance_leaving, 'omitnan');
max_distance_leaving = max(prey_distance_leaving);
min_distance_leaving = min(prey_distance_leaving);

prey_distance_escaping(:,2) = sonar_game.index_jerk(idx_prey_distance_escaping);
prey_distance_leaving(:,2) = sonar_game.index_jerk(idx_prey_distance_leaving);


% to compare prey size / prey reaction and prey distance / prey reaction
% create a csv file with distance, size, reaction type 

size_reaction = cat(1, prey_size_escaping(:,1), prey_size_leaving(:,1));
distance_reaction = cat (1, prey_distance_escaping, prey_distance_leaving);
size_distance_prey_reaction = cat(2, size_reaction, distance_reaction);
size_distance_prey_reaction = array2table(size_distance_prey_reaction,'variableNames',{'size','distance','pcanb'});
size_distance_prey_reaction.reaction(1:20) = ["escape"];
size_distance_prey_reaction.reaction(21:65) = ["leave"];

% filename  =['.../size_distance_prey_reaction_314a.csv'];
% writetable(size_distance_prey_reaction, filename);


%% find cluster for each PCA with reacting prey 
idx_pca_14a = size_distance_prey_reaction.pcanb;

% from now on, use idx_pca_14a in jerk file to index PCAs only among jerks

% get corresponding dive index
% get start and end time of jerks 
jerk_start = jerks.start_time(idx_pca_14a); % extract pca start time
jerk_end = jerks.end_time(idx_pca_14a); % extract pca end time
idx_dives = []; % initate dive index double
for ii = 1 : length(idx_pca_14a) % for all PCA
    k = idx_pca_14a(ii,1);
    idx = find(dives.start <= jerk_start(k) & dives.end >= jerk_end(k)); % get corresponding dive index 
    a = find(classif(:,2)==idx_pca_14a(ii,1));
    idx_dives =cat(1,idx_dives,idx); % save index of all dives with PCA in a double, /!\ including repeated dives if several PCAs
end 
idx_dives_pca_14a = unique(idx_dives); % here are corresponding dive index !!!!!!!
clear idx_dives
clear idx





