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_neuron : int

A number of neuronal responses. The default is 1.

firing_model : list

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_model : dictionary 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_model : dictionary 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.

D0 : float

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

boundary_mode : string

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.

Nv : int

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_param : dictionary

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:

xbegin : float

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

xend : float

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

method : dictionary

A dictionary that contains 2 keys:

name : string

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

gridsize : dictionary

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

BoundCond : dictionary

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_calc : str

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_method : str

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_mode : str

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.

verbose : bool

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.

N : int

Number of grid points

Np : int

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

Ne : int

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:

peq : numpy 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:

peq : numpy 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_type : str

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

options : dict

Availiable options:

results : dictionary

dictionary with the results (mandatory)

path : str

path to files. The default is empty str

name : str

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

sim_id : int

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

iter_start : int

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

iter_end : int

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

calc_F(peq)

Calculate force from the peq

Parameters:

peq : numpy array, dtype=float

Equilibirum probability distribution (peq), 1D array

Returns:

F : numpy array, dtype=float

Driving force, 1D array

calc_peq(F)

Calculate peq from the force

Parameters:

F : numpy array, dtype=float

Driving force, 1D array

Returns:

result : numpy array, dtype=float

Equilibirum probability distribution (peq), 1D array

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

Perform model fitting

Parameters:

data : numpy array (N,2), dtype=np.ndarray.

See docstring for score function for details.

optimizer : str

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

options : dictionary

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

Returns:

em : self

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:

deltaT : float

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

time_epoch : list

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_spike : bool

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

Returns:

data : numpy 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_bins : numpy 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.

x : numpy array (N,), dtype = np.ndarray

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

metadata : dictionary

A dictionary with two entries:

last_event_is_spike : bool

Equals to the input parameter with the same name

absorption_event : list (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:

data : numpy 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]).

metadata : dictionary

Metadata dictionary that supports the following options:

last_event_is_spike : bool 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_event : bool 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.

peq : numpy array, dtype=float

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

rho0 : numpy 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).

D : float

Noise intensity. The default is self.D_

fr : numpy 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:

likelihood : float

Negative loglikelihood of a data given a model

transform_spikes_to_isi(spikes, time_epoch, last_event_is_spike)

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

Parameters:

spikes : numpy 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_epoch : list 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_spike : bool

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

Returns:

data : numpy 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:

xbegin : float

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

xend : float

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

method : dictionary

A dictionary that contains 2 keys:

name : string

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

gridsize : dictionary

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

BoundCond : dictionary

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

grid_mode_calc : str

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_method : str

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_mode : str

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_d : numpy array (N,N), dtype=float

Integration matrix (only for SEM method).

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

Differentiation matrix (only for SEM method).

dx : float

Uniform grid step size (only for FD method).

N : int

A total number of the grid points.

Np : int

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

Ne : int

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

w_d : numpy 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:

f : numpy array, dtype=float

Function values evaluated on the grid

result : numpy 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:

BoundCond : dictionary

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:

peq : numpy array, dtype=float

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

D : float

Noise magnitude.

q : numpy array, dtype=float

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

w : numpy array, dtype=float

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

mode : str

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

fr : numpy 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.

Nv : int

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

Returns:

lQ : numpy array (Nv,), dtype=float

The least Nv eigenvalues for the eigenvalue problem of H0 operator.

QxOrig : numpy array (Nv,Nv), dtype=float

The corresponding scaled eigenvectors

Qx : numpy 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:

peqs : numpy array

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

grid_params : dictionary

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_index : int

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

stop_index : int

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:

data1 : dictionary

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.

data2 : dictionary

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_params : dictionary

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_thres : float

Threshold KL divergence for determining optimal feature complexity.

FC_radius : float

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

FC_final : float

Maximum feature complexity on shared feature complexity axis

FC_stride : float

Feature complexity resolution on shared feature complexity axis

Returns:

FC_shared : numpy array, dtype=float

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

KL : numpy array, dtype=float

1D array of KL divergences between the two models.

FC_opt_ind : int

Index of optimal KL divergence in FC_shared array

ind1_shared : numpy array, dtype=int

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

ind2_shared : numpy 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, arxiv 2020 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:

results : dictionary

Dictionary with the results (fitted models)

em : EnergyModel

An instance of the EnergyModel class.

iterations : numpy array

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

Returns:

FCs_array : numpy array

A 1D array of feature complexities.

lQ_array : numpy array

A 2D array of H0 eigenvalues.

Qx_array : numpy 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:

em : EnergyModel

An instance of the EnergyModel class.

peq : numpy array

Equilibirum probability distribution.

p0 : numpy array

Distribution of the initial latent states.

D : float

Noise magnitude

Returns:

FC : float

Feature complexity.

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

Calculate JS divergence between two distributions

Parameters:

p1 : numpy array

The first distribution, 1D array

p2 : numpy array

The second distribution, 1D array

weights: numpy array

1D array of SEM weights

mode : ENUM(‘normalized’,’unnormalized’)

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

Returns:

JS : float

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:

peq1 : numpy array

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

D1 : float

Noise magnitude for the first dynamics.

p01 : numpy array

1D array of p0 distribution for the first dynamics.

peq2 : numpy array

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

D2 : float

Noise magnitude for the second dynamics.

p02 : numpy array

1D array of p0 distribution for the second dynamics.

weights : numpy array

1D array of SEM weights

lQ1 : numpy array

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

Qx1 : numpy array

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

lQ2 : numpy array

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

Qx2 : numpy array

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

terminal_time : float, optional

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

number_of_samples : int, optional

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

Returns:

out : float

JS divergence between the two time-dependent distributions.