%% Calculate evoked single-unit timepoints in binned time intervals
% Version vb5_NEX_v3: add function to import ChirpWave detection points.
% Version vb5_NEX_v2: option to remove last trial of trigger (catchpoint)
% Version vb5: compatible with OpenEphys
% Raw data can be calculated with A_NSx_EP_20180810_v2b.m or v2a.m file
% version 4: change Y axis of mean plot to Hz 
% Sorted with BOSS first;
%load NEV file for sorted event
%load NS4 file 
clear all, close all;
%% Input Parameters %%
DataSource=1; % 1: BlockRock; 2: OpenEphys 
NEVtrigger=2; % if 2, load EP detection points from EPindex_ChirpWave.mat.
              % if 1, use sorted channel unit as EP trigger point,
              % if 0, use analog input channel as trigger point
  StimID=136; % define the channel of trigger (Sorted channel)
  StimUnitID=1; % define sorted unit ID for trigger
  NEVchannels=[1 2 3 4 5 6 9 9 12 13 14 15 16 136]; % Can be multiple channels, matrix length match NEVunits
  NEVunits=[1 1 1 1 1 1 1 2 1 1 1 1 1 1]; % Corresponding detected units of each channel, can be more than one unit in each channel
  ExtInt=[-50,200]; % in ms, Interval to extra, negative means before the catch point 
  Bin=10; % in ms, devided ExtInt into sections NOTE!!!(ExtInt(2)-ExtInt(1))/Bin should always be scalar.
RemoveLastTrial=0; % if==1 means remove last trial; keep =0 in most case
PlotEPdetection=1; % if==1, plot EP detection points in one EP graph
%%%%%%%%% Put input parameters into structure %%%%%%%%%%%%%
save temp.mat;
a_InputParameters=load('temp.mat');
delete temp.mat;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Main function%% 
%% open NEV file %%
NEV=F_openNEX;
%% NSx SourceFile identification %%
 % identify the Sourefile for NEV offline Sorting
if DataSource==1
NSx=openNSx;
SourceFile=[NSx.MetaTags.Filename NSx.MetaTags.FileExt];
a_InputParameters.SourceFile=SourceFile;
% check if the .nev & .nsx file matched
 if NSx.MetaTags.DataPoints-1<NEV.MetaTags.DataDuration<NSx.MetaTags.DataPoints+1 &...
     NEV.MetaTags.Filename==NSx.MetaTags.Filename;
     if NSx.MetaTags.FileExt=='.ns6';
     else % correct the error of NEV if not processed by .ns6 file.
        NEV.MetaTags.SampleRes=NSx.MetaTags.SamplingFreq;
        NEV.MetaTags.DataDurationSec=NEV.MetaTags.DataDuration/NEV.MetaTags.SampleRes;
     end
 else
    error('NEV & NSx file mismatch');
 end
end
if DataSource==2
   [GetFileName GetFilePath]=uigetfile('EphysOutputOFS*.mat'); % get filename and path
   ValidFileName=GetFileName;
   ValidFileName(end-3:end)=[];
   if NEV.MetaTags.Filename==ValidFileName
   else error('NEV & EphysOutputOFS file mismatch');
   end
end
fs=double(NEV.MetaTags.SampleRes);
%% Stimulation (trigger) identification
if NEVtrigger==1;
   A=find(NEV.Data.Spikes.Electrode==StimID); % find spikes index that match the selected ChannelID
   B=find(NEV.Data.Spikes.Unit(A)==StimUnitID); % find spikes index that match the sorted UnitID. 
   CatchPoint=NEV.Data.Spikes.TimeStamp(A(B)); % find spikes taht match the sorted channel-unit.
end

if NEVtrigger==2;
    a_ChirpParameters=load('EPindex_ChirpWave.mat','a_InputParameters');
    load('EPindex_ChirpWave.mat','EPindex');
    CatchPoint=EPindex;
    clear EPindex;
    if strcmp(a_ChirpParameters.a_InputParameters.SourceFile,SourceFile)==0;
        errordlg('FileName of EP and Chirp not matched','File Error');
        return
    end
end

%% remove last trial
if RemoveLastTrial==1;
    CatchPoint=CatchPoint(1:end-1);
end
%% Detect 
  ExtEdges(:,1)=CatchPoint'+ExtInt(1)*fs/1000;
  ExtEdges(:,2)=CatchPoint'+ExtInt(2)*fs/1000;
