brainpy / BrainPy

Brain Dynamics Programming in Python
https://brainpy.readthedocs.io/
GNU General Public License v3.0
537 stars 94 forks source link

"show_code" does not work when ode integration method is "exp_euler" #457

Closed CloudyDory closed 7 months ago

CloudyDory commented 1 year ago

In the following code:

import brainpy as bp
import brainpy.math as bm

Iext=10.;   ENa=50.;   EK=-77.;   EL=-54.387
C=1.0;      gNa=120.;  gK=36.;    gL=0.03

def dm(m, t, V):
    alpha = 0.1 * (V + 40) / (1 - bm.exp(-(V + 40) / 10))
    beta = 4.0 * bm.exp(-(V + 65) / 18)
    dmdt = alpha * (1 - m) - beta * m
    return dmdt

def dh(h, t, V):
    alpha = 0.07 * bm.exp(-(V + 65) / 20.)
    beta = 1 / (1 + bm.exp(-(V + 35) / 10))
    dhdt = alpha * (1 - h) - beta * h
    return dhdt

def dn(n, t, V):
    alpha = 0.01 * (V + 55) / (1 - bm.exp(-(V + 55) / 10))
    beta = 0.125 * bm.exp(-(V + 65) / 80)
    dndt = alpha * (1 - n) - beta * n
    return dndt

def dV(V, t, m, h, n, Iext):
    I_Na = (gNa * m ** 3.0 * h) * (V - ENa)
    I_K = (gK * n ** 4.0) * (V - EK)
    I_leak = gL * (V - EL)
    dVdt = (- I_Na - I_K - I_leak + Iext) / C
    return dVdt

hh_derivative = bp.JointEq([dV, dm, dh, dn])

integral_1 = bp.odeint(hh_derivative, method='rk2', show_code=True)
integral_2 = bp.odeint(hh_derivative, method='exp_euler', show_code=True)

Running the integral_1 = ... line outputs the following in console:

def brainpy_itg_of_ode7_joint_eq(V, m, h, n, t, Iext, dt=0.1):
  dV_k1, dm_k1, dh_k1, dn_k1 = f(V, m, h, n, t, Iext)
  k2_V_arg = V + dt * dV_k1 * 0.6666666666666666
  k2_m_arg = m + dt * dm_k1 * 0.6666666666666666
  k2_h_arg = h + dt * dh_k1 * 0.6666666666666666
  k2_n_arg = n + dt * dn_k1 * 0.6666666666666666
  k2_t_arg = t + dt * 0.6666666666666666
  dV_k2, dm_k2, dh_k2, dn_k2 = f(k2_V_arg, k2_m_arg, k2_h_arg, k2_n_arg, k2_t_arg, Iext)
  V_new = V + dV_k1 * dt * 0.25 + dV_k2 * dt * 0.75
  m_new = m + dm_k1 * dt * 0.25 + dm_k2 * dt * 0.75
  h_new = h + dh_k1 * dt * 0.25 + dh_k2 * dt * 0.75
  n_new = n + dn_k1 * dt * 0.25 + dn_k2 * dt * 0.75
  return V_new, m_new, h_new, n_new

However, running the integral_2 = ... line only outputs this:

{'f': <brainpy._src.integrators.joint_eq.JointEq object at 0x0000026BDD214610>}

If the formatted code is not shown, how can I determine the position of variables and parameters when calling the joint equation?

chaoming0625 commented 1 year ago

Thanks for your report.

Actually, exp_euler does not support show_code = True/False. The positions of variables and parameters are the same as you are using rk2. This is determined by the brainpy.JointEq, which sequentially adds variables and then parameters.

Indeed, documentation is needed for this functionality. Thanks.