Fixed slight bug with grid-search where circuit parametrizations were…

… overwritten.

documentation/Gast_2018_PyRates_benchmarks.py

@@ -11,15 +11,15 @@
 from numba import njit, config
 
 # threading configs
-config.THREADING_LAYER = 'tbb'
-os.environ["KMP_BLOCKTIME"] = '0'
-os.environ["KMP_SETTINGS"] = 'true'
-os.environ["KMP_AFFINITY"] = 'granularity=fine,verbose,compact,1,0'
-os.environ["OMP_NUM_THREADS"] = '1'
-tf.config.threading.set_inter_op_parallelism_threads(7)
-tf.config.threading.set_intra_op_parallelism_threads(1)
-tf.config.optimizer.set_jit(True)
-tf.config.experimental.set_synchronous_execution(False)
+# config.THREADING_LAYER = 'tbb'
+# os.environ["KMP_BLOCKTIME"] = '0'
+# os.environ["KMP_SETTINGS"] = 'true'
+# os.environ["KMP_AFFINITY"] = 'granularity=fine,verbose,compact,1,0'
+# os.environ["OMP_NUM_THREADS"] = '1'
+# tf.config.threading.set_inter_op_parallelism_threads(7)
+# tf.config.threading.set_intra_op_parallelism_threads(1)
+# tf.config.optimizer.set_jit(True)
+# tf.config.experimental.set_synchronous_execution(False)
 #tf.debugging.set_log_device_placement(True)
 
 
@@ -76,7 +76,7 @@ def benchmark(Ns, Ps, T, dt, init_kwargs, run_kwargs, disable_gpu=False):
             conns.index = [f'jrc_{idx}/PC/RPO_e_pc/m_in' for idx in range(n)]
 
             # define input
-            inp = 220 + np.asarray(np.random.randn(int(T / dt), n), dtype=np.float32) * 22.
+            inp = 220 + np.random.randn(int(T / dt), n) * 22.
 
             # set up template
             template = CircuitTemplate.from_yaml("model_templates.jansen_rit.simple_jansenrit.JRC")
@@ -107,8 +107,8 @@ def benchmark(Ns, Ps, T, dt, init_kwargs, run_kwargs, disable_gpu=False):
 dt = 1e-4                                       # integration step-size of the forward euler solver in s
 T = 1.0                                         # simulation time in s
 c = 1.                                          # global connection strength scaling
-N = np.round(2**np.arange(8))[::-1]             # network sizes, each of which will be run a benchmark for
-p = np.linspace(0.1, 1.0, 5)                    # global coupling probabilities to run benchmarks for
+N = np.round(2**np.arange(10))[::-1]            # network sizes, each of which will be run a benchmark for
+p = np.linspace(0.0, 1.0, 5)                    # global coupling probabilities to run benchmarks for
 use_gpu = False                                 # if false, benchmarks will be run on CPU
 n_reps = 1                                      # number of trials per benchmark
 

documentation/coupled_jrcs.py

@@ -1,5 +1,5 @@
 # pyrates imports
-from pyrates.utility import grid_search
+from pyrates.utility import grid_search, functional_connectivity
 from pyrates.utility.grid_search import linearize_grid
 
 # additional imports
@@ -10,29 +10,29 @@
 import matplotlib.pyplot as plt
 
 # threading configs
-config.THREADING_LAYER = 'tbb'
-os.environ["KMP_BLOCKTIME"] = '0'
-os.environ["KMP_SETTINGS"] = 'true'
-os.environ["KMP_AFFINITY"] = 'granularity=fine,verbose,compact,1,0'
-os.environ["OMP_NUM_THREADS"] = '2'
-tf.config.threading.set_inter_op_parallelism_threads(4)
-tf.config.threading.set_intra_op_parallelism_threads(2)
-tf.config.optimizer.set_jit(True)
-tf.config.experimental.set_synchronous_execution(False)
+# config.THREADING_LAYER = 'tbb'
+# os.environ["KMP_BLOCKTIME"] = '0'
+# os.environ["KMP_SETTINGS"] = 'true'
+# os.environ["KMP_AFFINITY"] = 'granularity=fine,verbose,compact,1,0'
+# os.environ["OMP_NUM_THREADS"] = '2'
+# tf.config.threading.set_inter_op_parallelism_threads(7)
+# tf.config.threading.set_intra_op_parallelism_threads(1)
+# tf.config.optimizer.set_jit(True)
+# tf.config.experimental.set_synchronous_execution(False)
 #tf.debugging.set_log_device_placement(True)
 
 # parameters
 ############
 
 # general parameters
 dt = 1e-4                                                       # integration step-size of the forward euler solver in s
