From 9cf948005e5d73d065d1dd7fdb3e662648a39a15 Mon Sep 17 00:00:00 2001 From: caglorithm Date: Wed, 30 Aug 2023 11:40:56 +0200 Subject: [PATCH] aln: fix index of input from i to i-1 --- neurolib/models/aln/timeIntegration.py | 46 ++++++++++++-------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/neurolib/models/aln/timeIntegration.py b/neurolib/models/aln/timeIntegration.py index c5b7346a..fb09dad2 100644 --- a/neurolib/models/aln/timeIntegration.py +++ b/neurolib/models/aln/timeIntegration.py @@ -7,7 +7,7 @@ def timeIntegration(params): """Sets up the parameters for time integration - + Return: rates_exc: N*L array : containing the exc. neuron rates in kHz time series of the N nodes rates_inh: N*L array : containing the inh. neuron rates in kHz time series of the N nodes @@ -397,12 +397,11 @@ def timeIntegration_njit_elementwise( noise_exc, noise_inh, ): - # squared Jee_max - sq_Jee_max = Jee_max ** 2 - sq_Jei_max = Jei_max ** 2 - sq_Jie_max = Jie_max ** 2 - sq_Jii_max = Jii_max ** 2 + sq_Jee_max = Jee_max**2 + sq_Jei_max = Jei_max**2 + sq_Jie_max = Jie_max**2 + sq_Jii_max = Jii_max**2 # initialize so we don't get an error when returning rd_exc_rhs = 0.0 @@ -416,7 +415,6 @@ def timeIntegration_njit_elementwise( ### integrate ODE system: for i in range(startind, startind + len(t)): - if not distr_delay: # Get the input from one node into another from the rates at time t - connection_delay - 1 # remark: assume Kie == Kee and Kei == Kii @@ -430,13 +428,12 @@ def timeIntegration_njit_elementwise( # loop through all the nodes for no in range(N): - # To save memory, noise is saved in the rates array noise_exc[no] = rates_exc[no, i] noise_inh[no] = rates_inh[no, i] - mue = Jee_max * seem[no] + Jei_max * seim[no] + mue_ou[no] + ext_exc_current[no, i] - mui = Jie_max * siem[no] + Jii_max * siim[no] + mui_ou[no] + ext_inh_current[no, i] + mue = Jee_max * seem[no] + Jei_max * seim[no] + mue_ou[no] + ext_exc_current[no, i - 1] + mui = Jie_max * siem[no] + Jii_max * siim[no] + mui_ou[no] + ext_inh_current[no, i - 1] # compute row sum of Cmat*rd_exc and Cmat**2*rd_exc rowsum = 0 @@ -447,33 +444,35 @@ def timeIntegration_njit_elementwise( # z1: weighted sum of delayed rates, weights=c*K z1ee = ( - cee * Ke * rd_exc[no, no] + c_gl * Ke_gl * rowsum + c_gl * Ke_gl * ext_exc_rate[no, i] + cee * Ke * rd_exc[no, no] + c_gl * Ke_gl * rowsum + c_gl * Ke_gl * ext_exc_rate[no, i - 1] ) # rate from other regions + exc_ext_rate z1ei = cei * Ki * rd_inh[no] z1ie = ( - cie * Ke * rd_exc[no, no] + c_gl * Ke_gl * ext_inh_rate[no, i] + cie * Ke * rd_exc[no, no] + c_gl * Ke_gl * ext_inh_rate[no, i - 1] ) # first test of external rate input to inh. population z1ii = cii * Ki * rd_inh[no] # z2: weighted sum of delayed rates, weights=c^2*K (see thesis last ch.) z2ee = ( - cee ** 2 * Ke * rd_exc[no, no] + c_gl ** 2 * Ke_gl * rowsumsq + c_gl ** 2 * Ke_gl * ext_exc_rate[no, i] + cee**2 * Ke * rd_exc[no, no] + + c_gl**2 * Ke_gl * rowsumsq + + c_gl**2 * Ke_gl * ext_exc_rate[no, i - 1] ) - z2ei = cei ** 2 * Ki * rd_inh[no] + z2ei = cei**2 * Ki * rd_inh[no] z2ie = ( - cie ** 2 * Ke * rd_exc[no, no] + c_gl ** 2 * Ke_gl * ext_inh_rate[no, i] + cie**2 * Ke * rd_exc[no, no] + c_gl**2 * Ke_gl * ext_inh_rate[no, i - 1] ) # external rate input to inh. population - z2ii = cii ** 2 * Ki * rd_inh[no] + z2ii = cii**2 * Ki * rd_inh[no] sigmae = np.sqrt( 2 * sq_Jee_max * seev[no] * tau_se * taum / ((1 + z1ee) * taum + tau_se) + 2 * sq_Jei_max * seiv[no] * tau_si * taum / ((1 + z1ei) * taum + tau_si) - + sigmae_ext ** 2 + + sigmae_ext**2 ) # mV/sqrt(ms) sigmai = np.sqrt( 2 * sq_Jie_max * siev[no] * tau_se * taum / ((1 + z1ie) * taum + tau_se) + 2 * sq_Jii_max * siiv[no] * tau_si * taum / ((1 + z1ii) * taum + tau_si) - + sigmai_ext ** 2 + + sigmai_ext**2 ) # mV/sqrt(ms) if not filter_sigma: @@ -531,10 +530,10 @@ def timeIntegration_njit_elementwise( seim_rhs = ((1 - seim[no]) * z1ei - seim[no]) / tau_si siem_rhs = ((1 - siem[no]) * z1ie - siem[no]) / tau_se siim_rhs = ((1 - siim[no]) * z1ii - siim[no]) / tau_si - seev_rhs = ((1 - seem[no]) ** 2 * z2ee + (z2ee - 2 * tau_se * (z1ee + 1)) * seev[no]) / tau_se ** 2 - seiv_rhs = ((1 - seim[no]) ** 2 * z2ei + (z2ei - 2 * tau_si * (z1ei + 1)) * seiv[no]) / tau_si ** 2 - siev_rhs = ((1 - siem[no]) ** 2 * z2ie + (z2ie - 2 * tau_se * (z1ie + 1)) * siev[no]) / tau_se ** 2 - siiv_rhs = ((1 - siim[no]) ** 2 * z2ii + (z2ii - 2 * tau_si * (z1ii + 1)) * siiv[no]) / tau_si ** 2 + seev_rhs = ((1 - seem[no]) ** 2 * z2ee + (z2ee - 2 * tau_se * (z1ee + 1)) * seev[no]) / tau_se**2 + seiv_rhs = ((1 - seim[no]) ** 2 * z2ei + (z2ei - 2 * tau_si * (z1ei + 1)) * seiv[no]) / tau_si**2 + siev_rhs = ((1 - siem[no]) ** 2 * z2ie + (z2ie - 2 * tau_se * (z1ie + 1)) * siev[no]) / tau_se**2 + siiv_rhs = ((1 - siim[no]) ** 2 * z2ii + (z2ii - 2 * tau_si * (z1ii + 1)) * siiv[no]) / tau_si**2 # -------------- integration -------------- @@ -596,7 +595,6 @@ def interpolate_values(table, xid1, yid1, dxid, dyid): @numba.njit(locals={"idxX": numba.int64, "idxY": numba.int64}) def lookup_no_interp(x, dx, xi, y, dy, yi): - """ Return the indices for the closest values for a look-up table Choose the closest point in the grid @@ -636,9 +634,9 @@ def lookup_no_interp(x, dx, xi, y, dy, yi): return idxX, idxY + @numba.njit(locals={"xid1": numba.int64, "yid1": numba.int64, "dxid": numba.float64, "dyid": numba.float64}) def fast_interp2_opt(x, dx, xi, y, dy, yi): - """ Returns the values needed for interpolation: - bilinear (2D) interpolation within ranges,