Implementation

Core classes

class model(peq, p0, D, fr, params_size, grid, peq_model=None, p0_model=None, fr_model=None, with_cuda=False)

Model class supports one or multiple Langevin dynamics defined by peq(x), p0(x), D, fr(x). For multiple dynamics, some of the parameters can be shared which is defined by params_size dictionary. For example, if we want to model 4 conditions, we might want a model with four Langevin dynamics.

Parameters:

peqnumpy array, (num_models, grid.N): peq function evaluated on a grid.
p0numpy array, (num_models, grid.N) or None: p0 function evaluated on a grid. If set to None, assume equilibirum model.
Dnumpy array, (num_models,): Diffusion, or noise magnitude.
frnumpy array, (num_models, grid.N, num_neuron): Firing rate function for each neuron.
params_sizedict: Number of distinct functions in the model. For each of the parameters (peq, p0, D, fr), it specifies whether the parameter is shared across the dynamics (in which case the value is 1), or non-shared (in which case the value is equal to num_models). For example, if there are 4 Langevin dynamics with distinct potentials, but shared p0, D, and fr, then params_size = {‘peq’: 4, ‘p0’: 1, ‘D’: 1, ‘fr’: 1}.
gridgrid.GLLgrid object: Initialized GLL grid. Should be initiatialized with with_cuda=True for GPU support.
peq_modeldict, optional: Specify peq model that was used to initialize peq. This is only needed to keep a record. The default is None.
p0_modeldict, optional: Specify p0 model that was used to initialize p0. This is only needed to keep a record. The default is None.
fr_modellist, optional: For each neuron this specifies a fr model. This is only needed to keep a record (and also can be used for data generation). The default is None.
with_cudabool, optional: Whether to include GPU support. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is False.

Methods

`FeatureComplexity`([model, model_num, ...])	Calculate feature complexity
`FeatureComplexityFderiv`(peq[, device])	Calculate variational derivative of equilirium model complexity w.r.t.
`Fr_from_fr`(fr[, device])	Calculates Fr (an auxiliary function for tuning function optimization) from fr (tuning function)
`compute_rho0`(p0, peq[, device, min_peq])	Compute rho0 from p0 and peq
`force_from_peq`(peq[, device])	Calculates driving force F from peq
`fr_from_Fr`(Fr[, C, device])	Calculates fr (tuning function) from Fr (an auxiliary function for tuning function optimization) and C (a scaling constant)
`get_fr_lambda`()	Extract firing rate function
`get_params`(model_num[, device])	Extract peq, p0, D, fr
`get_params_for_grad`(model_num[, device])	Extract peq, rho0, D, fr
`new_model`(peq_model, p0_model, D, fr_model, grid)	Initialize a model
`peq_from_force`(F[, device])	Calculates normalized peq from a given driving Force
`sync_model`([mode])	Sync GPU and CPU by coping a model over.

FeatureComplexity(model=None, model_num=0, pde_solver=None, device='CPU'): Calculate feature complexity

FeatureComplexityFderiv(peq, device='CPU'): Calculate variational derivative of equilirium model complexity w.r.t. force F. Used in Genkin, Engel 2020 for regularization.

Fr_from_fr(fr, device=None)

Calculates Fr (an auxiliary function for tuning function optimization) from fr (tuning function)

Parameters:

frnumpy array: An array that represents tuning function f(x).

Returns:

Frnumpy array: An array that represents an auxiliary function Fr(x) for tuning function optimization.

compute_rho0(p0, peq, device='CPU', min_peq=1e-20)

Compute rho0 from p0 and peq

Parameters:

p0numpy array: p0(x) distribution.
peqnumpy array: peq distribution.
devicestr, optional: “CPU” or “GPU”. The default is None, which is “CPU”.
min_peqfloat, optional: Min peq to avoid division by zero. The default is 10**-20.

Returns:

rho0numpy array: rho0 distribution.

force_from_peq(peq, device=None)

Calculates driving force F from peq

Parameters:

peqnumpy array: peq.

Returns:

Fnumpy array: Driving force.

fr_from_Fr(Fr, C=1, device=None)

Calculates fr (tuning function) from Fr (an auxiliary function for tuning function optimization) and C (a scaling constant)

Parameters:

Frnumpy array: Fr reprsented on the grid.
Cfloat, optional: Constant C. The default is 1.
devicestr, optional: “CPU” or “GPU”. The default is None, which is CPU.

Returns:

numpy array: Tuning function fr.

get_fr_lambda(): Extract firing rate function

get_params(model_num, device='CPU'): Extract peq, p0, D, fr

get_params_for_grad(model_num, device='CPU'): Extract peq, rho0, D, fr

classmethod new_model(peq_model, p0_model, D, fr_model, grid, params_size=None, with_cuda=False)

Initialize a model

Parameters:

peq_modelcallable or dict or np.array: peq function. Can be either callable lambda, or dictionary that specify peq from peq_models.py, or an array of peq values evaluated on the grid.
p0_modelcallable or dict or np.array: p0 function. Can be either callable lambda, or dictionary that specify peq from peq_models.py, or an array of peq values evaluated on the grid.
Dnp.array or float: Diffusion coefficint. Either an array that contains a value for each of the models, or a float.
fr_modellist or np.array: Firing rate functions. Either a list of callables/dict (one entry for each neuron), or np.array.
gridgrid.GLLgrid object: Initialized GLL grid. Should be initiatialized with with_cuda=True for GPU support.
params_sizedict, optional: Number of distinct functions in the model. For each of the parameters (peq, p0, D, fr), it specifies whether the parameter is shared across the dynamics (in which case the value is 1), or non-shared (in which case the value is equal to num_models). For example, if there are 4 Langevin dynamics with distinct potentials, but shared p0, D, and fr, then params_size = {‘peq’: 4, ‘p0’: 1, ‘D’: 1, ‘fr’: 1}.
with_cudabool, optional: Whether to include GPU support. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is False.

