Implementation

EnergyModel class

class EnergyModel(num_neuron=1, firing_model=[{'model': 'linear', 'params': {'r_slope': 50, 'r_bias': 60}}], peq_model={'model': 'linear_pot', 'params': {'slope': -2.65}}, p0_model=None, boundary_mode=None, D0=0.56, Nv=None, pde_solve_param={}, verbose=False)

Energy Model class.

Parameters

num_neuronint

A number of neuronal responses. The default is 1.

firing_modellist

For each neuron, this list contains the firing rate function. Each entry is either a function that returns an array of firing rate values, or a dictionary that specifies a model from firing_rate_models.py file. The default is [{“model”: “linear”, “params”: {“r_slope”: 50, “r_bias”: 60}}].

peq_modeldictionary or callable

Equilibrium probability distribution density (which defines the driving force and Langevin potential, see Supplementary Section 1.1 in Genkin et. al. 2020 paper). This quantitiy defines the Potential function via Phi(x)=-log(peq(x)). Can be specified as a function, or as a dictionary that specifies one of the models from peq_models.py file. The default is {“model”: “linear_pot”, “params”: {“slope”: -2.65}}.

p0_modeldictionary or callable or None

Initial probability distribution. Can be specified as a function, or as a dictionary that specifies one of the models from peq_models.py file. If set to None, it will be assumed to be equal to peq (the equilibirum/stationary inference). The default is None.

D0float

Diffusion, or noise magnitude value. The default is 0.56.

boundary_modestring

Specify boundary mode that will set boundary conditions for the optimization/data generation. Possbile options are “reflecting”, “absorbing”, None. If set to None, the boundary conditions can be specified in pde_solve_param dictionary. The default is None.

Nvint

A number of retained eigenvalues/eigenvectors of the operator H. Too low value can degrade the quality, while too high value can result in long computational times. The default value is N-2, where N is the total number of grid points(for SEM N=Ne*(Np-1)+1). Since there are two boundary conditions, this is the maximum possible number of the retained eigenvalues. Note, that for numberical stability it is preferable to use the default value due to possible parasitic oscillation that arise from absorption operator. This is because for the operator exp(-Hdt) the contributions from the large eigenvalues decays exponentially, while for the absorption operator A=H0 it decays linearly, and generally cannot be neglected.

pde_solve_paramdictionary

Dictionary of the parameters for PDE_Solve class, see also docstrings in PDE_Solve.py file. Optional, as every entry has the default value. Possible entries:

xbeginfloat

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

xendfloat

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

methoddictionary

A dictionary that contains 2 keys:

namestring

Specifies the method for the numerical solution of EV problem, can be either ‘FD’ or ‘SEM’ (forward differences or spectral elements method).

gridsizedictionary

Specifies number of grid size points N for ‘FD’ method, or Np and Ne for ‘SEM’ method (Ne is the number of elements, Np is the number of grid points in each element).

The default is {‘name’: ‘SEM’, ‘gridsize’: {‘Np’: 8, ‘Ne’: 256}}.

BoundConddictionary

A dictionary that specifies boundary conditions (Dirichlet, Neumann or Robin). If boundary mode is supplied, it will overwrite this variable. For this input pararmeter to be in effect, the boundary_mode should be set to None. The default is {‘leftB’: ‘Neumann’, ‘rightB’: ‘Neumann’}.

grid_mode_calcstr

Specify how to calculate SEM grid collocation points. Availiable options: ‘built_in’: using numpy.polynomial module ‘newton’: using Newton’s method to calculate zeros of Legendre polinomial for the GLL grid The default is ‘newton’.

BC_methodstr

Specify the method of boundary condition handling when transforming the EV problem into linear system of equations. Availiable options: ‘projection’: use projection operator. ‘bound_subst’: use boundary condition substitution into the first and the last equations of the associated linear system. The default is ‘projection’

int_modestr

Specify the integration mode. Availiable options:

‘full’ - use full integration matrix. ‘sparse’ - use sparse integration matrix with a bias vector.

The default is ‘full’. See Supplementary Information for M. Genkin, T. A. Engel, Nat Mach Intell 2, 674–683 (2020) for details.

verbosebool

