% This shows how to load data for a single cell and recreate the binary
% spatio-temporal stimulus.
% The data is organized as as a matrix of spike-counts versus stimulus
% frames and separated into trials. A trial here is 50 sec (1500 frames)
% for running noise (i.e. non-repeated sequences) and 10 sec (300 frames)
% for frozen noise.
% To recreate the stimulus, we use a mex file that implements the ran1
% random number generator from the Numerical Recipes library. (This is what
% we use in our stimulus generator program for the experiments.)
% The code below then provides the stimulus (as a 3D matrix, see below) for
% chunks of data (several trials combined in a block), and then you have to
% insert code to use these stimulus chunks. For us, working with the
% stimulus in chunks works better than trying to get the entire stimulus
% into memory.

% by J.K. Liu and T. Gollisch; 2017/05/17

clear;
close all;

RefreshRate = 30; % Hz; stimulus update rate
pStim.xPix = 800; % pixels on the monitor
pStim.yPix = 600;
pStim.stixelwidth = 4; % monitor pixels per stimulus pixel
pStim.stixelheight = 4;
pStim.Nx = ceil(pStim.xPix/pStim.stixelwidth);
pStim.Ny = ceil(pStim.yPix/pStim.stixelheight);

% The stimulus here is separated into sections of running noise (different
% for each trial) and frozen noise (same sequence repeated).
% The numbers below give the numbers of stimulus frames in each section
pStim.RunningFrames = 1500;
pStim.FrozenFrames = 300;
pStim.seed1 = -10000;  % the random generator seed for running noise
pStim.seed2 = -20000;  % the random generator seed for frozen noise

% As a convenient way to find out the number of trials for which you need
% to recreate the stimulus, you can read in the spike count data and
% determine the number of trials as below.
% Otherwise, Ntrial can also be set by hand, e.g., to a lower number to
% start by analyzing only the first few trials, e.g. "Ntrial = 3;".
load('cell_data_01_NC.mat'); % loads one of the data sets
Ntrial = size(spk1,3);       % number of trials for non-repeated noise

% regenerate running noise by chunks
for nn = 1:Ntrial
    
    fprintf(1, 'Runing ... %d/%d \n', nn, Ntrial );
    if nn == 1
        seed = pStim.seed1;
    end;
    
    % CB1: CheckerBoard Stimulus in black and white
    % ran1: the random number generator 
    [CB1 seed] = ran1(seed, pStim.Nx*pStim.Ny*pStim.RunningFrames);
    % We generate random numbers of the stimulus through the Numerical
    % Recipes routine "ran1". We have a mex file for this that is called
    % above, which is fairly fast for generating sequences of the random
    % numbers (rather than returning just 1 random number for each call).
    % The parameters of the ran1 function are the seed and the number of
    % random numbers to be obtained.
    
    % Below, the random numbers in the inverval (0,1) are turned into
    % contrast values of -1 and 1 just like we do in the stimulation
    % program.
    CB1(CB1>=0.5) = 1;
    CB1(CB1<0.5)  = 0;
    CB1 = 2*(CB1 - 0.5);
    CB1 = reshape(CB1,pStim.Ny,pStim.Nx,[]);
    % CB1 now contains a 3D matrix with values of -1 and 1.
    % The dimensions are spatial-y versus spatial-x versus time (in
    % frames). The number of frames depends on the trials per block.
    
    % Here now is the place to do something with the recreated stimulus,
    % e.g. use it for computing the receptive field or (after having
    % determined the receptive field center and size and thereby the
    % spatial region of interest) extracting the spike-triggered stimulus
    % ensemble. Or maybe save the recreated stimulus chunk to disk for
    % later use.
    % For memory reasons, we generally don't recreate the entire
    % spatio-temporal stimulus and save it to disk.
    % E.g. use the following to save the stimulus in successive .mat files
    filename = sprintf( 'nonrepeatedStim%03d.mat', nn );
    save( filename, 'CB1' );
    
end

% For recreating the frozen noise:
[CB2] = ran1(pStim.seed2, pStim.Nx*pStim.Ny*pStim.FrozenFrames);
CB2(CB2>=0.5) = 1;
CB2(CB2<0.5)  = 0;
CB2 = 2*(CB2 - 0.5);
CB2 = reshape(CB2,pStim.Ny,pStim.Nx,[]);

% Again, probably you just want to save the stimulus to file:
save( 'frozenStim.mat', 'CB2' );