I am implementing an ODE solver, where the user provides rates and coefficients as a string. ODE solver has to work with vectors. The best implementation I got so far is the following:
# k_1 = dt*dcdt(C0, dt);
# k_2 = dt*dcdt(C0+0.5*k_1, dt);
# k_3 = dt*dcdt(C0+0.5*k_2, dt);
# k_4 = dt*dcdt(C0+k_3, dt);
# C_new = C0 + (k_1+2*k_2+2*k_3+k_4)/6;
import numpy as np
import numexpr as ne
dt = 0.1
coef = {'k': 1}
rates = {'R1': 'k*y1*y2'}
dcdt = {'y1': '-4 * R1', 'y2': '-R1'}
C0 = {'y2': np.random.rand(400), 'y1': np.random.rand(400)}
def k_loop(conc):
rates_num = {}
for k in rates:
rates_num[k] = ne.evaluate(rates[k], {**coef, **conc})
dcdt_num = {}
for k in dcdt:
dcdt_num[k] = ne.evaluate(dcdt[k], rates_num)
Kn = {}
for k in C0:
Kn[k] = dt * dcdt_num[k]
return Kn
def sum_k(A, B, b):
C = {}
for k in A:
C[k] = A[k] + b * B[k] * dt
return C
def result(C_0, k_1, k_2, k_3, k_4):
C_new = {}
for k in C_0:
C_new[k] = C_0[k] + (k_1[k] + 2 * k_2[k] + 2 * k_3[k] + k_4[k]) / 6
return C_new
k1 = k_loop(C0)
k2 = k_loop(sum_k(C0, k1, 0.5))
k3 = k_loop(sum_k(C0, k2, 0.5))
k4 = k_loop(sum_k(C0, k3, 1))
C_new = result(C0, k1, k2, k3, k4)
This ODE solver I am planning to use inside of the Object-oriented PDE solver. The PDE solver will provide C0
vectors. Therefore, the performance is crucial. I used numexpr.evaluate() for simplified user input of the rates. Can it be implemented in vectorized form? If you know the better implementation of the solver with such simplified input, please, let me know.
Metrics:
eval():
1000 loops, best of 3: 372 µs per loop
ne.evaluate():
1000 loops, best of 3: 377 µs per loop
@Graipher example:
1000 loops, best of 3: 452 µs per loop
np.array
inC0
do anything for you? Why not just lists? I don't see any other use ofnumpy
operations. \$\endgroup\$ – hpaulj Jan 16 '17 at 8:07eval
is considered to be dangerous and slow.ast.literal_eval
is safer, though my impression is that it may be slower.sympy
is another option if you math expressions. It might help if you gave examples ofrates
. Just how general do you need to go? \$\endgroup\$ – hpaulj Jan 16 '17 at 17:18