A boolean specifying the verbose level. The true value displays more messages. The default is False.

Attributes

x_d_numpy array, dtype=float

Grid points. Will be calculated with PDE_Solve class.

w_d_numpy array, dtype=float

The corresponding SEM weights (only for SEM method).

dmat_d_numpy array, dtype=float

The matrix for numerical differentiation.

peq_numpy array, dtype=float

The equilibrium probability density evaluated on SEM grid

p0_numpy array, dtype=float

The initial distribution p0.

D_float

Noise intensity.

fr_: numpy array, dtype=float

An array with firing rate functions evaluated on SEM grid. The number of columns is equal to the number of neuronal responses.

Nint

Number of grid points

Npint

Number of grid points per element (only SEM method).

Neint

Number of elements (only SEM method).

Methods

FeatureComplexity([peq])	Calculate feature complexity from a givein peq for the case of equilibrium inference
FeatureComplexityFderiv([peq])	Variational derivative of feature complexity wrt force F for the case of equilibrium inference
SaveResults([results_type])	Save fitting results to a file
calc_F(peq)	Calculate force from the peq
calc_peq(F)	Calculate peq from the force
fit([optimizer, options])	Perform model fitting
generate_data([deltaT, time_epoch, ...])	Generate spike data and latent trajectories.
score(data[, metadata, peq, rho0, D, fr])	Evaluate negative log-likelihood for a given data and model
transform_spikes_to_isi(spikes, time_epoch)	Convert spike times to data array, which is a suitable format for optimization.

FeatureComplexity(peq=None)

Calculate feature complexity from a givein peq for the case of equilibrium inference

Parameters

peqnumpy array, dtype=float

The equilibrium probability density evaluated on SEM grid. The default is self.peq_

Returns

float

Feature complexity (only valid for equilibrium or stationary inference).

FeatureComplexityFderiv(peq=None)

Variational derivative of feature complexity wrt force F for the case of equilibrium inference

Parameters

peqnumpy array, dtype=float

The equilibrium probability density evaluated on SEM grid. The default is self.peq_

Returns

numpy array

The variational derivative of loglikelihood w.r.t. feature complexity (only valid for equilibrium or stationary inference).

SaveResults(results_type='iterations', **options)

Save fitting results to a file

Parameters

results_typestr

Type of the results to save. Supported types: ‘iterations’, which saves iteration results. The default is ‘iterations’

optionsdict

Availiable options:

resultsdictionary

dictionary with the results (mandatory)

pathstr

path to files. The default is empty str

namestr

enforce a particular file name, otherwise, the default name will be generated.

sim_idint

id of a simulation (appended to the name and saved inside the dictionary). The default is empty str.

iter_startint

Starting iteration used for automatic filename generation. The default is iter_num[0]

iter_endint

Terminal iteration used for automatic filename generation. The default is iter_num[-1]

calc_F(peq)

Calculate force from the peq

Parameters

peqnumpy array, dtype=float

Equilibirum probability distribution (peq), 1D array

Returns

Fnumpy array, dtype=float

Driving force, 1D array

calc_peq(F)

Calculate peq from the force

Parameters

Fnumpy array, dtype=float

Driving force, 1D array

Returns

resultnumpy array, dtype=float

Equilibirum probability distribution (peq), 1D array

fit(optimizer='GD', options=None)

Perform model fitting

Parameters

datanumpy array (N,2), dtype=np.ndarray.

See docstring for score function for details.

optimizerstr

Optimizier. Currently, the only availible option is ‘GD’. The default is ‘GD’.

optionsdictionary

see self._GD_optimization for availiable options of GD optimizer. The default is None

Returns

emself

A fitted EnergyModel object.

generate_data(deltaT=1e-05, time_epoch=[(0, 1)], last_event_is_spike=False)

Generate spike data and latent trajectories.

Parameters

deltaTfloat

Size of the time bin in seconds for the numerical integration of the Langevin equation. The default is 0.00001.

time_epochlist

List of N tuples, where N is the number of trials. Each tuple consists of start time and stop time in seconds. For the case of absorbing boundary, stop time will be the maximum allowed time for the trial to last (the trial will terminate before this time due to absorption, or at this time in an arbitrary latent state). The default is [(0,1)]. Example: We want to generate 100 trials that start at t=0 and end at t=2, in this case time_epoch=[(0,2)]*100

