[SciPy-User] Optimizing odeint without a for loop

[Moore, Eric (NIH/NIDDK) [F]]

From: Barrett B [mailto:barrett.n.b at gmail.com] 
Sent: Wednesday, June 04, 2014 12:18 AM
To: scipy-user at scipy.org
Subject: [SciPy-User] Optimizing odeint without a for loop

This is an initial attempt to model a neural network with two differential variables per cell--voltage, and the response variable n. I can go back and add more later as needed. I have a question about this code (obviously, all constants are given):

--------

#ODE 
def f(X, t):
   N = len(X)/2
   dV = np.zeros(N); dn = np.zeros(N)
   for i in range(N):
   E = X[i]; n = X[N+i] #dummy
      n_inf_E = 1/(1 + np.exp((E_half_n - E)/k_n))
      m_inf_E = 1/(1 + np.exp((E_half_m - E)/k_m))
      dV[i] = (I - g_L*(E - E_L) - g_Na*m_inf_E*(E - E_Na) - g_K*n*(E - E_K))/C
      for j in range(N):
          dV[i] += eps*connex[j, i]*X[j]
      dn[i] = (n_inf_E - n)/tau
    return np.concatenate((dV, dn))

connex = np.matrix([[-1,1,0],
[1,-2,1],
[0,1,-1]]) #connection matrix

t = np.arange(0, stopTime, timeStep)
V0 = np.array([0, -20, -50])
n0 = np.array([0.2, 0.4, 0.7]); N = len(n0)

soln = odeint(f, np.concatenate((V0, n0)), t)

-----------------
Is there a way to do this without the "for i in range (N)" loop; i.e., can I run through all the dV[i]'s in a more efficient method?
 
 
Untested:

def f(X, t):
    N = len(X)/2
    E = X[:N]
    n = X[1:N+1]
    n_inf_E = 1/(1 + np.exp((E_half_n - E)/k_n))
    m_inf_E = 1/(1 + np.exp((E_half_m - E)/k_m))
    dV = (I - g_L*(E - E_L) - g_Na*m_inf_E*(E - E_Na) - g_K*n*(E - E_K))/C
    dV += eps * np.dot(X[:N], connex)
    dn = (n_inf_E - n)/tau
    return np.concatenate((dV, dn))

The basic idea is to operate on full arrays at once rather than looping over them.  See, for instance, https://scipy-lectures.github.io/intro/numpy/operations.html

Eric