
%% Default values
p = inputParser;
addParameter(p,'basepath',pwd,@isdir);
addParameter(p,'rippleChannel',[],@isnumeric);
addParameter(p,'SWChannel',[],@isnumeric);
addParameter(p,'thresholds',[1.5 3.5],@isnumeric); 
addParameter(p,'SWthresholds',[-0.5 -2], @isnumeric);
addParameter(p,'durations',[30 100],@isnumeric); % 'anna:' minimum inter-ripple duration, maximum ripple duration , maybe change to 20 300
addParameter(p,'restrict',[],@isnumeric); 
addParameter(p,'frequency',[],@isnumeric);
addParameter(p,'stdev',[],@isnumeric);
addParameter(p,'show','off',@isstr);
addParameter(p,'noise',[],@ismatrix);
addParameter(p,'passband',[120 200],@isnumeric);
addParameter(p,'SWpassband',[2 10],@isnumeric);
addParameter(p,'EMGThresh',1,@isnumeric);
addParameter(p,'saveMat',true,@islogical);
addParameter(p,'minDuration',20,@isnumeric);
addParameter(p,'plotType',2,@isnumeric);
addParameter(p,'srLfp',[],@isnumeric); %'anna' human data lfp is 2048
addParameter(p,'rippleStats',true,@islogical);
addParameter(p,'debug',false,@islogical);
addParameter(p,'eventSpikeThreshold',1,@isnumeric);
addParameter(p,'force',false,@islogical);
addParameter(p,'removeRipplesStimulation',true,@islogical);
addParameter(p,'useCSD',false,@islogical);
addParameter(p,'stdThreshold',1.5,@isnumeric);

parse(p,varargin{:})

basepath = p.Results.basepath;
rippleChannel = p.Results.rippleChannel;
SWChannel = p.Results.SWChannel;
thresholds = p.Results.thresholds;
SWthresholds = p.Results.SWthresholds;
durations = p.Results.durations;
restrict = p.Results.restrict;
frequency = p.Results.frequency;
stdev = p.Results.stdev;
show = p.Results.show;
noise = p.Results.noise;
passband = p.Results.passband;
SWpassband = p.Results.SWpassband;
EMGThresh = p.Results.EMGThresh;
saveMat = p.Results.saveMat;
minDuration = p.Results.minDuration;
plotType = p.Results.plotType;
srLfp = p.Results.srLfp;
rippleStats = p.Results.rippleStats;
debug = p.Results.debug;
eventSpikeThreshold = p.Results.eventSpikeThreshold;
force = p.Results.force;
removeRipplesStimulation = p.Results.removeRipplesStimulation;
useCSD = p.Results.useCSD;
stdThreshold = p.Results.stdThreshold;

%% Load Session Metadata and several variables if not provided
% session = sessionTemplate(basepath,'showGUI',false);
session = loadSession(basepath);

if (exist([session.general.name '.ripples.events.mat'],'file') ...
        && ~force)
    disp(['Ripples already detected for ', session.general.name, '. Loading file.']);
    load([session.general.name '.ripples.events.mat']);
    return
end

% Ripple and SW Channel are loaded separately in case we want to provide
% only one of the
if isempty(rippleChannel)
    if ~isempty(dir([session.general.name,'.hippocampalLayers.channelinfo.mat']))
        file = dir([session.general.name,'.hippocampalLayers.channelinfo.mat']);
        load(file.name);
    else
        [hippocampalLayers] = getHippocampalLayers();
    end
    rippleChannel = hippocampalLayers.bestShankLayers.pyramidal;
end

if isempty(SWChannel)
    if ~isempty(dir([session.general.name,'.hippocampalLayers.channelinfo.mat']))
        file = dir([session.general.name,'.hippocampalLayers.channelinfo.mat']);
        load(file.name);
    else
        [hippocampalLayers] = getHippocampalLayers();
    end
    SWChannel = hippocampalLayers.bestShankLayers.radiatum;
end

if removeRipplesStimulation && ~isempty(dir('*optogeneticPulses.events.mat'))
    targetFile = dir('*optogeneticPulses.events.mat'); load(targetFile.name);
    restrict_temp = SubtractIntervals([0 Inf],optoPulses.stimulationEpochs);
    restrict =  ConsolidateIntervals([restrict; restrict_temp]);
end

if useCSD
    disp('Computing CSD...');
    rippleChannel = computeCSD([],'channels',rippleChannel);
end

if isempty(frequency)
    frequency = session.extracellular.srLfp;
end



%% The code starts here

cutwind = 50; % number of datapoints taken before and after the event

[Epi]= findRipples('rippleChannel', 11,'thresholds',[1.5 3.5],'passband', [120 200],...
    'EMGThresh',EMGThresh,'durations',[30 200], 'saveMat',false,'frequency',2048);
%humfindRipples_a(ripplech,'thresholds',[3 20],'passband',[50 80],'durations',[250 500], 'saveMat',false,'frequency',2048);

%durations 1 is th emax duration of the event, durations 2 is the min
%distance between 2 events