last_event_is_spikebool

If true, trial termination time will not be recorded. Otherwise, trial termination time will be recorded. The default is False.

Returns

datanumpy array (N,2), dtype=np.ndarray.

Spike data packed as numpy array of the size (N,2), where each elements is a 1D array. N is the number of trials, and for each trial the first column contains inter spike intervals (ISIs) in seconds for all neurons, and the second column contains the corresponding neuronal IDs (trial termination, if recorded, is indicated with -1). data[i][0] - 1D array, a sequence of ISIs of all neurons for the trial i. The last entry can be 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: neuron 0 spiked at times 0.12, 0.15, 0.25, and neuron 1 spiked at times 0.05, 0.2. Trial 0 started at t=0 and ended at t=0.28. In this case data[0][0]=np.array([0.05,0.07,0.03,0.05,0.05,0.03]), and data[0][1]=np.array([1,0,0,1,0,-1]).

time_binsnumpy array (N,), dtype=np.ndarray

For each trial contains times at which latent trajectory was recorded. N is the number of trials, and for each trial time is represented as 1D array of floats.

xnumpy array (N,), dtype = np.ndarray

Latent trajectories for each trial, N is the number of trials. Each entry is 1D array of floats.

metadatadictionary

A dictionary with two entries:

last_event_is_spikebool

Equals to the input parameter with the same name

absorption_eventlist (N,)

List of strings for each trial with the following entries: ‘absorbed’, if the trial terminated due to trajectory absorption, or ‘observation_ended’ if the trial terminated due to time out in an arbitrary latent state.

score(data, metadata=None, peq=None, rho0=None, D=None, fr=None)

Evaluate negative log-likelihood for a given data and model

Parameters

datanumpy array (N,2), dtype=np.ndarray.

Spike data packed as numpy array of the size (N,2), where each elements is a 1D array. N is the number of trials, and for each trial the first column contains inter spike intervals (ISIs) in seconds for all neurons, and the second column contains the corresponding neuronal IDs (trial termination, if recorded, is indicated with -1). data[i][0] - 1D array, a sequence of ISIs of all neurons for the trial i. The last entry can be 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: neuron 0 spiked at times 0.12, 0.15, 0.25, and neuron 1 spiked at times 0.05, 0.2. Trial 0 started at t=0 and ended at t=0.28. In this case data[0][0]=np.array([0.05,0.07,0.03,0.05,0.05,0.03]), and data[0][1]=np.array([1,0,0,1,0,-1]).

metadatadictionary

Metadata dictionary that supports the following options:

last_event_is_spikebool or list/array of bools

If true, trial termination time will be ignored (even if recorded). Otherwise, trial termination time will be used. The default is False. Can be specified for all trials, or for each trial separately.

absorption_eventbool or list/array of bools

If true, absorption operator will be applied in the end of the loglikelihood chain. The default is True. Can be specified for all trials, or for each trial separately.

The default is None.

peqnumpy array, dtype=float

The equilibrium probability density evaluated on SEM grid. The default is self.peq_

rho0numpy array, dtype=float

Scaled initial probaiblity distribution, normalized by sqrt(peq): rho0=p0/np.sqrt(peq). The default is self.p0_/np.sqrt(peq), if self.p0_ is not None, otherwise np.sqrt(peq).

Dfloat

Noise intensity. The default is self.D_

frnumpy array, (N,num_neuron), dtype=float

2D array that contains firing rate functions for each neuron evaluated on SEM grid. The default is self.fr_

Returns

likelihoodfloat

Negative loglikelihood of a data given a model

transform_spikes_to_isi(spikes, time_epoch, last_event_is_spike=False)

Convert spike times to data array, which is a suitable format for optimization.

Parameters

spikesnumpy array (num_neuron,N), dtype=np.ndarray

A sequence of spike times for each neuron on each trial. Each entry is 1D array of floats.

time_epochlist of tuples

List of N tuples, where N is the number of trials. Each tuple consists of the trial’s start time and end time in seconds. Note that the end time should be an actual end time, but not the timeout in the case of last_event_is_spike is True.

