error while using delay DE

[pyrates-neuroscience]

error while using delay DE

[Richert]

Dear Nina,

sorry for the delayed response, I was on vacation in Europe!
I believe your problem is caused by the choice of numerical solver you use.
If you want to use a numerical solver method with an adaptive step-size, you will have to implement a proper delay differential equation system, as explained in
this tutorial.
This requires that you define delayed versions of state-variables rather than use the edge method that you used so far.

If you define delays on edges, PyRates will attempt the following:

(1) If you also define the keyword spread, PyRates will approximate each edge delay via a convolution of the edge source variable with a gamma kernel. In this case, the delay keyword is treated as the mode of that gamma kernel rather than as a discrete offset in time.

(2) If no spreads are defined, PyRates will treat delay as a discrete offset of the source variable in time. If you choose an adaptive step-size solver, PyRates will try to represent the entire model as a delayed differential equation system. This will only work if the source variable of the edge is a state variable.

(3) If you want to couple non-state variables via delays that is certainly possible. However, you will have to choose the standard Euler solver for solving the differential equation system, because PyRates will build a buffer matrix for each delay and store the values of the delayed variables in there. This requires a discretization of the buffer matrix via the integration step-size. Thus, adaptive integration step-sizes break that method.

I hope that one of these options works in your case!

Best, Richard

[NinaOmejc]

Dear Richard,

thank you for the explanation!
I think the edge delays are the only option for me, as the model only includes delays between the nodes and I don't see a way how that can be represented as a state variable.

Do I understand correctly that if the Gamma kernel is used (by applying spread), then the offset is not descrete anymore but rather estimated over a continuous range? In that case, the adaptive step-size should work? but it seems I'm getting something wrong, as it doesnt work in my case.

I have tried to use discrete solvers (Euler and Heun) but I am having troubles using using them for simulation, when I add the input. I was able to reproduce the error on a simpler example, of which I attach the code. The simulation works for examples 1 (using RK solver, no input), example 2 (using Heun, no input) and example 4 (using LSODA, with input), but not for example 3 (Heun, with input), where I get the error inside "(3) Running the simulation ...".

`
Traceback (most recent call last):
File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\IPython\core\interactiveshell.py", line 3508, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "", line 1, in
model_simulation_results_3 = model.run(simulation_time=simulation_time, step_size=1e-3, outputs=outputs, inputs=noise_input)
File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\frontend\template\circuit.py", line 473, in run
outputs = net._ir.run(simulation_time=simulation_time, solver=solver, sampling_step_size=sampling_step_size,
File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\ir\circuit.py", line 1046, in run
results = self.graph.run(func=func, func_args=func_args, T=simulation_time, dt=self._dt, dts=sampling_step_size,
File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\backend\computegraph.py", line 472, in run
results, times = self.backend.run(func=func, func_args=func_args, T=T, dt=dt, dts=dts, solver=solver, **kwargs)
File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\backend\base\base_backend.py", line 343, in run
results = self._solve(solver=solver, func=func, args=func_args[2:], T=T, dt=dt, dts=dts, y0=y0, t0=t0,
File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\backend\base\base_backend.py", line 412, in _solve
results = self._solve_euler(func, args, T, dt, dts, y0, t0)
File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\backend\base\base_backend.py", line 448, in _solve_euler
rhs = func(step, y, *args)
File "D:\Experiments\modeling_EEG\pyrates_run.py", line 10, in vector_field
s_ext_timed_input[0] = s_ext_input[t]

IndexError: index 512 is out of bounds for axis 0 with size 512
`

Any thoughts about the error? I checked your examples online and it doesn't seem I am missing some settings? I have also tried to decrease the size of the input, but that didn't help.

Thank you very much.
Best,
Nina

discrete_solver_with_input_test.zip

[Richert]

Dear Nina,

it is correct that if you use the spread argument, each edge is modeled as a convolution of the edge source variable with a gamma kernel. So there is no discrete delay anymore. Since this is not part of the example you shared, I am not sure what the concrete issue is that you are having with it. But since the convolution is translated into a set of ODEs, adaptive step-size solvers should work without any issues.

Regarding the example you provided, I see that there is mismatch between the sampling_rate = 256 you defined and the integration step-size of step_size=1e-3. Extrinsic inputs have to be defined for each simulation time step, of which you have steps = int(simulation_time / step_size) = 2000. Your extrinsic input vectors have only 512 entries, however. So at step 513 the error you received is thrown.

Hope that helps!
Best, Richard

[NinaOmejc]

Dear Richard,

oh, that was a silly mistake. I actually knew (saw in your examples) about the size of the input requirement at some point but I guess I forgot or didnt focus on that anymore, because the simulation worked for the LSODA solver ( with the same short input ).

Actually, the delay also sometimes works now with the spread argument and adaptible step size solvers. I have a general issue of jobs failing (around 1/3 jobs fails) but I think I have to debug that on my own because I have some randomness involved in generation of the models..

Thank you again for all the help!

Best,
Nina