set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 12);
set(0,'defaultAxesFontSize',12);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
set(0, 'DefaultAxesFontName', 'Arial');
set(0, 'DefaultLineLineWidth', 2);

%%
clear all
clc

animalNames = {'OXYT025';
    'OXYT026';
    'OXYT027'};

%
fs = 20e3;
fs_dff = 100;

CNOTimes= [2 3]*3600;

RipPowerCNO = zeros([7200,20]);
thetaPowerCNO = zeros([7200,20]);
animalID = [];
sessionID =[];
filenamesID =[];
sessionCount = 0;
rippleDurs = {};

[Bf,Af] = butter(3,[100 200]/(1250/2));
for i = 1:length(animalNames)
    Names_CNO = dir(['CNO/',animalNames{i},'*CNO.mat']);
    filenamesCNO = {Names_CNO.name};
    foldernameCNO = {Names_CNO.folder};
    for cnID = 1:size(Names_CNO,1)
        sessionCount = sessionCount + 1
        load(fullfile(foldernameCNO{cnID},filenamesCNO{cnID}));
        animalID =[animalID; convertCharsToStrings(animalNames{i})];
        sessionID =[sessionID; "CNO"];
        filenamesID = [filenamesID;convertCharsToStrings(filenamesCNO{cnID}) ];
        if ~exist('OXT_CNO.lfp_fbtime','var')
            Tstart = 0;
        else
            Tstart = OXT_CNO.lfp_fbtime(1)/fs;
        end
        
        ripples= findImmobileRipple(OXT_CNO.ripples,OXT_CNO.acel);
        ripLFPfiltered = filtfilt(Bf,Af,OXT_CNO.RipLFP.data);
        [maps,data,stats] = bz_RippleStats(ripLFPfiltered,OXT_CNO.RipLFP.timestamps,OXT_CNO.ripples);
        rippleFrequency(sessionCount,1)  = mean(data.peakFrequency);

        ripple_times = ripples.peaks -Tstart;
        tCNOIds = find((ripple_times>CNOTimes(1)).*(ripple_times<CNOTimes(2)));
        ripp_diff = diff(ripple_times);
        IRI = ripp_diff(tCNOIds);
        immobilityDur(sessionCount,1) = sum(ripples.immobile(CNOTimes(1):CNOTimes(2)));
        IRImean(sessionCount,1) = mean(IRI);
        IRIstd(sessionCount,1) = std(IRI);
        ripplePeakPowerMean(sessionCount,1) = mean( ripples.peakNormedPower(tCNOIds));
        ripplePeakPowerMax(sessionCount,1) = max( ripples.peakNormedPower(tCNOIds));
        ripplePeakPowerStd(sessionCount,1) = std( ripples.peakNormedPower(tCNOIds));
        rippleNumbers(sessionCount,1) = length( ripples.peakNormedPower(tCNOIds));
        rippleDurations = ripples.durations;
        rippleDurMean(sessionCount,1) = mean(rippleDurations(tCNOIds));
        rippleDurStd(sessionCount,1) = std(rippleDurations(tCNOIds));
        rippleDursCNO{sessionCount} = rippleDurations(tCNOIds);

        LFPs = OXT_CNO.RipLFP.data;
        [spectrogrammt,tm,fm] = MTSpectrogram(LFPs,'window',2);
        Ntm = length(tm);
        bands = SpectrogramBands(spectrogrammt,fm,'theta', [5,10]);
        tCNOIds = find(((tm>CNOTimes(1)).*(tm<CNOTimes(2))).*(ripples.immobile(1:Ntm))');
        Ripp_power = (bands.ripples);
        %RipPowerCNO(:,sessionCount) = Ripp_power(tCNOIds);

        LFPs = OXT_CNO.THLFP.data;
        [spectrogrammt,tm,fm] = MTSpectrogram(LFPs,'window',2);
        bands = SpectrogramBands(spectrogrammt,fm,'theta', [5,10]);
        thetaPower = (bands.theta);
        %thetaPowerCNO(:,sessionCount) = thetaPower(tCNOIds);

        meanThetaPower(sessionCount,1) = mean(thetaPower(tCNOIds),1);
        meanRipplePower(sessionCount,1) = mean(Ripp_power(tCNOIds),1);



    end

    Names_Control = dir(['Control/',animalNames{i},'*CNO_Ctrl.mat']);
    filenamesCtrl = {Names_Control.name};
    foldernameCtrl = {Names_Control.folder};

    for ctID = 1:size(Names_Control,1)
        sessionCount = sessionCount + 1
        load(fullfile(foldernameCtrl{ctID},filenamesCtrl{ctID}));
        animalID =[animalID; convertCharsToStrings(animalNames{i})];
        sessionID =[sessionID; "Ctrl"];
        filenamesID = [filenamesID;convertCharsToStrings(filenamesCtrl{ctID}) ];
        if ~exist('OXT_CNO_Ctrl.lfp_fbtime','var')
            Tstart = 0;
        else
            Tstart = OXT_CNO.lfp_fbtime(1)/fs;
        end


        ripples= findImmobileRipple(OXT_CNO_Ctrl.ripples,OXT_CNO_Ctrl.acel);
        ripLFPfiltered = filtfilt(Bf,Af,OXT_CNO_Ctrl.RipLFP.data);
        [maps,data,stats] = bz_RippleStats(ripLFPfiltered,OXT_CNO_Ctrl.RipLFP.timestamps,OXT_CNO_Ctrl.ripples);
        rippleFrequency(sessionCount,1)  = mean(data.peakFrequency);
        ripple_times = ripples.peaks -Tstart;

        CtrlTimes= [0 length(OXT_CNO_Ctrl.Oxt_dFF)/100];
        tCtrlIds = find((ripple_times>CtrlTimes(1)).*(ripple_times<CtrlTimes(2)));

        IRI = diff(ripple_times(tCtrlIds));
        immobilityDur(sessionCount,1) = ripples.durImmobile;
        IRImean(sessionCount,1) = mean(IRI);
        IRIstd(sessionCount,1) = std(IRI);
        ripplePeakPowerMean(sessionCount,1) = mean( ripples.peakNormedPower(tCtrlIds));
        ripplePeakPowerMax(sessionCount,1) = max( ripples.peakNormedPower(tCtrlIds));
        ripplePeakPowerStd(sessionCount,1) = std( ripples.peakNormedPower(tCtrlIds));
        rippleNumbers(sessionCount,1) = length( ripples.peakNormedPower(tCtrlIds));
        rippleDurations = ripples.durations;
        rippleDurMean(sessionCount,1) = mean(rippleDurations(tCtrlIds));
        rippleDurStd(sessionCount,1) = std(rippleDurations(tCtrlIds));
        rippleDurs{sessionCount} = rippleDurations(tCtrlIds);

        LFPs = OXT_CNO_Ctrl.RipLFP.data;
        [spectrogrammt,tm,fm] = MTSpectrogram(LFPs,'window',2);
        bands = SpectrogramBands(spectrogrammt,fm,'theta', [5,10]);
        Ripp_power = (bands.ripples);
        %RipPowerCtrl(:,sessionCount) = Ripp_power(tCNOIds);

        LFPs = OXT_CNO_Ctrl.THLFP.data;
        [spectrogrammt,tm,fm] = MTSpectrogram(LFPs,'window',2);
        bands = SpectrogramBands(spectrogrammt,fm,'theta', [5,10]);
        thetaPower = (bands.theta);
        %thetaPowerCtrl(:,sessionCount) = thetaPower(tCNOIds);

        Ntm = length(tm);
        meanThetaPower(sessionCount,1) = mean(thetaPower(ripples.immobile(1:Ntm)),1);
        meanRipplePower(sessionCount,1) = mean(Ripp_power(ripples.immobile(1:Ntm)),1);
    end
end

%%
rippleRate = rippleNumbers./immobilityDur;
RippleCNO_Ctrl = table(animalID,sessionID,filenamesID,meanThetaPower,meanRipplePower,IRImean,IRIstd,ripplePeakPowerMean,...
    ripplePeakPowerMax,ripplePeakPowerStd,rippleNumbers,rippleDurMean,rippleDurStd,immobilityDur,rippleRate,rippleFrequency);
 writetable(RippleCNO_Ctrl,'RippleCNO_Ctrl.xlsx','Sheet',1);
%%
figure
for i =1:4
    subplot(4,1,i)
    histogram((RipPowerCNO(:,i)),200,'Normalization','pdf')
end

%%
unID= unique(RippleCNO_Ctrl.animalID);
markr = ['p','+','d','s'];
figure
subplot(1,4,1)
for i =1:length(unID)
    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(2*ones([ height(tfiltered),1]),tfiltered.rippleRate,100,markr(i),'MarkerEdgeColor',[0,0.7,0]); hold on

    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(1*ones([ height(tfiltered),1]),tfiltered.rippleRate,100,markr(i),'MarkerEdgeColor',[0,0.5,0]); hold on
    xlim([0.5,2.5])
    ylim([0.3,1.2])
    i
end
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') ,:);
boxchart(2*ones([ height(tfiltered),1]),(tfiltered.rippleRate)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.25)
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') ,:);
boxchart(1*ones([ height(tfiltered),1]),(tfiltered.rippleRate)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.5)
title('Ripple Rate')



subplot(1,4,2)
for i =1:length(unID)
    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(2*ones([ height(tfiltered),1]),tfiltered.meanRipplePower,100,markr(i),'MarkerEdgeColor',[0,0.7,0]); hold on

    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(1*ones([ height(tfiltered),1]),tfiltered.meanRipplePower,100,markr(i),'MarkerEdgeColor',[0,0.5,0]); hold on
    xlim([0.5,2.5])
    i
end
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') ,:);
boxchart(2*ones([ height(tfiltered),1]),(tfiltered.meanRipplePower)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.25)
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') ,:);
boxchart(1*ones([ height(tfiltered),1]),(tfiltered.meanRipplePower)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.5)
title('mean Ripple Power')

subplot(1,4,4)
for i =1:length(unID)
    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(2*ones([ height(tfiltered),1]),tfiltered.meanThetaPower,100,markr(i),'MarkerEdgeColor',[0,0.7,0]); hold on

    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(1*ones([ height(tfiltered),1]),tfiltered.meanThetaPower,100,markr(i),'MarkerEdgeColor',[0,0.5,0]); hold on
    xlim([0.5,2.5])
    i
end
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') ,:);
boxchart(2*ones([ height(tfiltered),1]),(tfiltered.meanThetaPower)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.25)
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') ,:);
boxchart(1*ones([ height(tfiltered),1]),(tfiltered.meanThetaPower)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.5)
title('mean Theta Power')

% subplot(1,4,4)
% for i =1:length(unID)
%     tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') & (RippleCNO_Ctrl.animalID==unID(i)),:);
%     scatter(2*ones([ height(tfiltered),1]),tfiltered.rippleDurMean,100,markr(i),'r'); hold on
%     tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') & (RippleCNO_Ctrl.animalID==unID(i)),:);
%     scatter(1*ones([ height(tfiltered),1]),tfiltered.rippleDurMean,100,markr(i),'b'); hold on
%     xlim([0.5,2.5])
%     i
% end
% tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') ,:);
% boxchart(2*ones([ height(tfiltered),1]),(tfiltered.rippleDurMean)')
% tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') ,:);
% boxchart(1*ones([ height(tfiltered),1]),(tfiltered.rippleDurMean)')
% title('ripple Dur Mean')


subplot(1,4,3)
for i =1:length(unID)
    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(2*ones([ height(tfiltered),1]),tfiltered.rippleFrequency,100,markr(i),'MarkerEdgeColor',[0,0.7,0]); hold on

    tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') & (RippleCNO_Ctrl.animalID==unID(i)),:);
    scatter(1*ones([ height(tfiltered),1]),tfiltered.rippleFrequency,100,markr(i),'MarkerEdgeColor',[0,0.5,0]); hold on
    xlim([0.5,2.5])
    i
end
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='CNO') ,:);
boxchart(2*ones([ height(tfiltered),1]),(tfiltered.rippleFrequency)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.25)
tfiltered = RippleCNO_Ctrl((RippleCNO_Ctrl.sessionID =='Ctrl') ,:);
boxchart(1*ones([ height(tfiltered),1]),(tfiltered.rippleFrequency)', 'BoxFaceColor',[0,0.5,0],'BoxFaceAlpha',0.5)
title('ripple Frequency')
exportgraphics(gcf,'OXT_CNORippleStats.pdf','BackgroundColor','none','ContentType','vector')

%%

ripDurCNO = [];
for i = 1:length(rippleDursCNO)
    ripDur = rippleDursCNO{i};
    if ~isempty(ripDur)
        ripDurCNO = [ripDurCNO;ripDur*1e3 ];
    end
end

edges = linspace(20,120,50);
bins = 0.5*(edges(1:end-1) + edges(2:end));
figure
N = smoothdata(histcounts(ripDurCNO,edges));
semilogx(bins,(N/sum(N)),'r'); hold on
xline(median(ripDurCNO),'--r')
ripDurCtrl = [];
for i = 1:length(rippleDurs)
    ripDur = rippleDurs{i};
    if ~isempty(ripDur)
        ripDurCtrl = [ripDurCtrl;ripDur*1e3 ];
    end
end

N = smoothdata(histcounts(ripDurCtrl,edges));
semilogx(bins,(N/sum(N)),'k')
xline(median(ripDurCtrl),'--k')
box off
legend('CNO','CNO median','Ctrl', 'Ctrl median')