last_event_is_spikebool

If true, trial termination time will not be recorded. Otherwise, trial termination time will be recorded.

Returns

datanumpy array (N,2),dtype=np.ndarray.

Spike data packed as numpy array of the size (N,2), where each elements is a 1D array of floats. N is the number of trials, and for each trial the first column contains the interspike intervals (ISIs), and the second column contains the corresponding neuronal indices.

PDESolve class

class PDESolve(xbegin=-1.0, xend=1.0, method={'gridsize': {'Ne': 256, 'Np': 8}, 'name': 'SEM'}, BoundCond={'leftB': 'Neumann', 'rightB': 'Neumann'}, grid_mode_calc='newton', BC_method='projection', int_mode='full')

Numerical solution of Stourm-Liouville problem

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.

methoddictionary

A dictionary that contains 2 keys:

namestring

Specifies the method for the numerical solution of EV problem, can be either ‘FD’ or ‘SEM’ (forward differences or spectral elements method).

gridsizedictionary

Specifies number of grid size points N for ‘FD’ method, or Np and Ne for ‘SEM’ method (Ne is the number of elements, Np is the number of grid points in each element).

The default is {‘name’: ‘SEM’, ‘gridsize’: {‘Np’: 8, ‘Ne’: 256}}.

BoundConddictionary

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

grid_mode_calcstr

Specify how to calculate SEM grid collocation points. Availiable options: ‘built_in’: using numpy.polynomial module ‘newton’: using Newton’s method to calculate zeros of Legendre polinomial for the GLL grid The default is ‘newton’.

BC_methodstr

Specify the method of boundary condition handling when transforming the EV problem into linear system of equations. Availiable options: ‘projection’: use projection operator. ‘bound_subst’: use boundary condition substitution into the first and the last equations of the associated linear system. The default is ‘projection’

int_modestr

Specify the integration mode. Availiable options:

‘full’ - use full integration matrix. ‘sparse’ - use sparse integration matrix with bias.

The default is ‘full’. See Supplementary Materials 2.3 from M. Genkin, T. A. Engel, Nat Mach Intell 2, 674–683 (2020) for details.

Attributes

AD_dnumpy array (N,N), dtype=float

Integration matrix (only for SEM method).

dmat_dnumpy array (N,N), dtype=float

Differentiation matrix (only for SEM method).

dxfloat

Uniform grid step size (only for FD method).

Nint

A total number of the grid points.

Npint

Degree of each element (number of grid points in each element, only for SEM method).

Neint

A number of SEM elements (only for SEM method).

w_dnumpy array (N,), dtype=float

Weights of the nodes (on the global grid).

x_d: numpy array (N,), dtype=float

Domain grid points.

Methods

Integrate(f[, result])	Takes an indefinite integral of a function f using integration matrix (and a cumulative correction for 'sparse' int_mode).
set_BoundCond(BoundCond)	Set new boundary conditions for the Stourm-Liouville probelem
solve_EV([peq, D, q, w, mode, fr, Nv])	Solve the Sturm-Liouville eigenvalue-eigenvector problem.

Integrate(f, result=None)

Takes an indefinite integral of a function f using integration matrix (and a cumulative correction for ‘sparse’ int_mode).

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.

set_BoundCond(BoundCond)

Set new boundary conditions for the Stourm-Liouville probelem

Parameters

BoundConddictionary

Specify boundary conditions 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}} The default is {‘leftB’: ‘Neumann’, ‘rightB’: ‘Neumann, ‘leftBCoeff’: {‘c1’: 1, ‘c2’: 2} }

solve_EV(peq=None, D=1, q=None, w=None, mode='hdark', fr=None, Nv=64)

Solve the Sturm-Liouville eigenvalue-eigenvector problem. The problem can be specified either by peq, q and w functions or by the precalculated stiffmat and massmat

Parameters

peqnumpy array, dtype=float

Equilibirum probabilioty distribution that determines potential Phi(x), see Suplementary Note 1.1. 1D array.

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 and H operator

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.

Nvint

A number of eigenvectors/eigenvalues returned. The default is 64.

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).

Utility functions

FC_stationary module

