Mismatch with results of other solvers

[Emma]

Dear Julien, 

 

I did a parameter optimization of a model I am using in ANNarchy. I optimized six parameters (tau, w, tau_a, b, mp0, a0) using scipy.least_squares(). 

 

When I used the same parameters in an ANNarchy simulation the results are completely off. I think I must be getting something fundamental very wrong. If you could give advice, I would apprechiate that so much. 

 

This is a minimal version oft he problem

 

# These are the parameters I want to use in the simulation:

 

from ANNarchy import *
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

 

# initial values
mp0 = 0.009433
a0 = 0.002228
y0 = [mp0, a0]

# Parameters
w = 1.205478
tau = 12.681473
tau_a = 1.071650
b = 12.189142

 

# Time-varying input (inp)
inp = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
t = np.arange(0,31,1)

 

 

# This is the solution with scipy ivp: 

 

def model(t, y, w, tau, tau_a, b, inp):
    mp, a = y

    # Find the index corresponding to the current time in the inp array
    max_idx = len(inp) - 1
    idx = min(int(np.floor(t)), max_idx)  # Ensure index is within bounds of inp

    # Differential equations
    dmp_dt = (-mp + w * inp[idx] - a * b) / tau
    da_dt = (-a + mp) / tau_a

    return dmp_dt, da_dt

t_span = (0,max(t))
sol = solve_ivp(model, t_span, y0, t_eval=t, args=(w, tau, tau_a, b, inp))
mp_sol_ivp = sol.y[0]

 

This is the solution in ANNarchy:

 

setup(dt=1.0)
setup(method='explicit') 

tau = Constant('tau',tau)
b = Constant('b',b)
tau_a = Constant('tau_a',tau_a)
w = Constant('w', w)

neuron_template = Neuron(
    equations=
    """
         dmp/dt = (-mp + sum(exc) -a * tau) / 2 : init=0.009433
         da/dt = (-a + mp) / tau_a : init=0.002228
         r = mp
     """)

# input to neuron
inp_neuron = TimedArray(rates=inp.reshape(-1, 1))

neuron = Population(geometry=1, neuron=neuron_template, name='neuron')
input_syn = Projection(pre=inp_neuron, post=neuron, target='exc').connect_one_to_one(w)

mon = Monitor(neuron,'r')

net = Network(everything=True)
net.compile()
net.simulate(max(t)+1)

neuron_r = net.get(mon).get('r')

 

 

This is the result:

 plt.figure()

plt.plot(mp_sol_ivp, label='SCiPy ivp')
plt.plot(neuron_r, label='ANNarchy')
plt.ylim(-0.05,0.2)
plt.legend()
plt.show()

 

 

I double-checked the scipy.ivp solution with other solvers. They all yield comparable results. Do you have any idea what I am doing wrong using ANNarchy? 

 

Thank you so much!

Emma 

[julien...@gmail.com]

Hi Emma,

the equation for mp in ANNarchy should be: dmp/dt = (-mp + sum(exc) -a * b) / tau, otherwise it is not comparable with the scipy code.

Furthermore, scipy uses rk45 as a solver by default (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html), a version of runge-kutta order 4, while the default "explicit" Euler method in ANNarchy is just order 1. If you set ANNarchy to use Runge-Kutta 4 ("rk4", quite close to rk45, but not entirely the same):

setup(dt=1.0, method='rk4')

you get  quite close results:

It is important to use higher-order methods when the time constants (tau_a) are close to the time step dt (1ms in both cases), otherwise you get strong numerical issues. You could also decrease dt, but that would lead to more computations.

Best regards

Julien 

[Emma]

Hi Julien, 

the equation for mp in ANNarchy should be: dmp/dt = (-mp + sum(exc) -a * b) / tau, otherwise it is not comparable with the scipy code.

Im so sorry: That got messed up while I was copy pasting around to put together the example. That was not actually the "real" issue. But thank you for you patience. 

Furthermore, scipy uses rk45 as a solver by default (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html), a version of runge-kutta order 4, while the default "explicit" Euler method in ANNarchy is just order 1. If you set ANNarchy to use Runge-Kutta 4 ("rk4", quite close to rk45, but not entirely the same):

Thank you so much. That is the perfect solution. I don't want to decrease dt, especially because I might need to monitor several variables at every step. 

Best, 

Emma