Returns:

model: Initialized model.

peq_from_force(F, device=None)

Calculates normalized peq from a given driving Force

Parameters:

Fnumpy array: Driving force.

Returns:

peqnumpy array: peq.

sync_model(mode='GPU_to_CPU'): Sync GPU and CPU by coping a model over.

class GLLgrid(xbegin=-1, xend=1, Np=8, Ne=64, with_cuda=False)

Grid class that calculates Gauss-Lobatto-Legendre grid and implements integration and differentiation on the grid.

Parameters:

xbeginfloat: The left boundary of the latent state. The default is -1.
xendfloat: The right boundary of the latent state. The default is 1.
Npint: The degree of Langrange interpolation polynomial, also the number of grid points at each element. The default is 8.
Neint: Number of SEM elements. The default is 64.
with_cudabool, optional: Whether to include GPU support. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is False.

Methods

`Differentiate`(f[, result, device])	Take a derivative of function f
`Integrate`(f[, result, device])	Indefinite integral of a function f using integration matrix.

Differentiate(f, result=None, device='CPU')

Take a derivative of function f

Parameters:

fnumpy array, dtype=float: Function values evaluated on the grid
resultnumpy array, dtype=float: A container for the results (to avoid additional allocation). If not provided, will return a result. The default is None.
devicestr, optional: ENUM(‘CPU’, ‘GPU’). The default is ‘CPU’.

Returns:

numpy array: If the result is not provided at the input, it will be returned.

Integrate(f, result=None, device='CPU')

Indefinite integral of a function f using integration matrix.

Parameters:

fnumpy array, dtype=float: Function values evaluated on the grid
resultnumpy array, dtype=float: A container for the results (to avoid additional allocation). If not provided, will return a result. The default is None.

Returns:

numpy array: If the result is not provided at the input, it will be returned.

class SpikeData(data, dformat='ISIs', time_epoch=None, num_neuron=None, with_trial_end=True, with_cuda=False)

A class for storing spike data in an ISI format and also for conversion between spiketimes and ISI formats

Parameters:

datanumpy array

Two formats are supported: 1) ISIs format: the data is numpy array of the size (N,2) of type np.ndarray, where each elements is a 1D array. N is the number of trials, and for each trial the first column contains all of the inter spike intervals (ISIs) in seconds, and the second column contains the corresponding neuronal IDs (trial termination, if recorded, is indicated with -1). Neuronal Ids start from zero. ISIs are represented as 1D arrays of floats, and neuronal indices as 1D array of integers.

data[i][0] - 1D array, ISIs of type float for the trial i.
The last entry can optionally be the time interval between the last spike and trial termination time.

data[i][1] - 1D array, neuronal IDs of type int64 for the trial
i. The last entry is -1 if the trial termination time is recorded.

Example: create a data with 3 neurons (id 0,1,2) and two trials.: Trial 0 started at time 0 s and ended at 1.55 s, where neuron 1 spiked at 0.05 s and neuron 2 spikes at 0.55 s. Trial 1 also started at 0 s, ended at 10.05 s, where neuron 0 spiked at 0.05 s, neuron 1 spiked at 3.05 s, and neuron 2 spiked at 1.05 and 6.05 s. ISIs = np.empty((2, 2), dtype=np.ndarray) ISIs[0][0] = np.array([0.05,0.5,1]) ISIs[0][1] = np.array([1,2,-1]) ISIs[1][0] = np.array([0.05,1,2,3,4]) ISIs[1][1] = np.array([0,2,1,2,-1])

2) spiketimes format, where the data is numpy array of size (num_neuron, N) of type object, N is the number of trials. Each of the entries is 1D array that specify spiketimes of each neuron on each trial. In this case time_epoch array specify each trial start and end times. For the example above, the spiketimes format would be the following:

spiketimes = np.array(

[
[np.array([], dtype=np.float64), np.array([0.05])], [np.array([0.05]), np.array([3.05])], [np.array([0.55]), np.array([1.05, 6.05])] ],

dtype=object )

timeepoch = [(0, 1.55), (0, 10.05)]

dformatstr, optional

ENUM(‘spiketimes’, ‘ISIs’). The default is ‘ISIs’.

time_epochlist, optional

For each trial, specify trial start time and trial end time as a tuple. Only needed for spiketimes format. The default is None.

num_neuronint, optional

Number of neurons in the data. If not provided, will be inferred. The default is None.

with_trial_endbool, optional

Whether trial end time is recorded. The default is True.

with_cudabool, optional

Whether to include GPU support. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is False.

Methods

`change_format`(new_format[, record_trial_end])	Convert the data between spiketimes and ISIs format.
`to_GPU`()	Copy data to GPU memory.
`transform_isis_to_spikes`(data)	Convert ISIs to spikes
`transform_spikes_to_isi`(spikes, time_epoch)	Convert spike times to ISI format, which is a suitable format for optimization.
`trial_average_fr`()	For each neuron, compute trial-average firing rate.

change_format(new_format, record_trial_end=True)

Convert the data between spiketimes and ISIs format.

Parameters:

new_formatstr: ENUM(‘spiketimes’, ‘ISIs’).
record_trial_endbool, optional: Whether to record trial end time in ISIs format. The default is True.

to_GPU(): Copy data to GPU memory.

static transform_isis_to_spikes(data): Convert ISIs to spikes

static transform_spikes_to_isi(spikes, time_epoch, record_trial_end=True): Convert spike times to ISI format, which is a suitable format for optimization.

trial_average_fr(): For each neuron, compute trial-average firing rate. This is needed to scale the initial guess of the firing rate function to match the average spike rate in the data.