Below are the complimentary functions for Model selection based on feature complexity analysis for stationary data. See Genkin, Engel, Nat. Mach. Intel. (2020) paper for details. This module is optional and should be imported separately: from neuralflow.utilities import FC_stationary.

FeatureComplexity_st(peqs, grid_params, start_index=0, stop_index=-1)

Calculate Feature complexities for a given array of peqs for stationary data.

Parameters

peqsnumpy array

2D array of peqs of size (N1,N2), where N1 - number of grid points, N2 - number of peqs.

grid_paramsdictionary

Contains two entries: ‘w’: SEM weights and ‘dmat’: SEM differentiation matrix. Can be extracted from Energymodel variable: grid_params = {‘w’: EnergyModel.w_d_, ‘dmat’: EnergyModel.dmat_d_}

start_indexint

Starting index that determines the number of points at the left boundary to be skipped for FC calculation. The default is 0.

stop_indexint

Stop index that determines the number of points at the right boundary to be skipped for FC calculation. The default is -1.

Returns

FCs: numpy array

1D array of shape (N2,) that contains FC for each of the provided peqs

FeatureConsistencyAnalysis_st(data1, data2, grid_params, KL_thres, FC_radius, FC_final, FC_stride)

Perform feature consistency analysis on two fitting results obtained from different data samples for stationary data. The algorithm works as follows:

1. Create a shared feature complexity axis.

2. For each entry in shared FC, select all models from data1 and data2 that differ from the current FC by less than FC_radius.

3. For all selected models, calculate each of the pairwise KL divergences.

4. Select and save minimal KL divergence and the corresponding indices in the original data1/data2 arrays.

Parameters

data1dictionary

Dictionary with two entries: ‘peqs’ and ‘FCs’ for fitting results on data sample 1. The first entry contains 2D array of peqs, and the second - 1D array of the corresponding feature complexities.

data2dictionary

Dictionary with two entries: ‘peqs’ and ‘FCs’ for fitting results on data sample 2. The first entry contains 2D array of peqs, and the second - 1D array of the corresponding feature complexities.

grid_paramsdictionary

Contains two entries: ‘w’: SEM weights and ‘dmat’: SEM differentiation matrix. Can be extracted from Energymodel variable: grid_params = {‘w’: EnergyModel..w_d_, ‘dmat’: EnergyModel..dmat_d_}

KL_thresfloat

Threshold KL divergence for determining optimal feature complexity.

FC_radiusfloat

Feature complexity radius that determines maximum slack in features complexities of two models.

FC_finalfloat

Maximum feature complexity on shared feature complexity axis

FC_stridefloat

Feature complexity resolution on shared feature complexity axis

Returns

FC_sharednumpy array, dtype=float

1D array of shared feature complexities (shared feature complexity axis)

KLnumpy array, dtype=float

1D array of KL divergences between the two models.

FC_opt_indint

Index of optimal KL divergence in FC_shared array

ind1_sharednumpy array, dtype=int

Indices of peqs/FCs in data1 array that correspond to each entry in FC_shared.

ind2_sharednumpy array, dtype=int

Indices of peqs/FCs in data2 array that correspond to each entry in FC_shared.

KL_divergence_st(peq1, peq2, weights, schema='pairwise')

Calculate KL divergence between two batches of distributions peq1 and peq2 for stationary data. KL = integral (p1*log(p1/p2))

Parameters

peq1: numpy array

2D array of the first batch of distributions (each column contains a distribution)

peq2: numpy array

2D array of the second batch of distributions (each column contains a distribution)

weights: numpy array

1D array of SEM weights

schema: ENUM(‘pairwise’,’all’)

‘pairwise’ mode calculates KL between each pair (peq1[:,0] peq2[:,0]), (peq1[:,1] peq2[:,1]), and so on ‘all’ mode calculates KL between all possible pairs from peq1 and peq2, e.g.

if peq1 contains 3 models and peq2 contains 4 models, 12 KL divergences will be calculated and packed in a 3x4 matrix.

Returns

numpy array

Array with KL divergences

FC_nonstationary module

