%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function iterates through a set of preprocessed traces, plots them
% as an overview, with zoom in on representative sections throughout the
% recordings: baseline, first hour after PTZ and second hour after PTZ. 
% Then it extracts the peri-stimulus time histogram for high amplitude and
% low amplitude events. 
%
% The preprocessed traces should be arranged in the following manner: in
% the current folder, each session in a separate folder, named Session 
% <session_number>, each trace in a separate file with the last two digits 
% representing the recording number.
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;

% TO START, create a file structure with the datasets:
%
% file_struct(file_number, feature_number) = [session_number recording_number
% baseline/response control/target dpf errors]
% 
% Note that:
% session_number = number of the session - should be present in the
% folder name;
% recording_number = number of the recording - should be present in the
% last digits of the recording name;
% baseline/response = 0 if the trace is recorded only during spontaneous
% activity; 1 if trace contains response to PTZ 
% control/target = 0 if the trace is recorded from a Mismatch control fish;
% 1 if trace is recorded from a DEPDC5 animal

file_struct(1,:) = [12 2 1 0 5 0];
file_struct(2,:) = [15 3 1 1 5 0];

% set the initial parameters
f_sample = 250;     % [Hz] = sampling frequency of the trace

% select datasets for mismatch control and DEPDC5
files_mis = file_struct(find(file_struct(:,4)==0),:);
files_dep = file_struct(find(file_struct(:,4)==1),:);

% set the moment of PTZ application (can be set automatically by detecting
% the artefact induced by PTZ application)
onset_mis = 3642; % seconds
onset_dep = 4100; % seconds

%% 
% select traces for mismatch control and DEPDC5
[V_base_mis] = select_traces(files_mis, f_sample);
[V_base_dep] = select_traces(files_dep, f_sample);

% plot full traces
figure;
subplot(2,1,1);
plot(V_base_mis.t_base-onset_mis, V_base_mis.trace, 'k');
ylim([-8 3]);
xlim([-2000 5500]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['MIS']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);
title('Typical traces MIS vs DEP');

subplot(2,1,2);
plot(V_base_dep.t_base-onset_dep, V_base_dep.trace, 'k');
ylim([-8 3]);
xlim([-2000 5500]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['DEP']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);


% ZOOM on baseline period
figure;
subplot(2,1,1);
plot(V_base_mis.t_base-onset_mis, V_base_mis.trace, 'k');
ylim([-1 1]);
xlim([-350 -250]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['MIS']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);
title('Zoom baseline traces MIS vs DEP');

subplot(2,1,2);
plot(V_base_dep.t_base-onset_dep, V_base_dep.trace, 'k');
ylim([-1 1]);
xlim([-350 -250]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['DEP']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);

% ZOOM on first PTZ response period
figure;
subplot(2,1,1);
plot(V_base_mis.t_base-onset_mis, V_base_mis.trace, 'k');
ylim([-6 3]);
xlim([1700 2200]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['MIS']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);
title('Zoom high amp traces MIS vs DEP');

subplot(2,1,2);
plot(V_base_dep.t_base-onset_dep, V_base_dep.trace, 'k');
ylim([-6 3]);
xlim([2200 2700]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['DEP']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);

% ZOOM on second PTZ response period
figure;
subplot(2,1,1);
plot(V_base_mis.t_base-onset_mis, V_base_mis.trace, 'k');
ylim([-1 1]);
xlim([3570 3670]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['MIS']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);
title('Zoom low amp traces MIS vs DEP');

subplot(2,1,2);
plot(V_base_dep.t_base-onset_dep, V_base_dep.trace, 'k');
ylim([-1 1]);
xlim([3500 3600]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['DEP']);
l = line([0 0],[-8,3]);
set(l, 'color', 'r', 'linewidth', 10);

%% PSTH
% set initial aparameters
binsize = 200;          % seconds = length of PSTH bin
threshold_amp = 1;      % mV = threshold value for separating high amplitude versus low amplitude events
direction = 'all';      % selects which events to count: 'low' / 'high' / 'all'
rec_length = 170*60;    % s = how long to compute the PSTH  

% compute PSTH for all events
[avg_psth_mis, bins_mis] = compute_count_time(threshold_amp, direction, f_sample, binsize, rec_length, files_mis); 
[avg_psth_dep, bins_dep] = compute_count_time(threshold_amp, direction, f_sample, binsize, rec_length, files_dep); 

% plot PSTH for all events
figure;
subplot(2,1,1);
bar(bins_mis-onset_mis, avg_psth_mis, 'k');
hold on;
ylim([0 10]);
xlim([-2000 5500]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['MIS']);
l = line([0 0],[0 15]);
set(l, 'color', 'r', 'linewidth', 10);
title('PSTH overlay MIS vs DEP'); 

subplot(2,1,2);
bar(bins_dep-onset_dep, avg_psth_dep, 'k');
hold on;
ylim([0 10]);
xlim([-2000 5500]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['DEP']);
l = line([0 0],[0 15]);
set(l, 'color', 'r', 'linewidth', 10);

% compute PSTH for high amplitude events
direction = 'high';
[avg_psth_mis, bins_mis] = compute_count_time(threshold_amp, direction, f_sample, binsize, rec_length, files_mis); 
[avg_psth_dep, bins_dep] = compute_count_time(threshold_amp, direction, f_sample, binsize, rec_length, files_dep); 

% overlay high amplitude PSTH over all events PSTH
subplot(2,1,1);
bar(bins_mis-onset_mis, avg_psth_mis,'FaceColor',[57 185 240]/255,'EdgeColor',[57 185 240]/255);
ylim([0 10]);
xlim([-2000 5500]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['MIS']);
l = line([0 0],[0 15]);
set(l, 'color', 'r', 'linewidth', 10);

subplot(2,1,2);
bar(bins_dep-onset_dep, avg_psth_dep ,'FaceColor',[57 185 240]/255,'EdgeColor',[57 185 240]/255);
ylim([0 10]);
xlim([-2000 5500]);
xlabel(['time [s]']);
ylabel(['Amplitude [mV]']);
title(['DEP']);
l = line([0 0],[0 15]);
set(l, 'color', 'r', 'linewidth', 10);