%figure; subplot(1,2,1); hist(Epi.peakNormedPower,100);
%  h = figure;
%         hist(Epi.peakNormedPower,100);
%         title('Select threshold with the mouse and press right click...');
%         [cutpoint,thr] = ginput(1);
%         hold on
%         line(cutpoint*[1 1],ylim,'Color','r');
%         pause(1);
%         close(h);
% validEvents = find(Epi.peakNormedPower>=cutpoint);
validEvents = find(Epi.peakNormedPower>=5);
Epi_Events.timestamps = Epi.timestamps(validEvents,:);
try 
   Epi_Events.peaks = Epi.peaks(validEvents,:);
   Epi_Events.peakNormedPower = Epi.peakNormedPower(validEvents,:);
end


%maxvalues = round(Epi_Events.peaks*1250); %maxvalues = (Epi_Events.peaks*1250);

%fprintf ('creating index...\n')
initial_index = zeros(1,length(Epi_Events.peaks));
single_events = zeros(2*cutwind+1,length(Epi_Events.peaks));
single_lowpass = zeros(2*cutwind+1,length(Epi_Events.peaks));
for yy = 1:length(Epi_Events.peaks)
initial_index(yy)=find(Epi.lfp.timestamps == Epi_Events.peaks(yy));
single_events(:,yy) = Epi.lfp.data((initial_index(yy)-cutwind):(initial_index(yy)+cutwind));
single_lowpass(:,yy) = Epi.signal((initial_index(yy)-cutwind):(initial_index(yy)+cutwind)); 
end

[M,I] = max(single_events); %maximum values and indices of single events
index = (initial_index - cutwind)+I-1; %creates an index of the peaks

% remove big distances between initially registered and real
h = figure;
hist(I,100);
title('Select threshold with the mouse and press right click...');
[cutind1,thr1] = ginput(1);
hold on
line(cutind1*[1 1],ylim,'Color','r');
title('Select threshold with the mouse and press right click...');
[cutind2,thr2] = ginput(1);
hold on
line(cutind2*[1 1],ylim,'Color','c');
pause(1);
close(h);
A = find(I>=cutind1 & I<=cutind2);
index = index(A);
initial_index = initial_index(A);nts.peakNormedPower = Epi.peakNormedPower(A,:);


Peaks_timestamps = Epi.lfp.timestamps(index);
Epi_Events.peaks = Epi_Events.peaks(A);
Epi_Events.timestamps = Epi_Events.timestamps(A,:);
Epi_Events.peakNormedPower = Epi_Events.peakNormedPower(A);
mkdir('Analysis')
% plot a 3x3 figure with individual events

% if size(A,2) >= 36 
%     z = 3;
% else z = fix(size(A,2)/9);
% end
% 
% for ii = 1:9:9*z
% figure;
% multiplot;
% title("session " + figure_title(end-1:end) + "")
% saveas(gcf,'Analysis\multiplot1.fig');  
% saveas(gcf,'Analysis\multiplot1.jpeg'); 
% end



for ii = 1:9:9*(fix(size(A,2)/9))
figure;
multiplot;
prompt = 'Do you want more? Y/N [Y][N]: ';
str = input(prompt, 's');
if isempty(str)
    str = 'Y';
else str = 'N';
end
    if str == 'Y';
      
        continue
    else
break
    end
end

%plotting the whole trace;

figure;
plot (Epi.lfp.timestamps, Epi.lfp.data);
hold on
plot (Epi.lfp.timestamps, Epi.signal);
plot(Peaks_timestamps,Epi.lfp.data(index),'o');
title("session " + figure_title(end-1:end) + "")
% scatter(Epi_Events.timestamps(:,1),Epi_Events.peakNormedPower);
% scatter(Epi_Events.timestamps(:,2),Epi_Events.peakNormedPower);
xlabel ('s');
axis tight

%mkdir('Analysis')
saveas(gcf,'Analysis\IED_sess.fig');  
saveas(gcf,'Analysis\IED_sess.jpeg');  

% wavelet;
lfp = Epi.lfp;
% if B;
%   [wavAvg, lfpAvg ] = bz_eventWavelet (lfp, Peaks_timestamps (B), 'twin',[.3 .3], 'frange',[10 1500], 'whitening',false);
% caxis([1000 20000]);  
% else
if size(Peaks_timestamps) >= 50;
[wavAvg, lfpAvg ] = bz_eventWavelet (lfp, Peaks_timestamps (1:50), 'twin',[.3 .3], 'frange',[6 90], 'whitening',false);
caxis([1000 30000])
else
[wavAvg, lfpAvg ] = bz_eventWavelet (lfp, Peaks_timestamps, 'twin',[.3 .3], 'frange',[10 90], 'whitening',false);
caxis([1000 20000])
end
title (" session " + figure_title(end-1:end) + " f_m_a_x " +wavAvg.maxFreq+ " pow_m_a_x " +wavAvg.maxPow+ "")
saveas(gcf,'Analysis\IED.fig');  
saveas(gcf,'Analysis\IED.jpeg');  

% saving
Epi_Events.peaks=[];
Epi_Events.peaks=Peaks_timestamps;

filename = split(pwd,filesep); filename = filename{end};
save([filename '.IED.mat'],'Epi_Events');