Optimization routines

class Optimization(dataTR, init_model, optimizer_name, opt_options, line_search_options={}, pde_solve_params={}, boundary_mode='absorbing', save_options={}, dataCV=None, device='CPU')

Parameters:

dataTRlist: List that contains neuralflow.spike_data.SpikeData object for each datasample in the data. For each datasample a separate model will be trained, however, it is allowed to have shared parameters across datasamples (e.g. allow different potentials but the same D, p0, and fr). Usually the number of datasamples is the number of conditions in the experimental data.
init_modelneuralflow.model.model: A model object. model.num_models should be equal to number of data samples.
optimizer_namestr: ENUM(‘ADAM’, ‘GD’).
opt_optionsdict: Optimization options (see init method of each of the optimizers).
line_search_optionsdict, optional: Line search options. The default is {}.
pde_solve_paramsdict, optional: Parameters for solving FPE. The default is {}.
boundary_modestr, optional: ‘absorbing’ or ‘reflecting’. The default is ‘absorbing’.
save_optionsdict, optional: Options for saving the results. The default is {}.
dataCVlist, optional: Validation data. The default is None.
devicestr, optional: ‘CPU’ or ‘GPU’. The default is ‘CPU’.

Methods

`compute_tr_loglik`(cur_samp, epoch)	Compute training negative loglik
`compute_val_loglik`(cur_samp, epoch)	Compute validation negative loglik
`line_search_score_C`(dataTR, model, ...)	Line search function for optimization
`line_search_score_D`(dataTR, model, ...)	Line search function for optimization
`run_optimization`()	Optimize the model on the data Inputs are defined at initialization, see __init__ docstring.

compute_tr_loglik(cur_samp, epoch): Compute training negative loglik

compute_val_loglik(cur_samp, epoch): Compute validation negative loglik

line_search_score_C(dataTR, model, model_num, fr, param_num, chi): Line search function for optimization

line_search_score_D(dataTR, model, model_num, D, param_num, chi): Line search function for optimization

run_optimization()

Optimize the model on the data Inputs are defined at initialization, see __init__ docstring.

Returns:

None.

class adam_opt(dataTR, init_model, opt_options, line_search_options={}, pde_solve_params={}, boundary_mode='absorbing', save_options={}, dataCV=None, device='CPU')

ADAM optimizer

Methods

`get_dataCV`(nsample)	Return nsample sample of validation data
`get_dataTR`(nsample)	Return nsample sample of training data
`initialize`(dataTR, init_model, opt_options, ...)	Initialize ADAM optimizer

classmethod initialize(dataTR, init_model, opt_options, line_search_options, pde_solve_params={}, boundary_mode='absorbing', save_options={}, dataCV=None, device='CPU')

Initialize ADAM optimizer

Parameters:

dataTRlist

List that contains neuralflow.spike_data.SpikeData object for each datasample in the data. For each datasample a separate model will be trained, however, it is allowed to have shared parameters across datasamples (e.g. allow different potentials but the same D, p0, and fr). Usually the number of datasamples is the number of conditions in the experimental data.

init_modelneuralflow.model.model

A model object. model.num_models should be equal to number of data samples.

opt_optionsdict

Optimization options dictionary:

learning_ratedict

ADAM learning rates:: alpha : the main ADAM learning rate. beta1 : ADAM weight for momentum. The default is 0.9. beta2 : ADAM weight for RMS prop. The default is 0.99. epsilon : ADAM epsilon. The default is 10**-8.

max_epochsint

Maximum number of optimization epochs. The default is 100.

mini_batch_numberint

Number of minibatches. On each epoch ADAM will update the model by iterating through each of the minibatches. Thus, the number of iterations per each epoch is equal to the product of number of data samples and mini_batch_number. Sometimes we cannot achieve the desired number of minibatches, in which case the actual number of minibatches will be smaller than mini_batch_number.

params_to_optlist

The parameters that will be optimized. The default is [‘F’, ‘F0’, ‘D’, ‘Fr’, ‘C’].

etaffloat

Regularization strength for potenial. Was used in 2020 paper. The default is zero.

line_search_optionsdict

Line search options dictionary:

C_optdictionary

Options for line search of C:

max_fun_evalint: Maximum number of function evaluation passed to scipy.optimize.minimize. The default is 3.
epoch_schedulenumpy array, dtype = int: Epochs on which line search will be performed. The default is np.array([]).
nSearchPerEpochint: Number of line searches per epoch. The default is 1.

D_opt: dictionary.

Options for line search of D. Same entries as in C_opt.

pde_solve_paramsdict

Parameters for solveing FP equation. By default these will be inherited from model.grid. Only provide if some unusual option is needed, e.g. using some unusual boundary conditions. Possible parameters:

xbeginfloat
The left boundary of the latent state. The default is -1.

xendfloat
The right boundary of the latent state. The default is 1.

Npint
The degree of Langrange interpolation polynomial, also the number of grid points at each element. The default is 8.

Neint
Number of SEM elements. The default is 64.

BoundConddict
Boundary conditions. Only include if you want to enforce unusual boundary conditions. Otherwise, use boundary_mode to specify boundary conditions.

Nvint, optional
Number of retained eigenvalues and eigenvectors of the operator H. If set to None, will be equal to grid.N-2, which is the maximum possible value. If Dirichlet BCs are used, it is stongly recommended to set this to None to avoid spurious high-freq oscillations in the fitted functions. The default is None.

boundary_modestr

Boundary mode, can be either absorbing or reflecting. The default is ‘absorbing’.

save_optionsdict

Options for saving the results:

pathstr: Local path for saving intermediate results. The default is None, in which case nothing will be saved (the results will only be in RAM).
strideint: Save intermediate files every stride number of epochs. The default is max_iteration (only save the final file).
schedulenumpy array, dtype=int: 1D array with epoch numbers on which the fitted models (peq, D, p0, fr) are saved to RAM. The default is np.arange(1-new_sim, max_iteration+1), which includes all of the epochs.
sim_startint: Starting epoch number (0 if new simulation, otherwise int > 0). The default is 0. Don’t change as loading the intermediate results and continuing optimization is not supported yet.

dataCVlist or None

List that contains neuralflow.spike_data.SpikeData object for each datasample in the data. If not provided, validation loglik will not be calculated. The default is None.

devicestr

Can be ‘CPU’ or ‘GPU’. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is ‘CPU’.

Returns:

self: Initialized optimizer object.

class gd_opt(dataTR, init_model, opt_options, line_search_options={}, pde_solve_params={}, boundary_mode='absorbing', save_options={}, dataCV=None, device='CPU')

Gradient-descent optimizer

Methods

`get_dataCV`(nsample)	Return nsample sample of validation data
`get_dataTR`(nsample)	Return nsample sample of training data
`initialize`(dataTR, init_model, opt_options, ...)	Initialize GD optimizer

classmethod initialize(dataTR, init_model, opt_options, line_search_options, pde_solve_params={}, boundary_mode='absorbing', save_options={}, dataCV=None, device='CPU')

Initialize GD optimizer

Parameters:

dataTRlist

List that contains neuralflow.spike_data.SpikeData object for each datasample in the data. For each datasample a separate model will be trained, however, it is allowed to have shared parameters across datasamples (e.g. allow different potentials but the same D, p0, and fr). Usually the number of datasamples is the number of conditions in the experimental data.

init_modelneuralflow.model.model

A model object. model.num_models should be equal to number of data samples.

opt_optionsdict

Optimization options dictionary:

learning_ratedict: Keys are each of the parameters to be optimized, and values are the learning rates, e.g. {‘F’: 0.01, ‘F0’: 0.02}.
max_epochsint: Maximum number of optimization epochs. The default is 100.
mini_batch_numberint: Number of minibatches. On each epoch ADAM will update the model by iterating through each of the minibatches. Thus, the number of iterations per each epoch is equal to the product of number of data samples and mini_batch_number.
params_to_optlist: The parameters that will be optimized. The default is [‘F’, ‘F0’, ‘D’, ‘Fr’, ‘C’].
etaffloat: Regularization strength for potenial. The default is zero.

line_search_optionsdict

Line search options dictionary:

C_optdictionary

Options for line search of C:

max_fun_evalint: Maximum number of function evaluation passed to scipy.optimize.minimize. The default is 3.
epoch_schedulenumpy array, dtype = int: Epochs on which line search will be performed. The default is np.array([]).
nSearchPerEpochint: Number of line searches per epoch. The default is 1.

D_opt: dictionary.

Options for line search of D. Same entries as in C_opt.

pde_solve_paramsdict

Parameters for solveing FP equation. By default these will be inherited from model.grid. Only provide if some unusual option is needed, e.g. using some unusual boundary conditions. Possible parameters:

xbeginfloat
The left boundary of the latent state. The default is -1.

xendfloat
The right boundary of the latent state. The default is 1.

Npint
The degree of Langrange interpolation polynomial, also the number of grid points at each element. The default is 8.

Neint
Number of SEM elements. The default is 64.

BoundConddict
Boundary conditions. Only include if you want to enforce unusual boundary conditions. Otherwise, use boundary_mode to specify boundary conditions.

Nvint, optional
Number of retained eigenvalues and eigenvectors of the operator H. If set to None, will be equal to grid.N-2, which is the maximum possible value. If Dirichlet BCs are used, it is stongly recommended to set this to None to avoid spurious high-freq oscillations in the fitted functions. The default is None.

boundary_modestr

Boundary mode, can be either absorbing or reflecting. The default is ‘absorbing’.

save_optionsdict

Options for saving the results:

pathstr: Local path for saving intermediate results. The default is None, in which case nothing will be saved (the results will only be in RAM).
strideint: Save intermediate files every stride number of epochs. The default is max_iteration (only save the final file).
schedulenumpy array, dtype=int: 1D array with epoch numbers on which the fitted models (peq, D, p0, fr) are saved to RAM. The default is np.arange(1-new_sim, max_iteration+1), which includes all of the epochs.
sim_startint: Starting epoch number (0 if new simulation, otherwise int > 0). The default is 0. Don’t change as loading the intermediate results and continuing optimization is not supported yet.

dataCVlist or None

List that contains neuralflow.spike_data.SpikeData object for each datasample in the data. If not provided, validation loglik will not be calculated. The default is None.

devicestr

Can be ‘CPU’ or ‘GPU’. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is ‘CPU’.

class Grads(pde_solve_params, boundary_mode='absorbing', grad_list=[], num_neuron=1, with_trial_end=True, device='CPU')

Grad class for loglikelihoods and gradient calculations.

Parameters:

pde_solve_paramsdict

Parameters for solveing FP equation:

xbeginfloat: The left boundary of the latent state. The default is -1.
xendfloat: The right boundary of the latent state. The default is 1.
Npint: The degree of Langrange interpolation polynomial, also the number of grid points at each element. The default is 8.
Neint: Number of SEM elements. The default is 64.
BoundConddict: Boundary conditions. Only include if you want to enforce unusual boundary conditions. Otherwise, use boundary_mode to specify boundary conditions.

boundary_modeENUM(‘absorbing’, ‘reflecting’), optional

Boundary mode. The default is ‘absorbing’.

grad_listlist, optional

List of paramters for gradient calculation. Can include ‘F’, ‘F0’, ‘D’, ‘C’, ‘Fr’. The default is [].

num_neuronint, optional

Number of neurons in the model/data. The default is 1.

