Added data files and script for recreating stimulus

ran1.cpp

@@ -0,0 +1,170 @@
+/*
+ * Copyright (c) 2008, Christian Mendl
+ * All rights reserved.
+ *
+ */
+
+// Pseudo random number generator 'ran1' used in the stimulation program.
+
+#include <mex.h>
+#include <memory.h>
+
+
+const long NTAB = 32;
+
+// seed also contains 'static' variables used in 'ran1'
+struct Seed
+{
+	long idum;
+	long iy;
+	long iv[NTAB];
+
+	Seed()
+	: idum(-1), iy(0) {
+	}
+};
+
+// function declarations
+double ran1(Seed& seed);	// random generator
+mxArray *Seed2Mat(const Seed& seed);
+void Mat2Seed(const mxArray *mat, Seed& seed);
+
+/* x = ran1(seed)
+ * x = ran1(seed,n)
+ * [x,seed] = ran1(seed)
+ * [x,seed] = ran1(seed,n)
+ * seed = ran1(seed,n,0)
+*/
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+	if (nrhs == 0 || nrhs > 3) mexErrMsgTxt("One, two or three input arguments required.");
+	if (nlhs > 2) mexErrMsgTxt("Wrong number of output arguments.");
+	if (nlhs > 1 && nrhs == 3) mexErrMsgTxt("One output argument expected given three input arguments.");
+
+	Seed seed;
+	if (mxIsDouble(prhs[0]))
+		seed.idum = (long)mxGetScalar(prhs[0]);
+	else if (mxIsInt32(prhs[0]))
+		memcpy(&seed.idum, mxGetData(prhs[0]), sizeof(long));
+	else if (mxIsStruct(prhs[0]))
+		Mat2Seed(prhs[0], seed);
+	else
+		mexErrMsgTxt("Invalid class of input argument 'seed'.");
+
+	// number of generated random numbers
+	int n;
+	if (nrhs < 2) n = 1;
+	else {
+		n = (int)mxGetScalar(prhs[1]);
+		if (n < 1) mexErrMsgTxt("'n' must be a positive integer.");
+	}
+
+	if (nrhs <= 2)
+	{
+		// generate random numbers
+		double *x = new double[n];
+		for (int i = 0; i < n; i++)
+			x[i] = ran1(seed);
+
+		// copy output: random numbers
+		plhs[0] = mxCreateDoubleMatrix(n, 1, mxREAL);
+		memcpy(mxGetPr(plhs[0]), x, n*sizeof(double));
+		delete [] x;
+		// seed
+		if (nlhs == 2) {
+			plhs[1] = Seed2Mat(seed);
+		}
+	}
+	else
+	{
+		// omit random numbers, just advance the seed value
+		for (int i = 0; i < n; i++)
+			ran1(seed);
+
+		// copy output: seed
+		plhs[0] = Seed2Mat(seed);
+	}
+}
+
+
+// convert Matlab 'mxArray' to the 'seed' structure and vice versa
+
+inline void Mat2Seed(const mxArray *mat, Seed& seed)
+{
+	// copy 'idum'
+	mxArray *idum = mxGetField(mat, 0, "idum");
+	if (!idum) mexErrMsgTxt("Structure field 'idum' not found.");
+	if (!mxIsInt32(idum)) mexErrMsgTxt("Structure field 'idum' should be of class 'int32'.");
+	memcpy(&seed.idum, mxGetData(idum), sizeof(seed.idum));
+	// copy 'iy'
+	mxArray *iy = mxGetField(mat, 0, "iy");
+	if (!iy) mexErrMsgTxt("Structure field 'iy' not found.");
+	if (!mxIsInt32(iy)) mexErrMsgTxt("Structure field 'iy' should be of class 'int32'.");
+	memcpy(&seed.iy, mxGetData(iy), sizeof(seed.iy));
+	// copy 'iv'
+	mxArray *iv = mxGetField(mat, 0, "iv");
+	if (!iv) mexErrMsgTxt("Structure field 'iv' not found.");
+	if (!mxIsInt32(iv)) mexErrMsgTxt("Structure field 'iv' should be of class 'int32'.");
+	if (!(mxGetM(iv) == NTAB && mxGetN(iv) == 1) && !(mxGetM(iv) == 1 && mxGetN(iv) == NTAB))
+		mexErrMsgTxt("Wrong dimension of structure field 'iv'.");
+	memcpy(&seed.iv, mxGetData(iv), sizeof(seed.iv));
+}
+
+inline mxArray *Seed2Mat(const Seed& seed)
+{
+	// create seed structure fields
+	mxArray *idum = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
+	memcpy(mxGetData(idum), &seed.idum, sizeof(seed.idum));
+	mxArray *iy = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
+	memcpy(mxGetData(iy), &seed.iy, sizeof(seed.iy));
+	mxArray *iv = mxCreateNumericMatrix(1, NTAB, mxINT32_CLASS, mxREAL);
+	memcpy(mxGetData(iv), &seed.iv, sizeof(seed.iv));
+
+	// create return value structure
+	const char *fieldnames[] = { "idum", "iy", "iv" };
+	mxArray *ret = mxCreateStructMatrix(1, 1, 3, fieldnames);
+	mxSetFieldByNumber(ret, 0, 0, idum);
+	mxSetFieldByNumber(ret, 0, 1, iy);
+	mxSetFieldByNumber(ret, 0, 2, iv);
+
+	return ret;
+}
+
+
+const long IA = 16807;
+const long IM = 2147483647;
+const double AM = 1.0/IM;
+const long IQ = 127773;
+const long IR = 2836;
+const long NDIV = 1+(IM-1)/NTAB;
+const double EPS = 1.2e-7;
+const double RNMX = 1.0-EPS;
+
+double ran1(Seed& seed)
+{
+	int j;
+	long k;
+	// static long iy=0;
+	// static long iv[NTAB];
+	double temp;
+
+	if (seed.idum <= 0 || !seed.iy) {
+		if (-seed.idum < 1) seed.idum=1;
+		else seed.idum = -seed.idum;
+		for (j=NTAB+7;j>=0;j--) {
+			k=seed.idum/IQ;
+			seed.idum=IA*(seed.idum-k*IQ)-IR*k;
+			if (seed.idum < 0) seed.idum += IM;
+			if (j < NTAB) seed.iv[j] = seed.idum;
+		}
+		seed.iy=seed.iv[0];
+	}
+	k=(seed.idum)/IQ;
+	seed.idum=IA*(seed.idum-k*IQ)-IR*k;
+	if (seed.idum < 0) seed.idum += IM;
+	j=seed.iy/NDIV;
+	seed.iy=seed.iv[j];
+	seed.iv[j] = seed.idum;
+	if ((temp=AM*seed.iy) > RNMX) return RNMX;
+	else return temp;
+}

recreate_spatiotemporal_white_noise_stimulus.m

@@ -0,0 +1,97 @@
+
+% 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' );