% Example code for loading and using the data used in
% "Predicting synchronous firing of large neural populations from sequential recordings"
% ArXiv: https://arxiv.org/pdf/1904.04544.pdf, 
% to appear in PLOS Computational Biology.
% These data where already published in 
% "Multiplexed computations in retinal ganglion cells of a single type"
% Deny et al., Nature communications 8 (1), 1-17, 2017
% Please cite both papers if the data have been useful for you

% Data contains multi-electrode array recordings of rat 
% retinal ganglion cells subject to visual stimulation.
% Populations of OFF Alpha cells have been identified.


%% DATA DESCRIPTION

% All datasets contain the same structures and variables.

% nPix: number of pixels for the side of a checkerboard square

% px_to_micron = linear size of pixels in microns

% Sp: cell array with the spiking times of all the cells

% Xell/Yell: X and Y coordinate for elleptical fit of spatial RFs

% dt: time bin

% df: timebin in units of sampling time (20000 * dt)

% nrep: number of repetion of the reapeted stimulus and also 
% number of unreapeted checkerboard segments

% repStart: first repetition to be considered

% rep_begin_index: frame numbers of the begin of checkerboard repetition

% rep_end_index: frame numbers of the end of checkerboard repetition

% OFF_CELLS: list of identified OFF Alpha cells

% distMat: matrix of distances between OFF cells in micron



%% EXPERIMENTS LIST


clear all
cd ~/Dropbox/dataMEA/data_copula_paper/

experiments={'2015_12_23_checkerboard' %1
    '2016_07_04_checkerboard'  %2
    '2016_07_11_checkerboard'  %3
    '2016_08_01_checkerboard'  %4
    '2017_01_18_checkerboard'  %5
    '2017_02_13_checkerboard'  %6
    '2017_02_21_checkerboard'  %7
    '2017_08_01_checkerboard'  %8
    '2014_07_03_checkerboard'  %9
    '2014_07_03_fullfield'     %10
    '2014_07_03_blackout'      %11
    '2014_07_03_twobars'       %12
    '2014_01_28_checkerboard'  %13
    '2014_01_28_fullfield'};   %14


exper = 14;    
load([ experiments{exper} '.mat']); % just an example, the others work the same apart from blackout data

N = numel(OFF_CELLS);

binrCheck=cell(N,1);
for cc = 1:N
    ST = Sp{OFF_CELLS(cc)};
    b1 = [];
    for rr=repStart:nrep
        T = Fr(rep_begin_index(rr)):df:Fr(rep_end_index(rr))+df; %time bins 16.6667
        T = T(ceil( (numel(T)+1)/2) :end);
        b1 = [b1; histcounts(ST, T)];
    end
    binrCheck{cc}= b1;
    
end

% SOME PLOTS

close all;

RFX =  Xell(:,OFF_CELLS(:));
RFY =  Yell(:,OFF_CELLS(:));

figure
subplot(1,3,1)
hold on
for cc = 1:N
    plot(RFX(:,cc)*px_to_micron+nPix, RFY(:,cc)*px_to_micron+nPix,'k')
end
plot(RFX(:,1)*px_to_micron+nPix, RFY(:,1)*px_to_micron+nPix,'b')
axis([min( RFX(:) ) max( RFX(:) ) min( RFY(:) ) max( RFY(:) ) ]*px_to_micron+nPix)
xlabel('X pos (micron)')
ylabel('Y pos (micron)')
plt1 = gca;

subplot(1,3,2)
hold on
plot((1:numel(T)-1)*dt,mean(binrCheck{3})/dt,'b')
plt2 = gca;
plt2.Position(3) = plt2.Position(3)*2+ (plt2.Position(1)- ( plt1.Position(1)+plt1.Position(3))  );
xlabel('Time (sec)')
ylabel('Firing rate (Hz)')

drawnow