with_trial_endbool, optional

Whether to take into account trial end time. The default is True.

deviceENUM(‘CPU’, ‘GPU’), optional

Device for the computations. The default is ‘CPU’. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is ‘CPU’.

Methods

get_grad_data(data, model[, model_num, ...])

Calculates loglikelihood and/or gradients

get_grad_data(data, model, model_num=0, mode='loglik', EV_solution=None)

Calculates loglikelihood and/or gradients

Parameters:

datadatanumpy array (num_trials, 2), dtype = np.ndarray.: Data in ISI format. See spike_data class for details.
model: model: An instance of neuralflow.model.
model_numint, optional: Model number to be used. The default is 0.
modeENUM(‘loglik’, ‘gradient’): Whether to only compute loglik, or also a gradient.
EV_solutiondicionary: Dictionary with the solution of the eigenvector-eigenvalue (EV) problem. The format is {‘lQ’: lQ, ‘QxOrig’: QxOrig, ‘Qx’: Qx, ‘lQd’:lQd, ‘Qd’: Qd}. If not provided, will be calculated. The default is None.

Returns:

results_alldictionary: Dictionary with the results. Possible entries depends on the grad_list, and can include ‘loglik’, ‘F’, ‘F0’, ‘D’, ‘Fr’, ‘C’.

class PDESolve(xbegin=-1.0, xend=1.0, Np=8, Ne=64, BoundCond={'leftB': 'Dirichlet', 'rightB': 'Dirichlet'}, Nv=None, with_cuda=False)

Numerical solution of Stourm-Liouville problem and Fokker-Planck.

Parameters:

xbeginfloat, optional

The left boundary of the latent state. The default is -1.

xendfloat: The right boundary of the latent state. The default is 1.
Npint: The degree of Langrange interpolation polynomial, also the number of grid points at each element. The default is 8.
Neint: Number of SEM elements. The default is 64.

BoundConddict

A dictionary that specifies boundary conditions (Dirichlet, Neumann or Robin). The default is {‘leftB’: ‘Dirichlet’,

‘rightB’: ‘Dirichlet’}.

In case of Robin, the dictionary must also contain ‘leftBCoeff’ and/or ‘rightBCoeff’, each are dictionaries with two entries ‘c1’ and ‘c2’ that specify BCs coefficients in the following format: c1*y[xbegin]+c2*y’[xbegin]=0 (left boundary) c1*y[xend]+c2*y’[xend]=0 (right boundary)

Nvint, optional: Number of retained eigenvalues and eigenvectors of the operator H. If set to None, will be equal to grid.N-2, which is the maximum possible value. If Dirichlet BCs are used, it is stongly recommended to set this to None to avoid spurious high-freq oscillations in the fitted functions. The default is None.
with_cudabool, optional: Whether to include GPU support. For GPU optimization, the platform has to be cuda-enabled, and cupy package has to be installed. The default is False.

Methods

`init_from_grid`(grid[, boundary_mode, BoundCond])	Init from grid object
`set_BoundCond`(BoundCond)	Set new boundary conditions for the Stourm-Liouville probelem
`solve_EV`([peq, D, q, w, mode, fr, device])	Solve the Sturm-Liouville/FP/FPE eigenvalue-eigenvector problem.

classmethod init_from_grid(grid, boundary_mode=None, BoundCond=None)

Init from grid object

Parameters:

gridneuralflow.grid

Grid object.

boundary_modestr, optional

Absorbing or reflecting. If None, BoundCond will be used for boundary conditions. The default is None.

BoundConddict

A dictionary that specifies boundary conditions (Dirichlet, Neumann or Robin). The default is {‘leftB’: ‘Dirichlet’,

‘rightB’: ‘Dirichlet’}.

In case of Robin, the dictionary must also contain ‘leftBCoeff’ and/or ‘rightBCoeff’, each are dictionaries with two entries ‘c1’ and ‘c2’ that specify BCs coefficients in the following format: c1*y[xbegin]+c2*y’[xbegin]=0 (left boundary) c1*y[xend]+c2*y’[xend]=0 (right boundary)

Returns:

self: Initialized PDESolve object.

set_BoundCond(BoundCond)

Set new boundary conditions for the Stourm-Liouville probelem

Parameters:

BoundConddictionary or str: If str, can be ‘absorbing’ or ‘reflecting’. Alternatively, specify boundary conditions as a dict: keys : ‘leftB’, ‘rightB’, (optionally: ‘leftBCoeff’, ‘rightBCoeff’) values : ‘Dirichlet’ ‘Neumann’ or ‘Robin’. If ‘Robin’, addionally specify coefficients as a dictionary with two keys: [c1,c2], consistent with the boundary condition of the form: c1*y(B)+c2*y’(B)=0. Example: {‘leftB’:’Robin’,’leftBCoeff’:{‘c1’=1, ‘c2’=2},

‘rightB’:’Robin’,’rightBCoeff’:{‘c1’=3, ‘c2’=4}}

solve_EV(peq=None, D=1, q=None, w=None, mode='hdark', fr=None, device='CPU')

Solve the Sturm-Liouville/FP/FPE eigenvalue-eigenvector problem.

Parameters:

peqnumpy array, dtype=float

Equilibirum probabilioty distribution that determines potential Phi(x), or a function p(x) in S-L problem.

Dfloat

Noise magnitude.

qnumpy array, dtype=float

A function q(x) in the S-L problem. The default value is None, in this case q(x)=0.

wnumpy array, dtype=float

A function w(x) in the S-L problem (non-negative). The default is None, in this case w(x)=1.

modestr

Specify mode. Availiable modes:: ‘normal’: solve Sturm-Liouville problem, ignore D and fr. ‘h0’: solve for eigenvalues and vectors of FP operator H0. ‘hdark’: solve for eigenvalues and vector of FP operator H.