-T = 1.0                                                         # overall simulation time in s
+T = 20.0                                                        # overall simulation time in s
 inp1 = np.random.uniform(120., 320., (int(T/dt), 1))            # white noise input to the pyramidal cells in Hz.
 inp2 = np.random.uniform(120., 320., (int(T/dt), 1))
 
 N = 10                                                          # grid-size
-C = np.linspace(0., 200., N)                                    # bi-directional connection strength
-D = np.linspace(0., 1e-2, N)                                    # bi-directional coupling delay
+C = np.linspace(0., 1., N)                                      # bi-directional connection strength
+D = np.linspace(0., 5e-2, N)                                    # bi-directional coupling delay
 
 # parameter grid
 params = {'C': C, 'D': D}
@@ -51,8 +51,7 @@
                                             "JRC2/PC/RPO_e_pc/u": np.asarray(inp2, dtype=np.float32)},
                                     outputs={"v": "all/PC/OBS/V"},
                                     dt=dt, simulation_time=T, permute_grid=True, sampling_step_size=1e-3,
-                                    init_kwargs={'solver': 'euler', 'backend': 'numpy'}, profile='t',
-                                    njit=True
+                                    init_kwargs={'solver': 'euler', 'backend': 'numpy'}, profile='t'
                                     )
 
 # tensorflow backend grid-search
@@ -65,5 +64,18 @@
 #                                     init_kwargs={'solver': 'euler', 'backend': 'tensorflow'}, profile='t',
 #                                     )
 
-results.plot()
+# coherence evaluation
+######################
+
+coherences = np.zeros((len(C), len(D)))
+for circuit_name in param_map.index:
+
+    results_tmp = results[circuit_name]
+    fc = functional_connectivity(results_tmp, metric='coh', fmin=8.0, fmax=12.0, faverage=True, tmin=1.0, verbose=False)
+    row = np.argwhere(C == param_map.loc[circuit_name, 'C'])
+    col = np.argwhere(D == param_map.loc[circuit_name, 'D'])
+    coherences[row, col] = fc[1, 0]
+
+plt.matshow(coherences)
+plt.tight_layout()
 plt.show()

pyrates/utility/data_analysis.py

@@ -213,7 +213,7 @@ def welch(data, tmin=0., tmax=None, **kwargs):
         return np.NaN, np.NaN
 
 
-def fft(data, tmin=0.):
+def fft(data, tmin=0., **kwargs):
     """
 
     Parameters
@@ -237,7 +237,7 @@ def fft(data, tmin=0.):
         data_tmp = data_tmp - np.mean(data_tmp)
 
         freqs = np.linspace(0, 1 / dt, n_two)
-        spec = np.fft.fft(data_tmp, n=n_two, axis=0)
+        spec = np.fft.fft(data_tmp, n=n_two, axis=0, **kwargs)
 
         # Cut of PSD and frequency arrays since its mirrored at N/2
         spec = np.abs(spec[:int(len(spec) / 2)])

pyrates/utility/grid_search.py

@@ -32,7 +32,8 @@
 # external imports
 import pandas as pd
 import numpy as np
-from typing import Optional, Union, Tuple
+from typing import Optional, Union
+from copy import deepcopy
 
 # system imports
 import os
@@ -41,7 +42,6 @@
 import time as t
 import glob
 import getpass
-from shutil import copy2
 from pathlib import Path
 from datetime import datetime
 from threading import Thread, currentThread, RLock
@@ -123,7 +123,7 @@ def grid_search(circuit_template: Union[CircuitTemplate, str], param_grid: Union
         new_params = {}
         for key in param_keys:
             new_params[key] = param_grid[key][idx]
-        circuit_tmp = circuit_template.apply()
+        circuit_tmp = deepcopy(circuit_template).apply()
         circuit_key = f'{circuit_tmp.label}_{idx}'
         circuit_tmp = adapt_circuit(circuit_tmp, new_params, param_map)
         circuit.add_circuit(circuit_key, circuit_tmp)