pyrates-neuroscience/pyrates

Open-source, graph-based Python code generator and analysis toolbox for dynamical systems (pre-implemented and custom models). Most pre-implemented models belong to the family of neural population models.

pyrates's Issues

Allow variable initialization

Right now I couldn't find a way how to initialize variables. By default, all variables (input, output, variable) are initialized at 0. When using a variable in equations as e.g.

tau_h = ...
h = (1/tau_h) ...

the compile method (I must add correctly) throws

ValueError: Result of operation (numpy_multiply) contains NaNs or infinite values.

It would be nice to give the users an option to initialize integration variables with some non-zero values.

Operator templates with identical equations cause KeyError in edge vectorization

We used a circuit that featured two edges with one operator each. These operators had each their own template with identical equations but different variable values. The working example is given in
model_templates/wilson_cowan/simple_wilsoncowan.yml (commit f246f0e )

As expected, both edge operators are mapped on the same operator, due to identical equations. However, something must be going wrong there, because the vectorization process later on fails to correctly map both edge operators on the same vectorized coupling_node. This is evident, because the process raises a KeyError, when processing the second operator, but looking for the the first operator (by name).

The problem can be fixed by making one operator a sub-template of the other, inheriting the equations and changing the variable values. The resulting template can be found in the linked example.

Expected behaviour is that two operators with identical equations (and possibly variable type definitions) are treated as two instances of the same template, even if the are not related by inheritance.

Continuous data generation

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

Add possibility to CircuitTemplate.apply() method to change constants

This is a suggestion on how to improve the flexibility of working with YAML templates for model definitions in PyRates.

I suggest to add an additional keyword argument to the pyrates.frontend.CircuitTemplate.apply() method. This keyword argument should be a dictionary, which allows to point to any variable on any operator on any node or edge in the graph via the standard slash-based notation and change its scalar value.

This would enable quick parameter changes within Python without the necessity of defining new yaml templates. Thus, interacting with the pre-implemented models in model_definitions becomes more intuitive, since users would not have to mess with the yaml templates in order to apply a simple parameter change. Finally, it would allow to change initial values of state variables more easily.

Merge ComputeGraph into CircuitIR

To improve separation between frontend and backend, I suggest to merge all ComputeGraph functionalities into the CircuitIR.

This includes:

• vectorization of graph nodes
• parsing of the graph nodes and edges into the backend
• a run method for forward simulations in time

To achieve this, I suggest the following:

• define a compile method for the CircuitIR and add all arguments of the ComputeGraph.__init__() method to this compile function
• add a vectorize keyword argument to the CircuitIR.__init__() which can be True or False for turning node vectorization of the graph on or off
• If vectorize is True, perform vectorization on the fly while adding the nodes to the network, by comparing new nodes to the already excisting ones and vectorizing those with identical operator graph structures and names
• To keep track of which of the original nodes resides in which vector, build a map dictionary that links the original node/op/var triplets to the respective vector indices of the new, vectorized nodes
• Define appropriate get_var and set_var methods that allow access to variables in the network. These should work based on the above described map. Also, since we allow a hierarchical node structure, these methods should allow to reference either all nodes or a specific one at each level of that hierarchy.

Remove unnecessary indexing operations and variable definitions from equation file

During the construction of the update equations file of any PyRates model (happens in the backend), it may happen that variables are extracted from the parameter or state variable vectors that are not required and thus cause unnecessary indexing operations. This should be fixed.

error while using delay DE

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

is there an option to get an ouput of the non-state variable?

Hi,

thank you for this tool, i find it very useful! I have one problem though, if you might help me out.

Below I attached the file with Kuramoto model as is more or less defined in your example script. When I run the model it works and it outputs the angles, if my outputs are:

outputs_dict = {f"{inode}_theta": f"{inode}/kmo_oper_basic/theta" for inode in model.nodes.keys()}
However, when I want to include the wrapper version as an ouput:

outputs_dict = {**{f"{inode}_theta": f"{inode}/kmo_oper_basic/theta" for inode in model.nodes.keys()},
         **{f"{inode}_wrap": f"{inode}/kmo_oper_sin_wrapper/wrap" for inode in model.nodes.keys()}}

i get the error, which is not really informative:

