

clear all;
close all;

%% importing data
data.fid = 'ex100928_2SpOT.1_DATA3'; %fid stand for file identification of the imported data
x = importdata([data.fid '_xcorr_bursttimes2.dat'],'\t');
data.info = x.textdata(2:size(x.textdata,1),1:2);

%headers = {'data.burst.start.left','data.burst.stop.left','data.burst.peak.left','data.burst.start.right','data.burst.stop.right','data.burst.peak.right'};
burstheaders = x.textdata(1,3:size(x.textdata,2)-1);

for i = 1:length(burstheaders)
    tmp = x.data(:,i);
    eval(['data.bursttimes.' burstheaders{i} ' = tmp(~isnan(tmp));']);
end
data.xcorr.spiketimes.values = x.data(:,i+1);

clearvars -except data burstheaders

%data.bursts.peaks.left = peakburst1(~isnan(peakburst1));
%peakburst1(isnan(peakburst1)) = zeros(size(peakburst1(isnan(peakburst1))));
%ts1 = timeseries(ones(size(peakburst1)),peakburst1);
%ts2 = timeseries(ones(size(peakburst2)),peakburst2);

%% generating counting matrices (zeros and ones)

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.' 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.' burstheaders{i} '= tmp_burstcount;']);
    clearvars tmp*
end

clearvars -except data burstheaders

%     
%     burstcount(round(t*1000)/1000==eval(['data.bursttimes.' eval(burstheaders{i}) '(
%     eval(['data.counts.' burstheaders(i) ' = zeros(size(t))']);
%     
% 
% for i = 1:length(headers)
% 
%     tmp = eval(headers{i});
%     
%     for i = 1:eval(['length(' headers{i} ')'])
%         counts1(round(t*1000)/1000==peakburst1(i))=1;
%     end
% 
% peakburst2 = peakburst2(~isnan(peakburst2));
% for i = 1:length(peakburst2)
%     counts2(round(t*1000)/1000==peakburst2(i))=1;
% end

%% calculating and plotting cross-correlations

%lag = input('lag (ms): ');
lag = 2000;
%sfactor = input('smoothing factor (ms): ');
%sfactor = [1 10 50];
sfactor = [10 50];
max_xcorr1 = zeros(3,3);
max_xcorr2 = zeros(3,3);

ind_zerolag = ceil(size(data.xcorr.spiketimes.values,1)/2);
data.xcorr.spiketimes.values = data.xcorr.spiketimes.values(ind_zerolag-2000:ind_zerolag+2000,1);
data.xcorr.lags = -lag:1:lag;
[data.xcorr.spiketimes.maxval tmpind] = max(data.xcorr.spiketimes.values,[],1);
if tmpind > lag-500 && tmpind < lag+500
    data.xcorr.spiketimes.background = mean([data.xcorr.spiketimes.values(lag-1000:lag-500,1);data.xcorr.spiketimes.values(lag+500:lag+1000,1)]);
    data.xcorr.spiketimes.maxvalnorm = (data.xcorr.spiketimes.maxval - data.xcorr.spiketimes.background) / data.xcorr.spiketimes.background;
else
    data.xcorr.spiketimes.background = NaN;
    data.xcorr.spiketimes.maxvalnorm = NaN;
end
data.xcorr.spiketimes.maxlag = tmpind - lag;

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.xcorr.lags/1000,data.xcorr.spiketimes.values,'-k','LineWidth',2.5); hold on;
plot(data.xcorr.spiketimes.maxlag/1000,data.xcorr.spiketimes.maxval,'r*','MarkerSize',10);
ylim([0 data.xcorr.spiketimes.maxval*1.1])
xlabel('lag (s)'); ylabel('cross correlation')

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

        for k = 1:2
            % calculating correlation matrices for smoothing factor 1 or 2
            data.xcorr.bursttimes.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.xcorr.bursttimes.values{i,j}(:,k));']);

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

clear tmp*
y_upperlim1 = max(max(max_xcorr1));
y_upperlim2 = max(max(max_xcorr2));

x = 1; y_tick1 = 0;
while y_tick1 == 0
    y_tick1 = floor(y_upperlim1*(10^x)/4)/(10^x);
    x = x+1;
end
x = 1; y_tick2 = 0;
while y_tick2 == 0
    y_tick2 = floor(y_upperlim2*(10^x)/4)/(10^x);
    x = x+1;
end

ylim1 = [y_upperlim1*(-0.1) y_upperlim1*1.10];
ylim2 = [y_upperlim2*(-0.1) max(max(max_xcorr2))*1.10];

% plotting correlations
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);       
        [ax h1 h2] = plotyy(data.xcorr.lags/1000,data.xcorr.bursttimes.values{i,j}(:,1),data.xcorr.lags/1000,data.xcorr.bursttimes.values{i,j}(:,2));

        %         if plot_index == 1
%             ylim1 = [0-max(data.xcorr.bursttimes.values{i,j}(:,1))*0.1 max(data.xcorr.bursttimes.values{i,j}(:,1))*1.40];
%             ylim2 = [0-max(data.xcorr.bursttimes.values{i,j}(:,2))*0.1 max(data.xcorr.bursttimes.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');
        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




%data.xcorr.lag = lag;
data.xcorr.smooth = sfactor;

%% saving files

print(f1,'-dpng',[data.fid '_xcorr_spiketimes2.png'])
print(f2,'-dpng',[data.fid '_xcorr_bursttimes2.png'])
save([data.fid '_xcorrbursts2.mat'],'data');