addpath(genpath('humans/functions/external/fooof_mat-main'))
addpath('humans/functions')


human_datadir = 'humans/data';

%% Load the data
load(fullfile(human_datadir, 'macrocontacts/ER20_macro1_HFO_snippets.mat'), 'HFOs', 'snippets', 'tt')
macroHFOs = HFOs; clear HFOs;
snippets_macro = snippets'; clear snippets;

load(fullfile(human_datadir, 'ripmap labels/ER20_macro1_ripmap_labels.mat'), 'labels')
macro_labels = labels; clear labels;
%%
load('humans/Figure 5/NY20_micro_all.mat')

%% microwire snippets around macro-ripples

win = [0.1 0.1];
srLFP = 2048;
freqs = 1:1:200;
threshold = 0.2;

all_chan = struct();
all_peaks = zeros(length(macroHFOs.peaks), 8);
all_amps = zeros(length(macroHFOs.peaks), 8);

for i = 1:8 % loop over channels
    this_lfp = micro_all(:,i);
    snippets = get_snippets (this_lfp, win, macroHFOs.peaks*srLFP,srLFP, 0)';

    [powspctrm] = pwelch(snippets, 100, [], freqs, srLFP);

    [stats] = ripple_stats_JS(powspctrm, freqs, threshold);
    rip_peaks = [stats.max];
    
    rip_amps = zeros(size(rip_peaks)); % amplitude of the rip_peaks
    for j = 1:length(stats)
        if rip_peaks(j) == 0
            rip_amps(j) = NaN;
        else
        amp = stats(j).amp;
        pk_freqs = stats(j).freqs;

        idx = ismember(pk_freqs, rip_peaks(j));
        rip_amps(j) = amp(idx);
        end
    end

    all_peaks(:,i) = rip_peaks';
    all_amps(:,i) = rip_amps';  

    all_chan(i).snippets = snippets;
    all_chan(i).powspctrm = powspctrm;
    all_chan(i).avg_all = mean(powspctrm,2);

end

%% # of microwire channels with ripple peak
rof = [60 150]; % 60-150 Hz

map = (all_peaks >= rof(1) & all_peaks < rof(2)); % ripple peak events

num_rip_chans = sum(map, 2);

num_coupled_events = sum(map,1);

% choose the max amplitude of ripple channel if any

max_amp_chan = zeros(size(all_amps,1),1);
for i = 1:size(all_amps,1)
    this_row_map = map(i,:);

    if sum(this_row_map) == 0 % no ripple chan
        [~ , max_amp_chan(i)] = max(all_amps(i,:));
    else % choose the strongest among ripple chans
        [~ , idx] = max(all_amps(i,this_row_map));
        rip_chans = find(this_row_map);
        max_amp_chan(i) = rip_chans(idx);
    end
end

max_chan_freqs = all_peaks(sub2ind(size(all_peaks), (1:length(all_peaks))', max_amp_chan));

max_chan_amps = all_amps(sub2ind(size(all_amps), (1:length(all_amps))', max_amp_chan));

% remove amps from non-ripple channels
max_chan_amps(num_rip_chans == 0) = NaN;

%% # of microwire channels with a ripple peak

figure(1), clf

histogram(num_rip_chans(macro_labels==1))

ylabel('Macro-Ripple #'); xlabel('# of coupled microwire channels')

%% average macroSWR waveform by number of microwire channels

figure(2), clf 
for i = 1:9
subplot(3,3,i)
plot(tt, mean(snippets_macro(:,(num_rip_chans == i-1) & macro_labels==1),2))
xlim([-0.1 0.1])
title({'# of ripple microwire chan: ', num2str(i-1)})

end

%% peak and amplitude distribution
[powspctrm_macro] = pwelch(snippets_macro, 100, [], freqs, srLFP);

threshold  = 0.2;
[stats_macro, corrected_macro] = ripple_stats_JS(powspctrm_macro, freqs, threshold);

pks = [stats_macro.max];

rip_amps = zeros(size(pks));
for j = 1:length(stats_macro)

    if pks(j) == 0
        rip_amps(j) = NaN;
    else
        amp = stats_macro(j).amp;
        pk_freqs = stats_macro(j).freqs;

        idx = ismember(pk_freqs, pks(j));
        rip_amps(j) = amp(idx);
    end

end

%% Peak Frequency histogram
figure(3), clf
for i = 1:9
    subplot(3,3,i)
    histogram(pks((num_rip_chans == i-1) & macro_labels==1),'BinWidth',5);

    xlim([10 160])
    xlabel('Frequency (Hz)'); ylabel('MacroSWR#')
    title({'# of ripple microwire chan: ', num2str(i-1)})
    ylim([0 400])
end


%% Amplitude histogram
figure(4),clf
for i = 1:9
    subplot(3,3,i)
    histogram(rip_amps((num_rip_chans == i-1) & macro_labels==1), 'BinWidth',0.1);

    xlabel('Amplitude'); ylabel('MacroSWR#')
    xlim([0 3]); ylim([0 200])
    title({'# of ripple microwire chan: ', num2str(i-1)})
    hold on

end