Continuous data generation

[ffriese]

Dear Pyrates Developers, I have a use-case in which I would like to initialize a Circuit Template once and then call the run-function multiple times (maybe even alter variables in between) to generate a continuous signal.

This was possible at least up to version 0.9.5 with a minimal fix to the numpy-backend, which used to incorrectly update the self.vars['y'] variable after each pass (see #e800a9ca), and appeared to be working flawlessly (even with variable alterations in between, i.e. I generated a continuous signal where the length of the inhibitory PSP changed over time, and the signal showed the expected characteristics changing over time).

Today I tried to update to the new version of PyRates and after I fixed the minor incompabilities in my code, I realized that the underlying implementation changed considerably (i.e. changes to the apply, compile and run mechanism, the new generate_run_func-mechanism, ...). Thus, I would like to ask whether there is a straightforward way to approach this problem, or if I should maybe just go back to version 0.9.5 to keep from having to deeply dig into the new mechanics and possibly cause different problems with my hacks.

Thanks in advance

[Richert]

Hi!

Indeed, some of the recent changes to PyRates made this implementation a bit more tricky. The main reason behind this is that PyRates is developing towards a software for model definition and code generation, mainly. To ensure that any update to a PyRates variable is also tracked properly by the code generator, no changes can be made to a model after apply has been called on the model instance. There are a few ways to do what you want to do though, I believe. As long as you do not work with edge delays, you can implement parameter jumps within a simulation as follows:

import matplotlib.pyplot as plt
from pyrates import OperatorTemplate, NodeTemplate, CircuitTemplate

# define simple rate operator
eqs = [
    "r' = (r0-r)/tau + tanh(r_in)"
]
variables = {
    "r": "output(0.0)",
    "r0": 1.0,
    "tau": 10.0,
    "r_in": "input(0.0)"
}

# set up node template
op = OperatorTemplate(name="rate_op", equations=eqs, variables=variables, path=None)
node = NodeTemplate(name="rate_node", operators=[op], path=None)

# define an edge
edges = [
    ("p/rate_op/r", "p/rate_op/r_in", None, {"weight": -2.0})
]

# define the circuit
net = CircuitTemplate(name="rate_pop", nodes={"p": node}, edges=edges, path=None)

# perform initial simulation
T = 10.0
dt = 1e-3
dts = 1e-2

# perform simulation
res = net.run(simulation_time=T, step_size=dt, sampling_step_size=dts, solver="scipy", method="RK45",
              outputs={"r": "p/rate_op/r"}, in_place=False)

# update the time constant tau
net.update_var(node_vars={"p/rate_op/tau": 5.0})

# update the initial value of r
net.update_var(node_vars={"p/rate_op/r": res.iloc[-1, 0]})

# perform second simulation
res2 = net.run(simulation_time=T, step_size=dt, sampling_step_size=dts, solver="scipy", method="RK45",
               outputs={"r": "p/rate_op/r"})

# combine results
res2.index += res.index[-1]
res = res.append(res2)

# visualize combined results
plt.plot(res)
plt.show()

In this minimal example, the important part is to use the keyword argument in_place=False during the simulation, which will ensure that the apply method is not used on the CircuitTemplate instance itself, but on a deepcopy of it instead. That way, you can call CircuitTemplate.run multiple times on the instance and use CircuitTemplate.update_var in between to change constants. Note that you will have to set the new initial values of the state variables manually though. There should be an easy way to implement a more convenient way to do the latter (i.e. adding a keyword argument y0 to the run method that could be used as initial state vector if provided). Feel free to come up with something like that and send a pull request! An alternative would be of course to just use an extrinsic input instead. The following code provides exactly the same result:

import matplotlib.pyplot as plt
from pyrates import OperatorTemplate, NodeTemplate, CircuitTemplate
import numpy as np

# define simple rate operator
eqs = [
    "r' = (r0-r)/(tau+tau_ext) + tanh(r_in)"
]
variables = {
    "r": "output(0.0)",
    "r0": 1.0,
    "tau": 10.0,
    "r_in": "input(0.0)",
    "tau_ext": "input(0.0)"
}

# set up node template
op = OperatorTemplate(name="rate_op", equations=eqs, variables=variables, path=None)
node = NodeTemplate(name="rate_node", operators=[op], path=None)

# define an edge
edges = [
    ("p/rate_op/r", "p/rate_op/r_in", None, {"weight": -2.0})
]

# define the circuit
net = CircuitTemplate(name="rate_pop", nodes={"p": node}, edges=edges, path=None)

# perform initial simulation
T = 20.0
dt = 1e-3
dts = 1e-2

inp = np.zeros((int(T/dt),))
inp[int(10/dt):] = -5.0

# perform simulation
res = net.run(simulation_time=T, step_size=dt, sampling_step_size=dts, solver="scipy", method="RK45",
              outputs={"r": "p/rate_op/r"}, inputs={"p/rate_op/tau_ext": inp})

# visualize combined results
plt.plot(res)
plt.show()

[ffriese]

Hi Richard, once again I am amazed at your quick and thorough response. I think your examples provide me with everything I need to reach my goal, I will let you know how it went tomorrow on monday! Thank you!

[ffriese]

Hi Richard, I tried to get my project running with your examples, and found your second example extremely helpful for my immediate needs. It allows me to easily change the parameter values continously as well (I used to change them in small, but somewhat more "discrete" steps every second only).

I might get to a point where I would want to keep the generation running indefinitely until I stop it explicitely - in which case I will try to adapt something from your first example. I will gladly send a pull request when/if I get there.

Thanks again!