Continuous data generation

[pyrates-neuroscience]

Continuous data generation

[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!