%% file information
%data.fid = 'ex100928_2SpOT.1_DATA3'; %fid stand for file identification of the imported data
%cd('D:\mingfaifong\experiments\ex100928_2SpOT.1\analyses\');

%% importing data

x = importdata([data.fid '_bursttimes.dat'],'\t');
y = importdata([data.fid '_spikecounts.dat'],'\t');

data.info = x.textdata(2:size(x.textdata,1),1:2);
burstheaders = x.textdata(1,3:size(x.textdata,2));
for i = 1:length(burstheaders)
    tmp = x.data(:,i);
    eval(['data.bursttimes.values.' burstheaders{i} ' = tmp(~isnan(tmp));']);
end
data.spikecounts.values = y.data;
data.spikecounts.headers = y.colheaders;

clearvars -except data burstheaders fidlist expdir expID triallist


%% generating counting matrices (zeros and ones) for bursts

dt = 0.001;
t = 0:dt:str2double(data.info(5,2));
c = zeros(size(t));

for i = 1:length(burstheaders)
    tmp_burstcount = c;                                             %initializing burst count array
    tmp_bursttime = eval(['data.bursttimes.values.' burstheaders{i} ';']);   %array of burst times
    for j = 1:length(tmp_bursttime)
        tmp_burstcount(round(t*1000)/1000==tmp_bursttime(j))=1;
    end
    eval(['data.burstcounts.values.' burstheaders{i} '= tmp_burstcount;']);
    clearvars tmp*
end

clearvars -except data burstheaders fidlist expdir expID triallist
%% computing cross-correlations for spikes
   % NOTE: The reference and target matrices are reversed compared to how
   % correlation is calculated in IGOR, so be careful about the signs for
   % maxlag

%lag = input('lag (ms): ');
lag = 2000;

data.spikecounts.xcorr.lags = [-lag:1:lag]';
data.spikecounts.xcorr.background = NaN;
data.spikecounts.xcorr.maxvalnorm = NaN;

data.spikecounts.xcorr.values = xcorr(data.spikecounts.values(:,1),data.spikecounts.values(:,2),lag,'coeff');
[data.spikecounts.xcorr.maxval(1) tmpind] = max(data.spikecounts.xcorr.values,[],1);
if tmpind > lag-500 && tmpind < lag+500
    data.spikecounts.xcorr.background = mean([data.spikecounts.xcorr.values(lag-1000:lag-500);data.spikecounts.xcorr.values(lag+500:lag+1000)]);
    data.spikecounts.xcorr.maxvalnorm = (data.spikecounts.xcorr.maxval - data.spikecounts.xcorr.background) / data.spikecounts.xcorr.background;
end
data.spikecounts.xcorr.maxlag = tmpind - lag;

clear tmp*

%% computing cross-correlation for bursts

%sfactor = input('smoothing factor (ms): ');
sfactor = [10 50 200];
max_xcorr1 = zeros(3,3);
max_xcorr2 = zeros(3,3);
data.burstcounts.xcorr.lags = [-lag:1:lag]';
data.burstcounts.xcorr.smoothing = sfactor;

% calculating correlation matrix
for i = 1:length(burstheaders)/2
    tmp_counts1 = eval(['data.burstcounts.values.' burstheaders{i} ';']);
    
    for j = 1:length(burstheaders)/2
        tmp_counts2 = eval(['data.burstcounts.values.' burstheaders{j+length(burstheaders)/2} ';']);

        for k = 1:length(sfactor)
            % calculating correlation matrices for smoothing factor 1 or 2
            data.burstcounts.xcorr.values{i,j}(:,k) = xcorr(smooth(tmp_counts1,sfactor(k)),smooth(tmp_counts2,sfactor(k)),lag,'coeff');
            eval(['max_xcorr' num2str(k) '(i,j) = max(data.burstcounts.xcorr.values{i,j}(:,k));']);

            % calculating max of correlation values and normalizing to background
            [data.burstcounts.xcorr.maxval{i,j}(k) tmp_xcorrmaxind(k)] = max(data.burstcounts.xcorr.values{i,j}(:,k),[],1);
            if tmp_xcorrmaxind(k) > lag-500 && tmp_xcorrmaxind(k) < lag+500
                data.burstcounts.xcorr.background{i,j}(k) = mean([data.burstcounts.xcorr.values{i,j}(lag-1000:lag-500,k);data.burstcounts.xcorr.values{i,j}(lag+500:lag+1000,k)]);
                data.burstcounts.xcorr.maxnorm{i,j}(k) = (data.burstcounts.xcorr.maxval{i,j}(k) - data.burstcounts.xcorr.background{i,j}(k))/data.burstcounts.xcorr.background{i,j}(k);
            else
                data.burstcounts.xcorr.background{i,j}(k) = NaN;
                data.burstcounts.xcorr.maxnorm{i,j}(k) = NaN;
            end
        end
        data.burstcounts.xcorr.maxlag{i,j} = tmp_xcorrmaxind - lag;
        data.burstcounts.xcorr.headers{i,j} = [burstheaders{i} '_vs_' burstheaders{j+length(burstheaders)/2}];
    end 
end

clear tmp*


%% plotting correlations for spikes

screen_size = get(0, 'ScreenSize');
f1 = figure;
set(f1, 'Position', [screen_size(3)*.2 screen_size(4)*.2 screen_size(3)*.5 screen_size(4)*.5 ] );
plot(data.spikecounts.xcorr.lags/1000,data.spikecounts.xcorr.values/data.spikecounts.xcorr.maxval,'Color',[.2 .2 .2],'LineWidth',2);
%plot(data.spikecounts.xcorr.maxlag/1000,data.spikecounts.xcorr.maxval,'r*'
%,'MarkerSize',10);
title('Normalized Cross Correlation of Spike Times Between Slices')
axis([-2 2 -0.1 1.1]);
xlabel('lag (s)'); ylabel('cross correlation')

%% plotting correlations for bursts

f2 = figure;
set(f2, 'Position', [0 0 screen_size(3)*.8 screen_size(4)*.8 ] );
plot_index = 1;

for i = 1:length(burstheaders)/2
    for j = 1:length(burstheaders)/2
        subplot(length(burstheaders)/2,length(burstheaders)/2,plot_index);
        for k = 1:length(sfactor)
            plot(data.burstcounts.xcorr.lags/1000,data.burstcounts.xcorr.values{i,j}(:,k)/data.burstcounts.xcorr.maxval{i,j}(:,k),'LineWidth',1.5);
            hold all;
        end
        set(gca,'FontSize',6,'Ylim',[-0.1 1.1]);
    %    [ax h1 h2] = plotyy(data.xcorr.lags/1000,data.burstcounts.xcorr.values{i,j}(:,1),data.xcorr.lags/1000,data.burstcounts.xcorr.values{i,j}(:,2));

        %         if plot_index == 1
%             ylim1 = [0-max(data.burstcounts.xcorr.values{i,j}(:,1))*0.1 max(data.burstcounts.xcorr.values{i,j}(:,1))*1.40];
%             ylim2 = [0-max(data.burstcounts.xcorr.values{i,j}(:,2))*0.1 max(data.burstcounts.xcorr.values{i,j}(:,2))*1.40];
%         end
        if plot_index == 3
            %legend([h1 h2],{[num2str(sfactor(1)) 'ms smoothing'],[num2str(sfactor(2)) 'ms smoothing']},'Location','NorthEast');
            legend({[num2str(sfactor(1)) 'ms smoothing'],[num2str(sfactor(2)) 'ms smoothing'],[num2str(sfactor(3)) 'ms smoothing']},'Location','NorthEast');
        end
%         set(ax(1),'FontSize',6,'Ylim',[y_upperlim1*(-0.1) y_upperlim1*1.20],'Ytick',[0:y_tick1:y_upperlim1*1.10]);
%         set(ax(2),'FontSize',6,'Ylim',[y_upperlim2*(-0.1)
%         y_upperlim2*1.20],'Ytick',[0:y_tick2:y_upperlim2*1.10]);
        title([burstheaders{i} ' vs ' burstheaders{j+length(burstheaders)/2}]);
        if plot_index > 6
            xlabel('lag (s)');
        end
        plot_index = plot_index + 1;
        
%         eval(['data.xcorr.headers' burstheaders{i} burstheaders{j+length(burstheaders)/2} '.values = tmp_xcorr;']);
%         eval(['data.xcorr.' burstheaders{i} burstheaders{j+length(burstheaders)/2} '.values = tmp_xcorr;']);
%        eval(['data.xcorr.' burstheaders{i} burstheaders{j+length(burstheaders)/2} '=' 'tmp_xcorr;']);

    end
end


%% saving files

print(f1,'-dpng',[data.fid '_xcorr_spiketimes.png'])
print(f2,'-dpng',[data.fid '_xcorr_bursttimes.png'])
save([data.fid '_xcorr.mat'],'data');