The default is ‘hdark’.

frnumpy array

The firing rate function (required for ‘hdark’ mode). This firing rate function is an elementwise sum of the firing rate functions of all the neuronal responses. The default is None.

devicestr

Can be ‘CPU’ or ‘GPU’.

Returns:

lQnumpy array (Nv,), dtype=float: The least Nv eigenvalues for the eigenvalue problem of H0 operator.
QxOrignumpy array (Nv,Nv), dtype=float: The corresponding scaled eigenvectors
Qxnumpy array (Nv,Nv), dtype=float: The eigenvectors of EV problem of H0 operator (only for ‘h0’ and ‘hdark’ modes).
lQd: numpy array (Nv,), dtype=float: The eigenvalues of H operator (only for ‘hdark’ mode).
Qd: numpy array (Nv,Nv), dtype=float: The corresponding eigenvectors in H0 basis (only for ‘hdark’ mode).

Predefined peq models

Template peq functions. Equilibirum probability distribution is related to potential U(x)=-log(peq).

cos_fourth_power(x, w, xbegin=-1.0, xend=1.0)

Cosine^4 peq model.

y, L –> peq(y, L) ~ cos( y*pi/L )**2, integral peq(x) dx =1: where y - x centered on the middle of the domain, L - domain length

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
xbegin: float: Left boundary of the x domain. The default is -1.
xend: float: Right boundary of the x-domain. The default is 1.

Returns:

peqnumpy array, dtype=float: Probability density distribution evaluated at grid ponits x.

cos_square(x, w, xbegin=-1.0, xend=1.0)

Cosine-squareed model.

y, L –> peq(y, L) ~ cos( y*pi/L )**2

integral peq(x) dx =1 y - x centered on the middle of the domain L - domain length

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
xbegin: float: Left boundary of the x domain. The default is -1.
xend: float: Right boundary of the x-domain. The default is 1.

Returns:

peqnumpy array (N,), dtype=float: Probability density distribution evaluated at grid ponits x.

custom(x, w, lambdafunc=None)

Custom model (for example, for real data fitting).

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
lambdafunc: function or None: Either supply a function, or leave None if the model of the generated data is not known, in which case x and w are ignored.

Returns:

peqnumpy array (N,), dtype=float: Probability density distribution evaluated at grid ponits x.

double_well(x, w, xmax=0, xmin=0.3, depth=4.0, sig=0)

Peq derived from a double-well potential.

peq(x) = exp[ -U(x) ]

y = x - xmax ymin = xmin - xmax

U(y) = depth * (y/ymin)**4 - 2*depth * (y/ymin)**2 +sig*y

integral peq(x) dx =1

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
xmaxfloat: position of the maximum of the potential in the symmetrix case (sig=0)
xminfloat: position of the two symmetric minima for sig=0 case. The default is 0.
depthfloat: depth of the potential wells (U_max - U_min) for the symmetric case (sig=0). The default is 4.
sigfloat: assymetry parameter, to bias one well over another. The default is 0.

Returns:

peqnumpy array (N,), dtype=float: Probability density distribution evaluated at grid ponits x.

linear_pot(x, w, slope=1)

Ramping model derived from a linear potential

V(x)=slope*x peq(x) = exp (-slope*x) / || exp (-slope*x) ||

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
slope: float: Slope of the potential function.

Returns:

peqnumpy array, dtype=float: Probability density distribution evaluated at grid ponits x.

manual(x, w, interp_x=[-1, -0.3, 0, 0.3, 1], interp_y=[0, 10, 5, 10, 0], bc_left=[(1, 0), (1, 0)], bc_right=[(1, 0), (1, 0)], stationary_points=None): Define arbitrary model by specifiying minima and maxima of the potential

mixture(x, w, theta=0.5, model1={'model': 'cos_square', 'params': {}}, model2={'model': 'double_well', 'params': {}})

Convex mixture of two models.

peq = theta*peq_model_1 + (1-theta)*peq*model_2

integral peq(x) dx =1

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
thetafloat: convex weight, 0<=theta<=1. The default is 0.5.
model1dictionary.: Dictionary that contains model name (one of the built-in models) and the corresponding parameters. The default is {‘model’: ‘cos_square’, ‘params’: {}}
model2dictionary.: Dictionary that contains model name (one of the built-in models) and the corresponding parameters. The default is {‘model’: ‘double_well’, ‘params’: {}}

Returns:

peqnumpy array (N,), dtype=float: Probability density distribution evaluated at grid ponits x.

single_well(x, w, xmin=0, miu=10.0, sig=0)

Peq derived from a single well quadratic potential.

peq(x) = exp[ -U(x) ]

y = x - xmin

U(y) = miu * y**2 + sig*y

integral peq(x) dx =1

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
xminfloat: position of the minimum for sig=0 case. The default is 0.
miufloat: curvature (steepness) of the potential. The default is 10.
sigfloat: assymetry parameter, to bias left vs. rights sides of the potential. The default is 0.

Returns:

peqnumpy array, dtype=float: Probability density distribution evaluated at grid ponits x.

stepping(x, w, interp_x=[-1, -0.5, 0, 0.3, 1], interp_y=[0, 4, 0, 1, -1], bc_left=[(1, 0), (1, 0)], bc_right=[(1, 0), (1, 0)])