%% Create Detection points of each unit %%
EPplotXY=cell(1,length(NEVunits));
for i=length(NEVunits):-1:1;
    SU_Name(i)={sprintf('Chan%d-U%d',NEVchannels(i),NEVunits(i))};
    SUpoints(i,:)=zeros(1,NEV.MetaTags.DataDuration);
    C=find(NEV.Data.Spikes.Electrode==NEVchannels(i)); % find spikes index that match the selected ChannelID
    D=find(NEV.Data.Spikes.Unit(C)==NEVunits(i)); % find spikes index that match the sorted UnitID. 
    E=NEV.Data.Spikes.TimeStamp(C(D)); % find spikes that match the sorted channel-unit.
    SUpoints(i,E)=1;
    for j=length(CatchPoint):-1:1;
        EP_SUpoints(j,:,i)=SUpoints(i,ExtEdges(j,1):ExtEdges(j,2));
    end
    [Y,X]=find(EP_SUpoints(:,:,i)==1); % Y:trials; X:timestamps
    X=ExtInt(1)+(X-1)*1000/fs; % X: timestamps in ms
    EPplotXY{i}=[X,Y]; 
end


if PlotEPdetection==1
  i=length(NEVunits)+1;
  SUpoints(i,CatchPoint)=1;
    for j=length(CatchPoint):-1:1;
        EP_SUpoints(j,:,i)=SUpoints(i,ExtEdges(j,1):ExtEdges(j,2));
    end
    [Y,X]=find(EP_SUpoints(:,:,i)==1); % Y:trials; X:timestamps
    X=ExtInt(1)+(X-1)*1000/fs; % X: timestamps in ms
    EPplotXY{i}=[X,Y]; 
    SU_Name{i}=['EPpoints'];
end
EP_SUmean=squeeze(mean(EP_SUpoints,1));

%Process binned data
BinEdges=1:Bin*fs/1000:([ExtInt(2)-ExtInt(1)]*fs/1000+1);
BinNumber=(ExtInt(2)-ExtInt(1))/Bin;
for k=BinNumber:-1:1;
    EP_SUmeanBinHz(k,:)=sum(EP_SUmean([BinEdges(k):BinEdges(k+1)-1],:))*1000/Bin;
end



%% Plot %%%%%%%%%%%%%%%%% 
%%Scatter plot for EP 
for i=1:length(NEVunits);
    
   AA=EPplotXY{i};
figure(i);
    h2=axes;
    h11=scatter(AA(:,1),AA(:,2));
    axis([ExtInt(1),ExtInt(2),0,inf]);
    set(h2,'ydir','reverse');
    set(h11,'Marker','.');
    xlabel('ms');
    ylabel('trials');
    title([sprintf('Chan%d-U%d  NEV-',NEVchannels(i),NEVunits(i)),SourceFile]);
    saveas(i,sprintf('EPrawNEV_Chan%d-U%d.fig',NEVchannels(i),NEVunits(i)));
    saveas(i,sprintf('EPrawNEV_Chan%d-U%d.jpg',NEVchannels(i),NEVunits(i)));
end

if PlotEPdetection==1
  i=length(NEVunits)+1;
  AA=EPplotXY{i};
figure(i);
    h2=axes;
    h11=scatter(AA(:,1),AA(:,2));
    axis([ExtInt(1),ExtInt(2),0,inf]);
    set(h2,'ydir','reverse');
    set(h11,'Marker','.');
    xlabel('ms');
    ylabel('trials');
    title(['EPpoints NEV-',SourceFile]);
    saveas(i,['EPrawNEV_EPpoints.fig']);
    saveas(i,['EPrawNEV_EPpoints.jpg']);
end

%%Mean plots
figure(i+1);
T_vec=ExtInt(1):1000/fs:ExtInt(2);
h21=plot(T_vec,EP_SUmean);
title(['EPmean NEV-' ,SourceFile]);
xlabel('Time (ms)'); 
ylabel('count');
axis([ExtInt(1),ExtInt(2),0,0.2]);
legend(SU_Name);
set(h21(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)
saveas(i+1,'EPmeanNEVraw.fig');
saveas(i+1,'EPmeanNEVraw.jpg');

%plot binned data
figure(i+2);
T_bin=ExtInt(1):Bin:[ExtInt(2)-Bin];
h22=plot(T_bin,EP_SUmeanBinHz);
title([sprintf('EP bin%dms NEV-',Bin) ,SourceFile]);
xlabel('Time (ms)'); 
ylabel('Frequency (Hz)');
%axis(EPaxis);
legend(SU_Name);
set(h22(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)
saveas(i+2,'EPmeanNEV.fig');
saveas(i+2,'EPmeanNEV.jpg');

%% Save data
a_EP_NEV.NEV=NEV;
a_EP_NEV.a_InputParameters=a_InputParameters;
a_EP_NEV.ExtEdges=ExtEdges;
a_EP_NEV.EP_SUpoints=EP_SUpoints;
a_EP_NEV.EP_SUmean=EP_SUmean;
a_EP_NEV.EP_SUmeanBinHz=EP_SUmeanBinHz;
a_EP_NEV.T_vec=T_vec;
a_EP_NEV.T_bin=T_bin;
save('EP_NEVresults.mat','a_EP_NEV','-v7.3');