freeman-lab/cascade

Library for fitting nonlinear neural encoding models

Cascade

Library for fitting neural encoding models

About

Imagine we measures the response of a neuron to several simultaneous stimulus and behavioral variables. This toolbox fits a family of cascade models with both linear and nonlinear stages, including NL, LN, and NLN, using maximum likelihood.

Use

The fitting routines require an input structure d, which contains at the minimum, the following fields:

t - time points per trial
k - number of trials (can be 1 if no trial structure)
u - upsampling factor for stimulus
c - number of stimuli (or behavioral variables)
samp - temporal sampling of stimulus / response ('same" or 'up')
S_ctk - matrix of stimuli / behaviors x time points x trials
R_ntk - matrix neurons x time points x trials

Fits are generated by calling:

fit = fitInit(d,'NL','loglik',20);
d = prepareStim(d,fit);
[train test] = prepareRoi(d,fit,1);
fit = fitDo(train,test,fit);

To-Do

incorporate an output nonlinearity (fitG)

cascade's Issues

Need to create d.rnd

Dear Freeman Lab,

Thank you for sharing the code! Just to see whether I am using it right, I'm (re)analyzing a bit of the original Peron-data, available at https://crcns.org/data-sets/ssc .

I was doing something like this:
`
%% Load data
load('../../an229719_2013_12_02_session.mat')
volume = 4;

trialvec = unique(s.timeSeriesArrayHash.value{volume}.trial);
Ntrial = length(trialvec);
[Nneuron, ~] = size(s.timeSeriesArrayHash.value{volume}.valueMatrix);

whiskermat = zeros(2,10000,Ntrial);
luminesmat = zeros(Nneuron,200,Ntrial);
nllmax = 0;
nwlmax = 0;
validvec = ones(size(trialvec));
for nt = 1:Ntrial
whiskernvec = find(s.timeSeriesArrayHash.value{1}.trial == trialvec(nt));
nwl = length(whiskernvec);
if nwl>nwlmax
nwlmax = length(whiskernvec);
end
if isempty(whiskernvec)
% trial not valid; no whisker data
validvec(nt)=0;
else
luminesnvec = find(s.timeSeriesArrayHash.value{volume}.trial == trialvec(nt));
nll = length(luminesnvec);
if nll>nllmax
nllmax = length(luminesnvec);
end
whiskermat(:,1:nwl,nt) = s.timeSeriesArrayHash.value{1}.valueMatrix(:,whiskernvec);
luminesmat(:,1:nll,nt) = s.timeSeriesArrayHash.value{volume}.valueMatrix(:,luminesnvec);
end
end
Ntrial = sum(validvec);
whiskermat = whiskermat(:,1:nwlmax,validvec == 1);
luminesmat = luminesmat(:,1:nllmax,validvec == 1);

%% Make input structure IPS
IPS.t = nllmax; % time points per trial
IPS.k = Ntrial; % number of trials (can be 1 if no trial structure)
IPS.u = (s.timeSeriesArrayHash.value{volume}.time(2)-s.timeSeriesArrayHash.value{volume}.time(1))/(s.timeSeriesArrayHash.value{1}.time(2)-s.timeSeriesArrayHash.value{1}.time(1));
IPS.c = 2; % number of stimuli (or behavioral variables)
IPS.samp = 'up'; % temporal sampling of stimulus / response ('same" or 'up')
IPS.S_ctk = whiskermat; % matrix of stimuli / behaviors x time points x trials
IPS.R_ntk = luminesmat; % matrix neurons x time points x trials

%% Fit cascade models

fit = fitInit(IPS,'NL','loglik',20);
IPS = prepareStim(IPS,fit);
[train test] = prepareRoi(IPS,fit,1);
fit = fitDo(train,test,fit);
`

However, I get the following error:

Reference to non-existent field 'rnd'.
Error in fitInit (line 24)
fit.rnd.n = d.rnd.n;

Apparently, my 'IPS' needs a field rnd. What is this field?

Parameters 'grad' and 'hess' are not set

In the file 'fitF', in the function 'fitF_err', parameters 'grad' and 'hess' are not set. I cannot leave them out (as they have to be defined for 'fminunc' later on in 'fitF' to work). So what should these be?