Stepping model.

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
interp_x: numpy array (5,), dtype=float, or list: The x positions of the boundaries, maixma and minima of the potential, sorted from the left to the right. Contains 5 values: the left boundary, the first maximum, the minimum, the second maximum and the right boundary. The default is [-1, -0.5, 0, 0.3, 1], which corresponds to the stepping potential used in M. Genkin, O. Hughes, T.A. Engel, Nat Commun 12, 5986 (2021).
interp_y: numpy array (5,), dtype=float, or list: The corresponding y-values of the potential at the points specified by interp_x. The default is [0, 4, 0, 1, -1], which corresponds to the stepping potential used in M. Genkin, O. Hughes, T.A. Engel paper, Nat Commun 12, 5986 (2021).
bc_left: list: A list that contains two tuples that specify boundary conditions for the potential on the left boundary. The format is the same as in bc_type argument of scipy.interpolate.CubicSpline function. The default is [(1, 0), (1, 0)], which corresponds to a zero-derivative (Neumann) BCs for the potential.
bc_right: list: Same as bc_left, but for a right boundary. The default is [(1, 0), (1, 0)], which corresponds to a zero-derivative (Neumann) BCs for the potential.

Returns:

peqnumpy array, dtype=float: Probability density distribution evaluated at grid ponits x.

uniform(x, w)

Uniform peq model, derived from a constant potential.

peq = const

integral peq(x) dx =1

Parameters:

xnumpy array (N,), dtype=float: Grid points in which the model will be evaluated. N is the number of grid points.
wnumpy array (N,), dtype=float: Weights used to evaluate integrals by the Gaussian quadrature.
Returns
——-
peqnumpy array (N,), dtype=float: Probability density distribution evaluated at grid ponits x.

Predefined firing rate models

Template firing rate functions.

cos_square(x, amp=1000, freq=1.5707963267948966, bias=0)

Cosine-squareed model with a given amplitude and frequency

Parameters:

xnumpy array, dtype=float: SEM grid points
ampfloat: Amplitude. The default is 1000.
freqfloat: Frequency. The default is np.pi/2.
biasfloat: Bias. The default is 0.

Returns:

numpy array: Firing rate function evaluated on SEM grid

custom(x, lambdafunc=None)

Custom fr model. Either supply a function, or leave None if the model of the generated data is not known (and define fr manually after).

Parameters:

xnumpy array, dtype=float: SEM grid points
lambdafuncfunction object.: The default is None.

Returns:

frnumpy array: Firing rate function evaluated on SEM grid

linear(x, slope=50.0, bias=2)

Linear firing rate model. r(x, slope, bias) = max[slope * x + bias, 0]

Parameters:

xnumpy array, dtype=float: SEM grid points
slopefloat: Firing-rate slope parameter. The default is 50.0.
biasfloat: Firing threshold parameter. The default is 2.

Returns:

numpy array: Firing rate function evaluated on SEM grid

peaks(A*exp((x-x0)^2/w^2))

Sum of gaussian peaks f= SUM(A*exp((x-x0)^2/w^2))

Parameters:

xnumpy array, dtype=float: SEM grid points
centernumpy array, dtype=float: Centers of Gaussian peaks. The default is np.array([-0.5,0.5]).
widthnumpy array, dtype=float: Widths of Gaussian peaks. The default is np.array([0.2,0.3]).
ampnumpy array, dtype=float: Magnitudes of Gaussian peaks. The default is np.array([100,80]).

Returns:

numpy array: Firing rate function evaluated on SEM grid

rectified_linear(x, slope=50.0, x_thresh=-1.0)

Rectified-linear firing rate model. r(x, slope, thresh) = max[slope*(x - x_thresh), 0]

Parameters:

xnumpy array, dtype=float: SEM grid points
slopefloat: Firing-rate slope parameter. The default is 50.0.
thresholdfloat: Firing threshold parameter. The default is -1.0.

Returns:

numpy array: Firing rate function evaluated on SEM grid

sinus(x, bias=100, amp=30, freq=3.141592653589793)

Rectified sinusoid

Parameters:

xnumpy array, dtype=float: SEM grid points
biasfloat: Bias. The default is 100.
ampfloat: Amplitude of the sinusoid. The default is 30.
freqfloat: Frequency of the sinusoid. The default is np.pi.

Returns:

numpy array: Firing rate function evaluated on SEM grid

Feature complexity analysis

class FC_tools(non_equilibrium, model, terminal_time=1, number_of_timesteps=10, boundary_mode='absorbing')

Feature consistency analysis for model selection

Parameters:

non_equilibriumbool: Whether or not the model is non_equilibrium.
modelneuralflow.model: Model class variable, which is used to extract those parameters that were not optimized. Also used for the integration/differention operations on GLL grid.
terminal_timefloat, optional: Terminal time for non-equilibrium JS divergence calculation in seconds. Only used for non-equilibrium analysis. Should set to a a typical reaction time for absorbing case. If working with reflective mode, set it to trial duration. Different trial durations for reflecting are not supported. The default is 1.
number_of_timestepsint, optional: Number of timesteps for non-equilibrium JS divergence calculation. Only used for non-equilibrium analysis. It is recommended to set this equal to ten times terminal_time. The default is 10.
boundary_modestr, optional: Boundary mode is needed for non-equilibrium analysis, as we need to solve FPE to compute feature complexity and JS divergence. The default is ‘absorbing’.

Methods

`FeatureConsistencyAnalysis`(data1, data2[, ...])	Perform feature consistency analysis
`JS_divergence_non_st`(peq1, peq2)	JS divergence between two distributions peq1 and peq2 that are not normalized.
`JS_divergence_st`(peq1, peq2[, normalized])	JS divergence between two stationary distributions peq1 and peq2.
`JS_divergence_timedep`(peq1, p01, lQ1, Qx1, ...)	Compute JS divergence between two non-stationary Langevin models
`NeedToReflect`(result1, result2)	Find if model2 has to be reflected around the axis y=0.
`compute_FC`([data, indices, invert])	Compute feature complexities for selected models in the data.
`moving_average`(x, w)	Moving average filter

FeatureConsistencyAnalysis(data1, data2, JS_thres=0.0015, FC_stride=5, smoothing_kernel=0, num_models=10000, epoch_offset=0)