Below are the complimentary functions for Model selection based on feature complexity analysis for non-stationary data. See Genkin, Hughes, Engel, Nat Commun 12, 5986 (2021) paper for details. This module is optional and should be imported separately: from neuralflow.utilities import FC_nonstationary.

FeatureComplexities(results, em, iterations)

Calculate feature complexities of the model at specified iterations, and also eigenvalues and eigenvectors of the H0 operator.

Parameters

resultsdictionary

Dictionary with the results (fitted models)

emEnergyModel

An instance of the EnergyModel class.

iterationsnumpy array

Iteration numbers on which feature complexities/EVs/EVVs should be calculated.

Returns

FCs_arraynumpy array

A 1D array of feature complexities.

lQ_arraynumpy array

A 2D array of H0 eigenvalues.

Qx_arraynumpy array

A 3D array of H0 eigenvectors.

FeatureComplexity(em, peq=None, p0=None, D=None)

Calculate feature complexity for a single model specified by (peq,D,p0).

Parameters

emEnergyModel

An instance of the EnergyModel class.

peqnumpy array

Equilibirum probability distribution.

p0numpy array

Distribution of the initial latent states.

Dfloat

Noise magnitude

Returns

FCfloat

Feature complexity.

FeatureConsistencyAnalysis(data1, data2, em, iteration_indices, FC_stride)

Calculate JS divergence and feature complexities for two sequences of models data1 and data2.

This function calculate JS(FC) curve for the case when only the model potential is optimized (and the rest of the parameters are fixed).

Parameters

data1dictionary

Dictionary with peqs fitted on datasample 1 (the 1st sequence of fitted models).

data2dictionary

Dictionary with peqs fitted on datasample 2 (the second sequence of fitted models).

emEnergyModel

An instance of energymodel class used for fitting (that contains common D, p0, and fr functions).

iteration_indicesnumpy array

Indices of peqs in data1/data2 arrays which will be used for the analsys (to accelerate the computation, one may only consider the models at a subset of iterations, rather then all availible models).

FC_strideint

This parameter determines how many models from the sequence 2 will be compared with the model from the sequence 1 at each level of feature complexity. The number of comparison at each level of feature complexity is 2*FC_stride+1.

Returns

FCnumpy array

Feature complexities.

JSnumpy array

JS divergences between the two models at each level of feature complexity.

inds1numpy array

Indices of the models at each level of feature complexity in data1 array.

inds2numpy array

Indices of the models at each level of feature complexity in data1 array.

JS_divergence(p1, p2, weights, mode='normalized')

Calculate JS divergence between two distributions

Parameters

p1numpy array

The first distribution, 1D array

p2numpy array

The second distribution, 1D array

weights: numpy array

1D array of SEM weights

modeENUM(‘normalized’,’unnormalized’)

If normalized, both distribution are assumed to be normalized (integrate to 1). Otherwise, boundary term will be accounted for.

Returns

JSfloat

JS divergence

JS_divergence_tdp(peq1, D1, p01, peq2, D2, p02, weights, lQ1, Qx1, lQ2, Qx2, terminal_time=0.5, number_of_samples=5)

Calculate JS divergence between two time-dependent probability distributions that comes from two different non-stationary Langevin dynamics.

Parameters

peq1numpy array

1D array of peq values (that defines Langevin potential) for the first dynamics.

D1float

Noise magnitude for the first dynamics.

p01numpy array

1D array of p0 distribution for the first dynamics.

peq2numpy array

1D array of peq values (that defines Langevin potential) for the second dynamics.

D2float

Noise magnitude for the second dynamics.

p02numpy array

1D array of p0 distribution for the second dynamics.

weightsnumpy array

1D array of SEM weights

lQ1numpy array

1D array of Eigenvalues of H0 operator for the first dynamics.

Qx1numpy array

2D array of Eigenfunctions of H0 operator for the first dynamics.

lQ2numpy array

1D array of Eigenvalues of H0 operator for the second dynamics.

Qx2numpy array

2D array of Eigenfunctions of H0 operator for the second dynamics.

terminal_timefloat, optional

Terminal time that defines upper limit of the time integral. The default is 0.5.

number_of_samplesint, optional

Number of time steps for numerical approxination of the time integral. The default is 5.

Returns

outfloat

JS divergence between the two time-dependent distributions.