error while using delay DE

[NinaOmejc]

Hi,

I have been using PyRates for modeling and came across a question which I don't know how to resolve and I would really appreaciate some help.

I've been using PyRates to model coupled kuramoto oscillators, as shown in the appended file. I run the model with such command:

phases = model.run(outputs=outputs_dict, inputs={**noise_inputs, **burst_timing_inputs}, simulation_time=sett['simulation_time'], step_size=solver_step_size, sampling_step_size=sampling_step_size, solver=solver, method=solver_method, clear=True, in_place=False, backend='default')

where delays is a matrix of (number_nodes * number_nodes). When the delays is a matrix of only zeros, the simulations run without error, but when I instead use the nonzero matrix, then I get an error:

/// Compilation Progress

(1) Translating the circuit template into a networkx graph representation... ...finished. (2) Preprocessing edge transmission operations... C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\ir\circuit.py:555: PyRatesWarning: PyRates detected an edge definition that implies to represent the model as a delayed differential equation system. PyRates will thus attempt to access the history of the source variable s of operator kmo_oper_coupling on node kmo_edge. Note that this requires s to be a state-variable, i.e. a variable defined by a differential equation. warn(PyRatesWarning(f'PyRates detected an edge definition that implies to represent the model as a ' ...finished. (3) Parsing the model equations into a compute graph... Traceback (most recent call last):

....

File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\backend\computegraph.py", line 897, in _generate_unique_label return self._generate_unique_label(new_label) [Previous line repeated 2970 more times] File "C:\Users\NinaO\AppData\Local\Programs\Python\Python310\lib\site-packages\pyrates\backend\computegraph.py", line 891, in _generate_unique_label if label in self.nodes: RecursionError: maximum recursion depth exceeded ///

I guess part in the warning message that I bolded gives an answer but I do not know how to actually approach this issue, as s is not actually a state variable. Also, as far as I was looking at your tutorial on delay coupling DE, the circuit in your case stayed the same, and only the delay variable was added to the edge. What am I missing?

Thank you again for your time! Best, Nina

kuramoto_model.zip

[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