Perform feature consistency analysis

Parameters:

data1Dictionary: Dictionary with the fitting results on datasample 1.
data2Dictionary: Dictionary with the fitting results on datasample 2.
JS_thresfloat, optional: Threshold JS at which models are considered different. The default is 0.0015.
FC_strideint, optional: Slack in the second sequence of models. For each Feature complexity calculated from data1, we consider 2*FC_stride + 1 models from data2 and choose the one that minimizes JS. The default is 5.
smoothing_kernelint, optional: Smooth JS with moving average before thresholding. The default is 0.
num_modelsint, optional: Number of model on which feature complexities will be computed. The code will pick the models from data1/data2 on logarithmic scale, so that more models on early epochs will be selected. Smaller number can accelerate the code, but can be inaccurate. The default is 10000.
epoch_offsetint, optional.: Ignore threshold breaks for the first epochs_offset epochs. The optimization speed is usually faster initially, so JS might diverge because of that. The default is 0.

Returns:

FCs1numpy array: Feature complexities for data1. Also serves as a “common” feature complexity axis.
min_inds_1numpy array: Indices of models in data1 for each FC1. To get peq of a model that corresponds to FC1[i], we take data1[‘peq’][min_inds_1[i]].
FCs2numpy array: Feature complexities for data2.
min_inds_2numpy array: Indices of models in data2 for each FC2. To get peq of a model that corresponds to FC2[i], we take data2[‘peq’][min_inds_2[i]].
JSnumpy array: Jensen-Shannon divergence between models from data1 and data2.
FC_opt_indint: The index of an optimal model in FC1/JS.

JS_divergence_non_st(peq1, peq2): JS divergence between two distributions peq1 and peq2 that are not normalized.

JS_divergence_st(peq1, peq2, normalized=False): JS divergence between two stationary distributions peq1 and peq2. Note that in 2020 paper a slightly different metric was used: we used symmetrized KL divergence, which is not the same as JS divergence

JS_divergence_timedep(peq1, p01, lQ1, Qx1, peq2, p02, lQ2, Qx2): Compute JS divergence between two non-stationary Langevin models

NeedToReflect(result1, result2): Find if model2 has to be reflected around the axis y=0.

compute_FC(data={}, indices=None, invert=False)

Compute feature complexities for selected models in the data.

Parameters:

datadictionary, optional: Dictionary with the fitting results. The default is {}.
indicesarray, optional: Indices of models in data dictionary on which FCs will be computed. If None, will compute FC on all indices. The default is None.
invertbool, optional: Whether to invert potenial and p0. The default is False.

Returns:

FCs_array: numpy array: Evaluated feature complexities.
lQ_array: numpy array: Eigenvalues of H0 operator needed for further analysis. Only for non-equilibirum mode.
Qx_array: numpy array: Eigenvectors of H0 operator needed for further analysis. Only for non-equilibirum mode.

static moving_average(x, w)

Moving average filter

Parameters:

xnumpy array: signal.
wint: Kernel width.

Returns:

outnumpy array: Filtered array.

Data generation from Langevin dynamics

class SyntheticData(model, boundary_mode, dt=0.0001, record_trial_end=True)

Spike data generation from a Langevin model

Parameters:

modelmodel: An instance of neuralflow.model.
boundary_modeENUM(‘absorbing’, ‘reflecting’): Boundary behavior.
dtfloat, optional: Time bin size for Langevin ODE integration. The default is 0.0001.
record_trial_endbool, optional: Whether or not to include trial end. Usually True is the best choice. The default is True.

Methods

generate_data([trial_start, trial_end, ...])

Generate spike data and latent trajectories.

Notes

No GPU support. Data generation is usually pretty fast, so CPU is sufficient for this purpose.

generate_data(trial_start=0, trial_end=1, num_trials=2, model_num=0)

Generate spike data and latent trajectories.

Parameters:

trial_startfloat or list, optional: Trial start time. To specify the same trial start time for all trials, provide a single float, or provide a list for each trial. The default is 0.
trial_endfloat or list, optional: Trial end time. To specify the same trial end time for all trials, provide a single float, or provide a list for each trial. Note that for absorbing boundary mode this serves as a timeout and the actual trial end time can be smaller if the boundary is reached before this time. The default is 1.
num_trialsint, optional: Number of trials to be generated. The default is 2.
model_numint, optional: Which model to use for data generation. The default is 0.

Returns:

datanumpy array (num_trials, 2), dtype = np.ndarray.: Data in ISI format. See spike_data class for details.
time_binsnumpy array (num_trials,), dtype = np.ndarray: For each trial contains times at which latent trajectory was recorded. Each entry is 1D numpy array of floats.
xnumpy array (num_trials,), dtype = np.ndarray: Latent trajectories for each trial. Each entry is 1D array of floats. Same shape as time_bins

Other modules and classes

Viterbi algorithm. Only supported on CPU

class Viterbi(grad)

Viterbi algorithm

Parameters:

gradneuralflow.Grads: Gradient object is needed for solving FPE.

Methods

run_viterbi(data, model[, model_num])

Viterbi algorithm Finds the latent path X that maximizes joint probability P(X,Y).

Returns:

None.

run_viterbi(data, model, model_num=0)

Viterbi algorithm Finds the latent path X that maximizes joint probability P(X,Y). The path X is sampled at trial start/trial end time, and at the time of each of spike.

Parameters:

dataneuralflow.spike_data.SpikeData.: SpikeData object in ISI format.
modelneuralflow.model.model: A model object.
model_nummodel number to be used. The default is 0.

Returns:

trajectoriesnumpy array, dtype=np.ndarray: Latent trajectory sampled at spike times for each trial in the data.
state_indsindeces of trajectory latent states on the gird.