Merge pull request #2 from ins-amu/develop

Develop

Author: Ziaeemehr

.github/workflows/tests.yml

@@ -0,0 +1,24 @@
+name: Test
+
+on: [push]
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: "3.10"
+    - name: Install dependencies
+      run: |
+        python -m pip install .
+
+    - name: Compile C++ code
+      run: |
+        cd vbi/models/cpp/_src
+        make
+    - name: Run tests
+      run: |
+        python -m pytest

pyproject.toml

@@ -39,7 +39,8 @@ dependencies = [
     "parameterized",
     "scikit-learn",
     "pycatch22",
-    "pytest"
+    "pytest",
+    "swig"
 ]
 classifiers = [
     "Programming Language :: Python :: 3",

vbi/models/cpp/mpr.py

@@ -1,8 +1,12 @@
 import numpy as np
 from typing import Union
 from copy import deepcopy
-from vbi.models.cpp._src.mpr_sde import MPR_sde as _MPR_sde
-from vbi.models.cpp._src.mpr_sde import BoldParams as _BoldParams
+
+try:
+    from vbi.models.cpp._src.mpr_sde import MPR_sde as _MPR_sde
+    from vbi.models.cpp._src.mpr_sde import BoldParams as _BoldParams
+except ImportError as e:
+    print(f"Could not import modules: {e}, probably C++ code is not compiled.")
 
 
 class MPR_sde:

vbi/models/cupy/bold.py

@@ -0,0 +1,111 @@
+import numpy as np
+
+
+class BoldStephan2008:
+
+    def __init__(self, par: dict = {}) -> None:
+
+        self._par = self.get_default_parameters()
+        self.valid_parameters = list(self._par.keys())
+        self.check_parameters(par)
+        self._par.update(par)
+
+        for key, value in self._par.items():
+            setattr(self, key, value)
+
+    def _prepare(self, nn, ns, xp, n_steps, bold_decimate):
+        s = xp.zeros((2, nn, ns), dtype=self.dtype)
+        f = xp.zeros((2, nn, ns), dtype=self.dtype)
+        ftilde = xp.zeros((2, nn, ns), dtype=self.dtype)
+        vtilde = xp.zeros((2, nn, ns), dtype=self.dtype)
+        qtilde = xp.zeros((2, nn, ns), dtype=self.dtype)
+        v = xp.zeros((2, nn, ns), dtype=self.dtype)
+        q = xp.zeros((2, nn, ns), dtype=self.dtype)
+        vv = np.zeros((n_steps // bold_decimate, nn, ns), dtype="f")
+        qq = np.zeros((n_steps // bold_decimate, nn, ns), dtype="f")
+        s[0] = 1
+        f[0] = 1
+        v[0] = 1
+        q[0] = 1
+        ftilde[0] = 0
+        vtilde[0] = 0
+        qtilde[0] = 0
+
+        return {
+            "s": s,
+            "f": f,
+            "ftilde": ftilde,
+            "vtilde": vtilde,
+            "qtilde": qtilde,
+            "v": v,
+            "q": q,
+            "vv": vv,
+            "qq": qq,
+        }
+
+    def check_parameters(self, par):
+        for key in par.keys():
+            if key not in self.valid_parameters:
+                raise ValueError(f"Invalid parameter {key:s} provided.")
+
+    def get_default_parameters(self):
+
+        theta0 = 41.0
+        Eo = 0.42
+        TE = 0.05
+        epsilon = 0.36
+        r0 = 26.0
+        k1 = 4.3 * theta0 * Eo * TE
+        k2 = epsilon * r0 * Eo * TE
+        k3 = 1 - epsilon
+
+        par = {
+            "kappa": 0.7,
+            "gamma": 0.5,
+            "tau": 1.0,
+            "alpha": 0.35,
+            "epsilon": epsilon,
+            "Eo": Eo,
+            "TE": TE,
+            "vo": 0.09,
+            "r0": r0,
+            "theta0": theta0,
+            "rtol": 1e-6,
+            "atol": 1e-9,
+            "k1": k1,
+            "k2": k2,
+            "k3": k3,
+        }
+        return par
+
+    def bold_step(self, r_in, s, f, ftilde, vtilde, qtilde, v, q, dt, P):
+
+        kappa, gamma, alpha, tau, Eo = P
+        ialpha = 1 / alpha
+
+        s[1] = s[0] + dt * (r_in - kappa * s[0] - gamma * (f[0] - 1))
+        f[0] = np.clip(f[0], 1, None)
+        ftilde[1] = ftilde[0] + dt * (s[0] / f[0])
+        fv = v[0] ** ialpha  # outflow
+        vtilde[1] = vtilde[0] + dt * ((f[0] - fv) / (tau * v[0]))
+        q[0] = np.clip(q[0], 0.01, None)
+        ff = (1 - (1 - Eo) ** (1 / f[0])) / Eo  # oxygen extraction
+        qtilde[1] = qtilde[0] + dt * ((f[0] * ff - fv * q[0] / v[0]) / (tau * q[0]))
+
+        f[1] = np.exp(ftilde[1])
+        v[1] = np.exp(vtilde[1])
+        q[1] = np.exp(qtilde[1])
+
+        f[0] = f[1]
+        s[0] = s[1]
+        ftilde[0] = ftilde[1]
+        vtilde[0] = vtilde[1]
+        qtilde[0] = qtilde[1]
+        v[0] = v[1]
+        q[0] = q[1]
+
+
+class BoldTVB:
+
+    def __init__(self):
+        pass

vbi/models/cupy/utils.py

@@ -1,10 +1,10 @@
 import numpy as np
-from numpy.matlib import repmat
+# from numpy.matlib import repmat
 
 try:
     import cupy as cp
 except:
-    pass
+    cp = None
 
 
 def get_module(engine="gpu"):
@@ -20,7 +20,7 @@ def get_module(engine="gpu"):
 
 def tohost(x):
     '''
-    move data to cpu
+    move data to cpu if it is on gpu
 
     Parameters
     ----------
@@ -32,7 +32,9 @@ def tohost(x):
     array
         data moved to cpu
     '''
-    return cp.asnumpy(x)
+    if cp is not None and isinstance(x, cp.ndarray):
+        return cp.asnumpy(x)
+    return x
 
 
 def todevice(x):
@@ -79,7 +81,7 @@ def repmat_vec(vec, ns, engine):
         repeated vector
 
     '''
-    vec = repmat(vec, ns, 1).T
+    vec = np.tile(vec, (ns, 1)).T
     vec = move_data(vec, engine)
     return vec
 

vbi/tests/test_module2.py

@@ -2,39 +2,43 @@
 import unittest
 import numpy as np
 import networkx as nx
-from vbi.models.cpp.mpr import MPR_sde
 
-seed = 2
-np.random.seed(seed)
-torch.manual_seed(seed)
+MPR_AVAILABLE = True
+try:
+    from vbi.models.cpp.mpr import MPR_sde
+except ImportError:
+    MPR_AVAILABLE = False
+
+
+SEED = 2
+np.random.seed(SEED)
+torch.manual_seed(SEED)
 
 nn = 3
 g = nx.complete_graph(nn)
-sc = nx.to_numpy_array(g)/ 10.0
+sc = nx.to_numpy_array(g) / 10.0
 
 
+@unittest.skipIf(not MPR_AVAILABLE, "vbi.models.cpp.mpr.MPR_sde module not available")
 class testMPRSDE(unittest.TestCase):
-    
+
     mpr = MPR_sde()
     p = mpr.get_default_parameters()
-    p['weights'] = sc
-    p['seed'] = seed
-    p['t_cut'] = 0.01 * 60 * 1000
-    p['t_end'] = 0.02 * 60 * 1000
-    
+    p["weights"] = sc
+    p["seed"] = SEED
+    p["t_cut"] = 0.01 * 60 * 1000
+    p["t_end"] = 0.02 * 60 * 1000
+
     def test_invalid_parameter_raises_value_error(self):
         invalid_params = {"invalid_param": 42}
         with self.assertRaises(ValueError):
             MPR_sde(par=invalid_params)
 
     def test_run(self):
-        
+
         control = {"G": 0.1, "eta": -4.7}
         mpr = MPR_sde(self.p)
         sol = mpr.run(par=control)
-        x = sol["x"]
-        t = sol["t"]
+        x = sol["bold_d"]
+        t = sol["bold_t"]
         self.assertEqual(x.shape[0], nn)
-        
-        
-    

vbi/tests/test_mpr_cpp.py

@@ -3,7 +3,12 @@
 import numpy as np
 import networkx as nx
 from copy import deepcopy
-from vbi.models.cupy.mpr import MPR_sde
+
+MPR_AVAILABLE = True
+try:
+    from vbi.models.cupy.mpr import MPR_sde
+except ImportError:
+    MPR_AVAILABLE = False
 
 
 seed = 2
@@ -15,6 +20,7 @@
 sc = nx.to_numpy_array(g) / 10.0
 
 
+@unittest.skipIf(not MPR_AVAILABLE, "vbi.models.cupy.mpr.MPR_sde module not available")
 class testMPRSDE(unittest.TestCase):
 
     mpr = MPR_sde()