%plot CSD profiles of DS
%show histogram of sinks
%allows selection of sinks using mouse
%2021 Dino Dvorak, New York University, neuraldino@gmail.com
%Dentate spikes and external control of hippocampal function,
%Cell Reports, Volume 36, Issue 5, 2021,109497,ISSN 2211-1247,
%https://doi.org/10.1016/j.celrep.2021.109497

clear; close all; clc;

drCSD = 'DS_CSD/';

%input LFP file (channels x samples)
fnEEG = '2016_09_30_PP12_334_PRETRAIN_BASE_eeg.mat';

load([drCSD fnEEG])

figure('Position',[10 100 1000 800])

plotAboveBelow = 1; %1=above,-1=below

minmax = [-1e4 2e4]; %for CSD plot
%minmax = [-3e4 4e4];

%set sinks to nan for first run
%sinks = nan(1,4);

%edit sinks when you localize sink bands
sinks = [104,118,177,190];
    
subplot(1,6,1:3);
imagesc(resultsCSD,minmax)

colormap(jet);

h = title(fnEEG);
set(h,'Interpreter','none');

maxes = [];
mins = [];
    
for i = 1:size(resultsCSD,2)

    x = resultsCSD(:,i)';

    %sources
    kMax = find((x > [x(1) x(1:(end-1))]) & (x >= [x(2:end) x(end)])); 
    kMax = kMax(kMax > 1 & kMax < 200);
    xMax = x(kMax);

    %sinks
    kMin = find((x < [x(1) x(1:(end-1))]) & (x <= [x(2:end) x(end)]));
    kMin = kMin(kMin > 1 & kMin < 200);
    xMin = x(kMin);

    ma = cat(1,kMax,xMax); %concatenate
    mi = cat(1,kMin,xMin); %concatenate

    maxes = cat(2, maxes, ma);
    mins = cat(2, mins, mi);

end

edges = 0.5:1:200.5;

hma = histc(maxes(1,:),edges); %hist of locations
hmi = histc(mins(1,:),edges); %hist of locations

subplot(1,6,4);
hold on;

%plot(hma,edges+0.5,'r');
plot(hmi,edges+0.5,'b');

for i = 1:length(sinks)
    x = sinks(i);
    if isnan(x); continue; end;
    plot(hmi(x),x,'r.','MarkerSize',20);
end

axis ij tight


%CSD average based on top 2 sinks
if ~isnan(sinks(1)) && ~isnan(sinks(2))    
    kKeep = false(2,size(resultsCSD,2));
    for i = 1:size(resultsCSD,2)
        x = resultsCSD(:,i)';
        %sinks
        kMin = find((x < [x(1) x(1:(end-1))]) & (x <= [x(2:end) x(end)]));
        kMin = kMin(kMin > 1 & kMin < 200);
        if sum(kMin > sinks(1)-4 & kMin < sinks(1)+4) > 0 %DS1
            kKeep(1,i) = 1;
        end
        if sum(kMin > sinks(2)-4 & kMin < sinks(2)+4) > 0 %DS2
            kKeep(2,i) = 2;
        end
    end

    subplot(1,6,5); hold on;
    %DS1
    k = kKeep(1,:);
    csd = resultsCSD(:,k);
    csd = nanmean(csd,2);
    plot(csd,1:200,'r');
    %DS2
    k = kKeep(2,:);
    csd = resultsCSD(:,k);
    csd = nanmean(csd,2);
    plot(csd,1:200,'k');

    axis ij tight;
end %top inks

%CSD average based on bottom 2 sinks
if ~isnan(sinks(3)) && ~isnan(sinks(4))    
    kKeep = false(2,size(resultsCSD,2));
    for i = 1:size(resultsCSD,2)
        x = resultsCSD(:,i)';
        %sinks
        kMin = find((x < [x(1) x(1:(end-1))]) & (x <= [x(2:end) x(end)]));
        kMin = kMin(kMin > 1 & kMin < 200);
        if sum(kMin > sinks(3)-4 & kMin < sinks(3)+4) > 0 %DS1
            kKeep(1,i) = 1;
        end
        if sum(kMin > sinks(4)-4 & kMin < sinks(4)+4) > 0 %DS2
            kKeep(2,i) = 2;
        end
    end

    subplot(1,6,6); hold on;
    %DS1
    k = kKeep(2,:);
    csd = resultsCSD(:,k);
    csd = nanmean(csd,2);
    plot(csd,1:200,'r');
    %DS2
    k = kKeep(1,:);
    csd = resultsCSD(:,k);
    csd = nanmean(csd,2);
    plot(csd,1:200,'k');

    axis ij tight;
end %top inks