% This script experiments with sensitivity to ISI jitter with varying correlations
% in firing time.  
% Correlation is varied in the following ways: 
% 1) Epoch-based correlation is varied via the ratio of firing rate to
%    epoch rate. 
% 2) Template-based correlation is varied via the threshold. 
% 3) Epoch-based and uncorrelated trains are mixed.
% 4) Template-based and uncorrelated trains are mixed. 

clear all

jitter = [0 .001 .002 .004 .008 .016];

%correlated rate, epoch rate, coarse correlated rate, threshold,
%uncorrelated rate
cases = [ ... 
         40 40 0 0 0 70; ...
        40 60 0 0 0 70; ...
        40 100 0 0 0 70; ...
        40 300 0 0 0 70; ...       
         0 75 40 .9 0 10; ...
        0 75 40 .5 0 10; ...
        0 75 40 0 0 10; ...
        0 75 40 -1 0 10; ...
        0 75 40 -2 0 10; ...        
         0 75 40 .9 0 22; ...         
        0 75 40 .5 0 22; ...
        0 75 40 0 0 22; ...
        0 75 40 -1 0 22; ...
        0 75 40 -2 0 22; ...        
         0 75 40 .9 0 55; ...         
        0 75 40 .5 0 55; ...
        0 75 40 0 0 55; ...
        0 75 40 -1 0 55; ...
        0 75 40 -2 0 55; ...        
%          0 75 40 .9 0 30; ...
%         0 75 40 .5 0 30; ...
%         0 75 40 0 0 30; ...
%         0 75 40 -1 0 30; ...
%         0 75 40 -2 0 30; ...        
%          0 75 40 .9 0 40; ...
%         0 75 40 .5 0 40; ...
%         0 75 40 0 0 40; ...
%         0 75 40 -1 0 40; ...
%         0 75 40 -2 0 40; ...        
%          0 75 40 .9 0 50; ...
%         0 75 40 .5 0 50; ...
%         0 75 40 0 0 50; ...
%         0 75 40 -1 0 50; ...
%         0 75 40 -2 0 50; ...        
%          0 75 40 .9 0 60; ...
%         0 75 40 .5 0 60; ...
%         0 75 40 0 0 60; ...
%         0 75 40 -1 0 60; ...
%         0 75 40 -2 0 60; ...        
%          0 75 40 .9 0 70; ...
%         0 75 40 .5 0 70; ...
%         0 75 40 0 0 70; ...
%         0 75 40 -1 0 70; ...
%         0 75 40 -2 0 70; ...        
%         30 30 0 0 10 70; ...
%         20 20 0 0 20 70; ...
%         10 10 0 0 30 70; ...
%         0 75 0 0 40 70; ...
%         0 75 30 .9 10 70; ...
%         0 75 20 .9 20 70; ...
%         0 75 10 .9 30 70; ...
        ]

% row indices for various groups of above cases 
epochRateIndices = [1; 2; 3; 4];
thresholdIndices10 = [5; 6; 7; 8; 9];
thresholdIndices22 = [10; 11; 12; 13; 14];
thresholdIndices55 = [15; 16; 17; 17; 19];

x = load('signals_figure4.mat');
signals = x.signals;

dt = .0002; %time step
T = .3; 
corrT = 20; % longer simulations used for cross-correlation calculations to reduce high peaks due to random variation

meanPeakCorrelation = [];
meanAreaCorrelation = [];
meanNormAreaCorrelation = [];
meanCOV = [];
meanE = [];
sdE = [];
for i = 1:size(cases,1)
    templateFreq = 2*pi*cases(i,6);
    template = sin([dt:dt:T]*templateFreq);
    corrTemplate = sin([dt:dt:corrT]*templateFreq);    
    
    tic 
    [spikes cov] = genSynthetic(20, corrT, dt, cases(i,1), cases(i,2), .003, cases(i,3), corrTemplate, cases(i,4), cases(i,5), [1 0 0]);
    xc = ensembleCrossCorr(spikes, .001, corrT);
    meanPeakCorrelation = [meanPeakCorrelation; mean(xc(:,1)) std(xc(:,1))];
    meanAreaCorrelation = [meanAreaCorrelation; mean(xc(:,2)) std(xc(:,2))];
    meanNormAreaCorrelation = [meanNormAreaCorrelation; mean(xc(:,3)) std(xc(:,3))];
    meanCOV = [meanCOV; mean(cov) std(cov)];
    toc
    
    for j = 1:length(jitter) 
        for k = 1:size(signals,1)
            signal = signals(k,:);
            [spikes cov] = genSynthetic(500, T, dt, cases(i,1), cases(i,2), .003, cases(i,3), template, cases(i,4), cases(i,5), [1 0 0]);
            [weights, err, t] = decode(signal, dt, spikes, [jitter(j) 0 0], [0 0 0], [0 10], 0, 32, 5, 1); t
            meanE(i,j,k) = mean(err)
            sdE(i,j,k) = std(err);
        end
    end
    save 'data_correlation.mat'
end

meanPeakCorrelation
meanE