

%% importing data
clear all;
data.fid = 'ex100902_2SpOT.2_DATA01'; %fid stand for file identification of the imported data
x = importdata([data.fid '_bursttimes_info.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));

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

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 time series histograms

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];
plot_index = 1;

figure;

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}]);
        subplot(length(burstheaders)/2,length(burstheaders)/2,plot_index);
        [tmp_xcorr(:,1) tmp_lags(:,1)] = xcorr(smooth(tmp_counts1,sfactor(1)),smooth(tmp_counts2,sfactor(1)),lag,'coeff');
        [tmp_xcorr(:,2) tmp_lags(:,2)] = xcorr(smooth(tmp_counts1,sfactor(2)),smooth(tmp_counts2,sfactor(2)),lag,'coeff');
%            plot(tmp_lags,tmp_xcorr);
             %hold all;
%            [tmp_xcorrmaxval tmp_xcorrmaxlag] = max(tmp_xcorr);
%             if tmp_xcorrmaxlag > lag-500 && tmp_xcorrmaxlag < lag+500
%                 tmp_bg = mean([tmp_xcorr(lag-1000:lag-500);tmp_xcorr(lag+500:lag+1000)]);
%                 tmp_xcorrmaxnorm = (tmp_xcorrmaxval - tmp_bg)/tmp_bg;
%             else
%                 tmp_bg = NaN;
%                 tmp_xcorrmaxnorm = NaN;
%             end
%             [data.xcorr.values
%        eval(['data.xcorr.' burstheaders{i} burstheaders{j+length(burstheaders)/2} '=' 'tmp_xcorr;']);
        [ax h1 h2] = plotyy(tmp_lags(:,1)/1000,tmp_xcorr(:,1),tmp_lags(:,2)/1000,tmp_xcorr(:,2));
        %set(h1,'Color',[0.7 0.1 0.3]);
        %set(h2,'Color',[0.5 0.1 0.8]);
        %set(get(ax(1),'Ylabel'),'String',['Smoothing Factor = ',num2str(sfactor(1))]) ;
        %set(get(ax(2),'Ylabel'),'String',['Smoothing Factor = ',num2str(sfactor(2))]) ;
        if plot_index == 1
            ylim1 = [0-max(tmp_xcorr(:,1))*0.1 max(tmp_xcorr(:,1))*1.20];
            ylim2 = [0-max(tmp_xcorr(:,2))*0.1 max(tmp_xcorr(:,2))*1.20];
            legend([h1 h2],{['smoothing = ' num2str(sfactor(1)) 'ms'],['smoothing = ' num2str(sfactor(2)) 'ms']},'Location','NorthWest');
        end
        set(ax(1),'Ylim',ylim1,'Ytick',[0:0.01:ylim1(2)]);
        set(ax(2),'Ylim',ylim2,'Ytick',[0:0.01:ylim2(2)]);
        title([burstheaders{i} ' vs ' burstheaders{j+length(burstheaders)/2}]);
        xlabel('lag (s)');
        plot_index = plot_index + 1;
    end
end

% subplot(3,3,1); ylabel('startburst1'); title('startburst2'); xlabel('');
% subplot(3,3,2); title('stopburst2');
% subplot(3,3,3); title('peakburst2');
% subplot(3,3,4); ylabel('stopburst1'); xlabel('');
% subplot(3,3,7); ylabel('peakburst1'); xlabel('');
% 
%eval(['data.xcorr.lags' '=' 'tmp_lags;']);

% for i = 1:length(sfactor)
%     [data.xcorr.x{i} data.xcorr.lags{i}] = xcorr(smooth(counts1,sfactor(i)),smooth(counts2,sfactor(i)),lag,'coeff');
%     subplot(length(sfactor),1,i)
%     plot(data.xcorr.lags{i},data.xcorr.x{i},'r-');
% end


%% normalizing peak correlation

corr_maxnorm = zeros(size(sfactor));

for i = 1:length(sfactor)
    [data.xcorr.max.value(i) data.xcorr.max.lag(i)] = max(data.xcorr.x{i});
    if data.xcorr.max.lag(i) > lag-500 && data.xcorr.max.lag(i) < lag+500
        corr_bg = mean([data.xcorr.x{i}(lag-1000:lag-500);data.xcorr.x{i}(lag-500:lag+1000)]);
        data.xcorr.max.valuenorm(i) = (data.xcorr.max.value(i) - corr_bg)/corr_bg;
    else
        data.xcorr.max.valuenorm(i) = NaN;
    end
end

data.smooth.factor = sfactor;