_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 29, in
phases = model.run(outputs=outputs_dict, simulation_time=simulation_time, step_size=step_size, sampling_step_size=sampling_step_size,
File "D:\Experiments\simulation_EEG\venv\lib\site-packages\pyrates\frontend\template\circuit.py", line 468, in run
outputs = net._ir.run(simulation_time=simulation_time, solver=solver, sampling_step_size=sampling_step_size,
File "D:\Experiments\simulation_EEG\venv\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 "D:\Experiments\simulation_EEG\venv\lib\site-packages\pyrates\backend\computegraph.py", line 464, in run
outputs[key] = self.state_var_indices[var]
KeyError: 'wrap'

I think the error is because my variable 'wrap' is not a state variable. Is there a way to nevertheless get the output of 'wrap', so the angle wrapped inside sin function? Or is there another problem?

Thank you for your help!

kuramoto_model.zip

Found test problem caused by changes of working directory in backend

While trying to understand differences between tests run on Windows and Linux I found that some problems (at least) are cause by the backend changing directories. @Richert : Can we change that?

issues when running the example codes

Hi, I was trying to run the codes from the documentation example:

from pyrates.frontend import CircuitTemplate
jrc = CircuitTemplate.from_yaml("model_templates.neural_mass_models.jansenrit.JRC")
results = jrc.run(simulation_time=1.0,
step_size=1e-2,
sampling_step_size=1e-2,
outputs={'V_pce': 'pc/rpo_e_in/Z',
'V_pci': 'pc/rpo_i/Z'},
backend='default',
solver='scipy')

It yields value error.

I run in Jupiter notebook with python 3.8 kernel. The spicy version is 1.7.3. The detail errors are shown below:

ValueError Traceback (most recent call last)

  1 from pyrates.frontend import CircuitTemplate
  2 jrc = CircuitTemplate.from_yaml("model_templates.neural_mass_models.jansenrit.JRC")

----> 3 results = jrc.run(simulation_time=1.0,
4 step_size=1e-2,
5 sampling_step_size=1e-2,

in pyrates/frontend/template/circuit.py in run(self, simulation_time, step_size, inputs, outputs, sampling_step_size, cutoff, solver, backend, vectorize, verbose, clear, in_place, **kwargs)
322
323 # apply template (translate into compute graph, optional vectorization process)
--> 324 net.apply(adaptive_steps=adaptive_steps, vectorize=vectorize, verbose=verbose, backend=backend,
325 step_size=step_size, **kwargs)
326

in pyrates/frontend/template/circuit.py in apply(self, adaptive_steps, label, node_values, edge_values, vectorize, verbose, **kwargs)
517
518 # group edges that should be vectorized
--> 519 old_edges = self.collect_edges(delay_info=True)
520 edge_col = self._group_edges(edges=old_edges)
521

in pyrates/frontend/template/circuit.py in collect_edges(self, delay_info)
799 edges.append((f"{c_scope}/{svar}", f"{c_scope}/{tvar}", template, edge_dict))
800 if delay_info:
--> 801 for i, (svar, tvar, template, edge) in enumerate(edges):
802 delayed = True if 'delay' in edge and edge['delay'] else False
803 edges[i] = (svar, tvar, template, edge, delayed)

ValueError: too many values to unpack (expected 4)

Perhaps no edge is generated? What are your thoughts?

Inheriting from an operator more than once causes operator key errors

I realized that there exists a bug in how operator keys are assigned, when multiple operators inherit from a common base operator:

base_op:
  base: OperatorTemplate
  equations:
    - "d/dt * R = -R/tau + J*I"
  variables:
    R:
      default: output
    I:
      default: input
    tau:
      default: 1.0
    J:
      default: 10.0

child_1:
  base: base_op
  variables:
    tau:
      default: 3.0

child_2
  base: base_op
  variables:
    tau:
      default: 5.0
    J:
      default: 20.0

In that case, only the first operator to inherit from the base operator receives the correct key ("child_1") when used in a CircuitTemplate, whereas the second operator will be called "child_1" as well.
I think that should be fixed and either all operators receive the key of the base operator, or their own.

Add variable setter method to CircuitIR

Similar to the CircuitIR.get_node_var(key) method, we should add a CircuitIR.set_node_var(key, val) method that allows to change values of constants and variables in the circuit even after the circuit has been vectorized or simulations have been performed.

After a change was made to a constant, the backend should check whether a recompilation of the system function is necessary. This might be the case if the respective constant is not an argument of the right-hand side function created during compilation, but its value was directly written to file.

The feature would make it easier to change parametrizations of circuits on the fly without having to re-initiate and -compile the whole circuit from scratch again.

issue with time-varying edge variables when multiple nodes are used

Dear Richard,

I've been trying to use your software to do some simulations of coupled Kuramoto oscillators. However, I'm having some issues.

I attach a file, with which I am running three coupled Kuramoto oscillators. As I want to have time varying coupling, I am trying to use
the 'trick' you mention in one of your github issues, that a variable can be given as an input. It seems this works for two oscillators, but for three, something seems weird. In the code, first oscillator is a slow one, trying to be coupled with the other faster two. I coupled 0-1 and 0-2 oscillators at the same time, with the same weights but when I plot the results, the coupling is seen in only 0-1 oscillators, but not between the 0-2 oscillators. I have tried to vary several properties, such as different timing of the second coupling, different frequencies, different edge combinations, either putting the coupling weight to the "weight" or to my time varying coupling variable "Ks" ... nothing worked the way as I would like it to... to have coupling between 0-1 oscillators, and then 0-2 oscillators (either during the same or during different time window).

Can you see the problem? I would really appreaciate your help.

Thank you and kind regards,
Nina
test_kuramoto.zip

Cache strings (equations, labels, etc.)

Python built-in support to cache strings with the intern function. By default in the CPython implementation strings are cached that qualify as names in Python. Equations and labels with slashes, however, do not fall into this rule. It makes sense though, to cache these strings, as they are used repetitively.

It might be possible to count references on these strings that can be used to speed up vectorization.

Implement logging

Summary

For better control about the amount of feedback the software gives and where it is written to, it would be helpful to use logging. This would replace any current print statements. We can flag outputs as INFO, DEBUG, or WARNING. Users can then control which level of information they would like to receive. Furthermore, it is possible to direct the feedback stream to various targets, including the standard output/terminal or a log file.

Details

• Possibly introduce some general setup/rc interface to set up logging globally for the user?

Instead of specifying file path in YAML templates, scan entire directory and then reference by names only?

I suggest to change how YAML templates are found and loaded.

Current implementation

• A YAML template is referenced by a path, either an actual path (/ or \\) or a Python-like import path (using .).
• In Python you can specify this path in a load function. Individual YAML templates can reference other templates (either in base attribute or when combining low-level templates (e.g. Operators) into higher-level templates (e.g. Nodes)).
• Template names are implicitly assumed to be unique throughout one session.

Problems with current implementation

• Paths can become rather long. If you want to use a path multiple times within the same template, it needs to be re-referenced or the whole thing becomes unreadable.
• Automatic refactoring does not work with YAML files (at least not in PyCharm). Moving templates or template files around can thus break your code. With paths being rather long, they also become annoying to maintain.
• Using the dot-notation invokes the Python import functionality, which only works if the template is placed in a python package. Using slash-notation relative or absolute paths can be specified. Absolute paths should be avoided in shared projects. Relative paths only work, if you always start from the same working directory.
• The slash notation collides with the internal slash notation used to walk through the model hierarchy. Template references must thus not be mixed with internal references (e.g. variables of operators of nodes, when defining edges).
• When the "base" template is the corresponding python class, it is referenced by name only which would imply that the template exists in the same file. This is inconsistent. (One could use the YAML syntax for defining external/Python objects, but that makes the file less readable again).

Suggested change

• Since templates names are assumed to be unique anyway, we can just use the plain names as unique identifiers: base.subdir.filename.TemplateName → TemplateName
• Instead of specifying a path for every single template, we can use an index to look these paths up. The index can be compiled automatically by recursively scanning a specified directory (or multiple d
  directories).
• Base classes (e.g. OperatorTemplate) are referenced by name as-is.

Advantages

• This way templates are easier to read.
• Template can be stored in an arbitrary folder that does not need to be related to the current project.

Disadvantages

• Templates are not self-contained anymore. You always need to specify a base directory before loading templates.
• The syntax rather reminds of Matlab than of Python.

Make IR objects (immutable) hashable

Making IR object immutable could produce a considerable performance boost, but also change a couple of interfaces.

Suggested changes

• CircuitIR, NodeIR, EdgeIR, OperatorGraphandOperatorIR` become (pseudo-)immutable, rather hashable. Changing something about them produces a new instance.
• IR objects implement custom __hash__ methods.
  • hash is computed from hash of contained objects
  • hash is computed upon initialization
• IR objects implement __eq__ method.
  • IR objects can be compared by simply using the == operator.
  • Equality is defined by same class and same hash. (Class could actually be part of the hash already.)
• Instead of copying IR objects from an object cache, identical objects could actually be references to the same object.
  • Changes to cached objects might need to produce new cache entries.

Advantages

• Simplifying syntax for comparing IR objects (a lot).
• Increase of performance in vectorization: Do not need to compare every little detail, just compare hash.
• Possibly simplify vectorization syntax a lot.
• Improve performance and memory usage of IR objects by referencing cached objects multiple times instead of copying.

Disadvantages

• Changes API of IR layer (only the part that faces the backend, though)
• Either backend needs to be changed to allow for immutable IR objects - or IR objects need to be implemented in such a way, that they can be changed naturally by silently producing new objects. That may cost performance, though.
• Possibly loss of performance due to complicated hashing (although this is only done once per object).

Edits: relaxed the immutability requirement to a hashable object that tries to prevent you from altering it.

Improving user friendliness and stability of PyRates

This is just a list of improvements that may be added to the user interface, documentation, testing suite etc.

Documentation

• FortranBackend
• CircuitIR methods
• readthedocs layout
• readthedocs how-to section
• known bugs, problems and unsupported usage of PyRates

pytests

• Tests for different types of edge delay conversions
• Tests for FortranBackend
• Tests for scipy solver
• Tests for multi-level circuit hierarchies
• Tests for new edge syntax in yaml files
• Tests for all methods on CircuitIR
• Tests for proper variable name translation into backend (including solving identical variable names)

convenience features

• variable setter method
• variable getter/setter methods that works with edge variables as well
• i/o methods for saving/loading different objects along the frontend-to-backend chain
• add more prints of ongoing processes to different stages of frontend-to-backend chain and make them mutable via verbose argument
• add more model templates for different descriptions of neural population and synaptic dynamics

Adding edge with any operator yields `KeyError` when running the simulation

I'd showcase on an example:
consider your Jansen-Rit Circuit implementation from model_templates:

...
nodes:
    PC: PC
    EIN: EIN
    IIN: IIN
  edges:
    - [PC/PRO/m_out, IIN/RPO_e/m_in, DummyEdge, {weight: 33.75}]
    - [PC/PRO/m_out, EIN/RPO_e/m_in, DummyEdge, {weight: 135.}]
    - [EIN/PRO/m_out, PC/RPO_e_pc/m_in, DummyEdge, {weight: 108.}]
    - [IIN/PRO/m_out, PC/RPO_i/m_in, DummyEdge, {weight: 33.75}]

Now, for the purpose of building more complex edges in the future, let just create linear coupling operator, assign it to the edge and scale the edge not via weight but rather via coupling value (the results should be the same):

LinearCoupling:
  base: OperatorTemplate
  equations: "V_in = c * V_out"
  variables:
    V_in:
      default: output
    c:
      default: 1.
    V_out:
      default: input

LinearEdge:
  base: EdgeTemplate
  operators:
    - LinearCoupling

and build the JRC with this edge template as:

 nodes:
    PC: PC
    EIN: EIN
    IIN: IIN
  edges:
    - [PC/PRO/m_out, IIN/RPO_e/m_in, LinearEdge, { LinearCoupling/c: 33.75, weight: 1. }]
    - [PC/PRO/m_out, EIN/RPO_e/m_in, LinearEdge, { LinearCoupling/c: 135., weight: 1. }]
    - [EIN/PRO/m_out, PC/RPO_e_pc/m_in, LinearEdge, { LinearCoupling/c: 108., weight: 1. }]
    - [IIN/PRO/m_out, PC/RPO_i/m_in, LinearEdge, { LinearCoupling/c: 33.75, weight: 1. }]

.

As you can see, the coupling parameter c takes the value of the edge weight and with weight being 1., this JRC circuit should have the same output as the original one from model template.

However, running this circuit (in the tutorial notebook) yields an error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~/.virtualenvs/sleepy3/lib/python3.7/site-packages/pyrates-0.8.0-py3.7.egg/pyrates/ir/circuit.py in get_node_var(self, key, apply_idx)
    843         try:
--> 844             return self[key]
    845         except KeyError:

~/.virtualenvs/sleepy3/lib/python3.7/site-packages/pyrates-0.8.0-py3.7.egg/pyrates/ir/abc.py in __getitem__(self, key)
     77             else:
---> 78                 raise e
     79 

~/.virtualenvs/sleepy3/lib/python3.7/site-packages/pyrates-0.8.0-py3.7.egg/pyrates/ir/abc.py in __getitem__(self, key)
     70             key = next(key_iter)
---> 71             item = self.getitem_from_iterator(key, key_iter)
     72             for key in key_iter:

~/.virtualenvs/sleepy3/lib/python3.7/site-packages/pyrates-0.8.0-py3.7.egg/pyrates/ir/circuit.py in getitem_from_iterator(self, key, key_iter)
    379         else:
--> 380             item = self.graph.nodes[key]["node"]
    381 

~/.virtualenvs/sleepy3/lib/python3.7/site-packages/networkx/classes/reportviews.py in __getitem__(self, n)
    177     def __getitem__(self, n):
--> 178         return self._nodes[n]
    179 

KeyError: 'PC'

During handling of the above exception, another exception occurred:

AttributeError                            Traceback (most recent call last)
<ipython-input-12-ddc332121c49> in <module>
      5                                      'V_EIN': 'EIN/RPO_e/PSP',
      6                                      'V_IIN': 'IIN/RPO_e/PSP'}, 
----> 7                             inputs={'PC/RPO_e_pc/u': ext_input})

~/.virtualenvs/sleepy3/lib/python3.7/site-packages/pyrates-0.8.0-py3.7.egg/pyrates/ir/circuit.py in run(self, simulation_time, inputs, outputs, sampling_step_size, out_dir, verbose, profile, **kwargs)
    966 
    967                 # extract respective output variables from the network and store their information
--> 968                 for var_key, var_info in self.get_node_var(val).items():
    969                     var_shape = tuple(var_info['var'].shape)
    970                     if var_shape in output_shapes:

~/.virtualenvs/sleepy3/lib/python3.7/site-packages/pyrates-0.8.0-py3.7.egg/pyrates/ir/circuit.py in get_node_var(self, key, apply_idx)
    849 
    850             # split original node keys
--> 851             node_keys = [key.split('/') for key in self.label_map]
    852 
    853             # remove all nodes from original node keys that are not referred to

~/.virtualenvs/sleepy3/lib/python3.7/site-packages/pyrates-0.8.0-py3.7.egg/pyrates/ir/circuit.py in <listcomp>(.0)
    849 
    850             # split original node keys
--> 851             node_keys = [key.split('/') for key in self.label_map]
    852 
    853             # remove all nodes from original node keys that are not referred to

AttributeError: 'tuple' object has no attribute 'split'

The label_map property of CircuitIR contains nodes' labels (PC, EIN, IIN) and edges in the form (PC, EIN, 0) and that is exactly the culprit - it fails on split() when encountering edge.

Allow additional (modulating) source nodes in edges

Follow-up to #14
Request by @Richert

Suggested syntax:

MyCircuit:
...
edges:
  -
    - [main_source: {...}, additional_source1: {...}, ...]
    - target/...
    - edge_template
    - {<attribute dict>}

which is identical to but clearer to read than:

MyCircuit:
...
edges:
  - [[main_source: {...}, additional_source1: {...}, ...], target/..., edge_template, {<attribute dict>}]

Within the networkx representation, the main source would be the actual "source" of the edge. Additional sources would only be known from within the edge.

Allow multiple source variables in edges

The goal is to enable dynamics in edges that need more than one input variable from the same node but do not make sense to move into nodes for logical or performance reasons. This feature was originally requested by @Richert .

Suggested YAML syntax:

MyCircuit: 
  ...
  edges:
    # single source variable as before
    # [source, target, edge_template, values]
    - [source_node/op/var, target_node/op/var, MyEdge, {weight: 1, ...}]
    # multiple source variables need to map source variables to edge variables
    # Option 1: keep syntax and repeat source node
    - [{source_node/op1/var1: edge_input_var1, 
         source_node/op2/var2: edge_input_var2, 
         ...}, 
       target_node/op/var, MyMultiEdge, {<variables>}]
    # Option 2: put source node in front and only reference ops/vars
    - [source_node: {op1/var1: edge_input_var1,
                               op2/var2: edge_input_var2, 
                               ...}, 
       target_node/op/var, MyMultiEdge, {<variables>}]

The previous syntax does not need to be changed, because single input variables can be uniquely identified. However, multiple input variables need to be mapped directly to source variables which complicates the syntax. Keeping the mapping logic of source -> target I prefer mapping <source_var>: <edge_var> over the alternative.

Option 1 is consistent with the previous syntax but repeating the source node is (1) more verbose which may be cumbersome and (2) suggests that more than one source node may exist which breaks the template. Option 2 makes it clear that only one source node exists but introduces a new syntax that is inconsistent with the single-variable case and the target variable syntax.

I prefer Option 2 for its clarity and comparative brevity.

The Python API should be updated the same way as the YAML syntax.

Note that variables used as input variables in the edge definition of a circuit also need to be defined as input variables in operators in the edge template. Furthermore the number of variables flagged as input in an edge template should always be the same as the assigned source variables in order to be consistent with current syntax. However, we could implement the system in such a way that a circuit can update variable definitions so the same template can be used with a different number of inputs as required.

Code changes that need to be made:

• EdgeIR: allow more than one input variable to exist in input property → remove error message
• CircuitIR (creation): Check if edge has multiple source variables and handle them correctly in
  • add_edge
  • add_edges_from
  • add_edges_from_matrix
• CircuitIR (vectorization): separate edges with multiple inputs into separate edges pointing to the same vectorized coupling node.

Note: The entire frontend code is surprisingly ignorant as to the content of source/target definitions and number of inputs.

Add to_file/from_file methods

To make it easier to work with large, densely connected networks, it would be helpful to have to_file(filename) and from_file(filename) methods implemented at the level of the CircuitIR class.

The file format could just be a simple pickle file in my opinion, since the purpose is not to just save a definition or configuration file, but to actually save progress that has been made due to the compilation/vectorization of the network, which can take very long for large systems.

This method should check, whether the circuit has already been compiled, i.e. whether the backend is already fully set up. In this case, pickling cannot be done right away, currently. The reason for this is that the backend creates and imports functions to simulate the network behavior, that cannot be pickled.

To solve this, I suggest to implement NumpyBackend.to_file() and .from_file methods as well that are called by the CircuitIR that write/load each of those functions to/from a python file. The strings for this are already generated in the backend.

Rework delay-coupling implementations

Currently, PyRates allows to define edge delays either within the framework of ordinary differential equations (ODEs) as convolutions with gamma kernels or within a delayed differential equation (DDE) framework via scalar delays.
While the former works just fine the way it is currently implemented, the latter comes with some restrictions.

It is currently implemented via the creation of local buffer variables for the delayed state variables.
This implementation only allows for the usage of the simple forward Euler algorithm for solving the DDE system, since higher-order methods would require PyRates to add interpolation schemes to the user-supplied equations that have not yet been implemented.
Also, the latter would require the definition of 2D interpolation methods for any of the different backends.
While that would technically be possible, it would be challenging to ensure that an appropriate interpolation algorithm is used, given the solver algorithm chosen by the user.
If this is not ensured, users might think that they solve their DDE system with a 5th order solver algorithm, even though the delay buffer interpolation is only performed at third order, for example.

I thus suggest that we implement an interface to commomly available DDE solvers instead, via the get_run_func method.
Such DDE solvers are available via Matlab's dde-biftool, Julia's DifferentialEquations.jl, or jiTCDDE in Python, for example.
, and they implement numerical solvers specifically dedicated to solving DDEs.
I suggest the following new feature for PyRates that would allow for such an interface:

• add a new backend function past(x,tau) that takes two positional arguments: a state variable x and a delay tau.
• given a choice of backend, translate this function call into the required syntax of a particular DDE solver
• if required, translate the state vector variable y from a 1D vector into a 2D matrix or callable (see requirements of the different DDE solvers listed above)

This implementation would probably work best, if we implement specific backends or backend keyword arguments for each DDE solver, since they all require slightly different implementations of the DDE system. This could be done similarly to the Auto-07p interface of the Fortran backend.

Help with model implementation / possible feature request

Dear PyRates developers,

I started playing with PyRates today and I must say - it is very intuitive and nicely written! Kudos to you!
I started implementing a model I need to play with (thalamocortical neural mass model for simulating sleep rhythms) due to Cona et al (link here).
Basically, it connects 6 nodes into a circuit (two thalamic and 4 cortical populations), where each node exhibit Jansen-Rit-like dynamics. I managed to create Nodes with correct Operators per each node in the model, however, I am struggling with the following:

1. some nodes (3 of them to be precise) have external input I_ext - not sure how I should implement. In the paper, the base equations is

and as far as I understand PyRates internals and Jansen-Rit example notebook the operators are taking care of z(t) and y(t). The whole summation over the network is taken care at the circuit level, so I am not sure how to add external input. The Jansen-Rit example has external input in the potential-to-rate operator, but to be honest, I am not sure whether that is applicable here.

2. as seen in the equations above, some of the edges (in particular the edges between cortical and thalamic nodes) experience a delay (to be precise a delay of 1ms). Is there a possibility to include delay in the edges? If not, might this be considered a feature request?

3. the peculiarity of the model is that self-connection of pyramidal neuron population (node) is subject to synaptic de-potentiation, i.e. when they fire at a high rate for a long time, they weaken as

I am not sure how to implement this temporal dependence of weights into the model. I noticed there is an `EdgeTemplate` that takes an operator (which would include this equation), however not sure how to set it up. Is this even possible? If not, might this be considered a feature request?

4. the thalamic populations of this model can fire in two modelities: tonacally or with bursts of action potentials. That means that their potential-to-rate operator has actually two parts:

where the first summand represents bursts, while the second is responsible for tonic firing. I would need to "record" as an output both rates separately. What would you suggest as the easiest way how to achieve this? I was thinking: set up three outputs, one for bursts, one for tonic firing and one as their sum, their sum will be considered input to other nodes and the former two just record as the outputs?

I would be very grateful for your help.
If I manage to implement the model correctly, I could contribute my code to PyRates docs to showcase some functionality not yet documented :)

Thanks a lot!
N.

Delay in edges does not work

Defining a delay in an edge yield an error of shape mismatch.

My implementation:

- [PYR/PyramidalCPO/V_out, TCR/ThalamicPRO/V_in, *LE, {*c : 1., weight: 1., delay: 0.001}]

hence I want to introduce a delay of 1ms (all my temporal variables are defined in seconds, hence delay=1e-3).
Both backends yield the following error:

ValueError: Invalid index shape. Operation has shape (12,), but received an index of length 2

using dt=1e-4. When using dt=1e-5 the shape mismatch becomes (102,) so I guess the discretization of delays or setting the correct length of history are the possible culprits.

I could try to fix it myself, if you'd point me in the right direction - the backend codebase is pretty large and on the first look, I couldn't see the associated methods.

Cheers,
N.

P.S. if I am not mistaken, also the associated test (tests/test_compute_graph.py) fails, I tested on MacOS python3.7 and Ubuntu18, python3.6

Unexpected behavior of a zero-weighted edge with delay

Dear PyRates developers,

after successfully reproducing Figure 5 from the David and Friston paper with @Richert's great help in this github issue, I ran into another problem while reproducing Figure 6 from the same paper.

I created two areas, each made up of the subpopulation model from Figure 5, added some variables to the RPO_e_pc Operator:

{
    'ext': {'default': 'input'},
    'n': {'default': 'variable'},
    'v_d1': {'default': 'variable'},
    'v_mean': {'default': 'variable'},
    'v_d2': {'default': 'variable'},
    'v_m2': {'default': 'variable'},
    'ext_variance': {'default': 'variable'},
    'k': {'default': 0.0}
}

added additional edges between the areas:

edges.append((f"PC0/PRO_pc/m_out", f"PC1/RPO_e_pc/ext", None, {'weight': 1.0, 'delay': 0.001, 'spread': 0.0005}))
edges.append((f"PC1/PRO_pc/m_out", f"PC0/RPO_e_pc/ext", None, {'weight': 1.0, 'delay': 0.001, 'spread': 0.0005}))

replaced the input equation in the RPO_e_pc Operator to weigh the inputs according to the connectivity parameters k (and the derived k*, as defined in the paper):

equations={'replace': {'m_in': '(m_in+(ext*(22.0*sqrt(2*k-k**2))/max(sqrt(ext_variance),1))+(u*(1-k)))'}}

and implemented an algorithm to iteratively calculate the sample variance of the external input ext (which is necessary for determining k*) in an incredibly hacky but apparently working way (any suggestions on how to do this differently, preferably without all those helper variables, is welcome)

equations = [
    'n = n + 1',
    'v_d1 = ext - v_mean',
    'v_mean = v_mean + v_d1 / (max(n,2) - 1)',
    'v_d2 = ext - v_mean',
    'v_m2 = v_m2 + v_d1*v_d2',
    'ext_variance = v_m2 / (max(n,3) - 2)',
]

This appeared to be working to some extent, as I received the expected synchronized signals (after giving the system a couple of seconds to converge) when connecting the areas bi-directionally with k1=k2=0.5, and increased the edge weights to values around 3.5 (I did not find any specific values for this anywhere, but weight 1.0 did not seem to work) .

However, when testing conditions with uni-directional (k1=0.5, k2=0.0) or zero-weighted (k1=k2=0.0) connections,
I received strange results, where areas that had their external input multiplied by zero started to behave differently than they would on their own. Also for some reason the order in which the areas were defined appeared to have an impact on the results.

I figured out that this behavior only emerged when using the delay parameter in the edge, and that it did not matter whether the variable that the connection was fed into was actually used.

In the table below, I plotted 6 conditions, all of which would be expected to produce results identical to each area being simulated on its own.
The additional equations and variables described above were removed to single out the problem, except for the edge-input ext,
which was however not replaced into the m_in calculation, but received a dummy equation ext = ext, removing (in my expectaion) its influence on the system.

The first column of the table shows how the two areas behave without any edge defined between them, producing the expected signals and spectra.
The bottom row shows the results when the areas are flipped, and the top right graph illustrates the main problem arising when using the edge delay (while the flipped condition in the bottom right corner works roughly as expected for some reason in this particular case).

..	Without Edges	Edges (weight 0.0, delay 0.000)	Edges (weight 0.0, delay 0.001)
Fast Area first
Slow Area first

In my understanding this is not the expected behavior of a zero-weighted edge with a delay.
Do I have any misconceptions on that or is it possible that I am using the edges incorrectly?

Thanks in advance for any insights on the matter.

error running jansenrit.py (jrc.compile) on MacOS

Hi, I have an issue when running the example jansenrit.py.
I'm using a conda environment and MacOS. Any ideas how to resolve the error?

Traceback (most recent call last):
  File "/Users/julian/Simulations/pyrates/jansenrit.py", line 113, in <module>
    jrc_compiled = jrc.compile(backend='numpy', step_size=1e-4, solver='scipy')
  File "/Users/julian/opt/miniconda3/envs/ADL/lib/python3.7/site-packages/pyrates/ir/circuit.py", line 1491, in compile
    squeeze=squeeze_vars)
  File "/Users/julian/opt/miniconda3/envs/ADL/lib/python3.7/site-packages/pyrates/backend/parser.py", line 800, in parse_equations
    instantaneous=instantaneous, **kwargs.copy())
TypeError: Can't instantiate abstract class ExpressionParser with abstract methods _generateDefaultName

Thank you

Removing parts of pyrates from sys.modules in the backend causes troubles in other parts of the pipeline

While debugging strange errors, I discovered that in lines 949-951 of pyrates.backend.numpy_backend the current collection known python modules is modified (sys.modules). In fact, all modules that contain the name of the current object instance (in this case "circuit") are removed. This state is kept across tests and thus leads to problems in other tests. However, I expect this to cause problems in productive code as well, because this is unexpected behaviour.

Commenting these lines out removed the problem and all (!) tests completed successfully. Is there a reason to keep these lines? @Richert

Trouble reproducing Jansen-Rit model dynamics in current version

Dear PyRates developers,
first off, I want to say that I really like PyRates for its versatility and simple usage- so thanks a lot for creating it!

That being said, I have been trying to reproduce findings from David and Friston (2003), but ran into trouble getting the expected results with the newest version of PyRates.

This Jupyter Notebook included with PyRates (Simulation 2: Jansen-Rit model sensitivity to changes in the synaptic timescales), appears to reproduce Fig. 4 from the paper quite well (as also stated in the PyRates paper), but when I try to run it with the current master (after adapting a few lines to make it execute without complaining - more on that below), it produces different output:

Figures

Output-Figure included in .ipynb-File	Result when running with current master
reproduces David and Friston (2003) Fig. 4	produces different data

Necessary changes

The changes required to execute the notebook-code with the current version were limited to the following:

results, params = grid_search(circuit_template="model_templates.jansen_rit.simple_jansenrit.JRC",
                             param_grid=params, param_map=param_map,
                             inputs={"PC/RPO_e_pc/u": inp}, outputs={"v": "PC/OBS/V"},
-                            dt=dt,
+                            step_size=dt, 
                             simulation_time=T, permute_grid=False, sampling_step_size=dts)

as well as

- freqs, power = fft(results[key], tmin=cut_off)
+ from pandas import DataFrame
+ freqs, power = fft(DataFrame(results[('v', key, 'PC')]), tmin=cut_off)

and, to make the graph more comparable (y-axis flip)

- tau_i = np.arange(1e-3, 63e-3, 5e-3)
+ tau_i = np.flip(np.arange(1e-3, 63e-3, 5e-3)) 

Steps taken so far

To investigate whether (and if so, at what point in the commit-history) the framework did actually produce the given plot, and to rule out that any updated dependencies may be the root of the problem, I checked out the newest commit that included the string notebook in its commit message.

This happend to be c97e8ba, and indeed (after applying the same np.flip adaption) it easily produces the following graph:

Result when running with c97e8ba
appears to be working correctly (minimal differences due to random input)

This must mean that somewhere since this commit, a change has been introduced that ever so slightly changes the behavior of the simple Jansen-Rit Model Circuit.
Interestingly enough, the properties of the signals that are produced within the usual ranges (τ_i = 20ms, τ_e = 10ms) do not change enough to cause noticable problems when reproducing other widely analysed features of the Jansen-Rit Model (e.g. the first part of the Jupyter Notebook).

As I have only recently started to dig into the deeper sections of the PyRates Code, I do not have any idea what kind of alteration may have caused such a subtle change, so I will probably start to try to single out the commit that introduces this change.

In the meantime, I would be glad if anybody could point me into the right direction.

Additional (possibly consequent?) problems

This change in behavior might also be responsible for the issue that made me aware of the problem in the first place:

When trying to reproduce Fig. 5 from the David and Friston paper, which required some minor changes to the model as well, I get the following result, which shows clear peaks in all conditions, while only the first and last are expected and all the other conditions should produce much more noisy, wider distributions:

Attempt to reproduce Fig. 5 from David and Friston

While it is entirely possible that there is an issue in my modification to the model that causes this, I have not yet gotten around to fix the compatibilty-errors that prevent me from running my newer implementations on the c97e8ba version and testing whether it works as expected there.

To keep this issue fairly simple, I wanted to put the focus on the more obvious problem first (which, on solution, might even fix this alltogether) - but again, I would be happy to get any insights on this as well, and will gladly provide and explain my code if anybody is interested in figuring this out with me.

Thanks in advance

Enable purely scalar network models by disabling vectorization

When calling the CircuitIR.compile() method, users can choose to either set vectorization=True or vectorization=False. However, currently this has no effect, since vectorization is performed in any case.

I suggest to implement the possibility to have a purely scalar network, without any vectorization being performed before the equations are parsed into the backend. This way, the compatibility with the FortranBackend will be improved, since it cannot handle non-scalar networks, as of now. Furthermore, this allows to create a PyAuto version of any network, since only scalar-valued differential equations are supported by auto-07p.

Find a way to count references on structurally unique IR objects for vectorization.

Vectorization could be sped up a lot by counting references on unique objects (like structurally identical nodes or operators).

Can we use Python's object reference counting for this? (Might not work due to missing back links)

If the above is not possible, it might make sense to save these references in a database-like structure (table) for quick access. Instead we could also just save back-references of the unique objects to the places they are referenced from as a list in the respective object.