From 8d9ab5ae69ccfdf6a75579d9f6a59ef5cd8d47d6 Mon Sep 17 00:00:00 2001
From: Remi Gau <remi_gau@hotmail.com>
Date: Sun, 18 Aug 2024 11:32:41 +0200
Subject: [PATCH] Add 'tutorial' (#23)

* add tutorial

* apply precommit
---
 .github/workflows/run_examples.yml |  2 +-
 examples/tutorial.py               | 81 ++++++++++++++++++++++++++++++
 neurodesign/classes.py             |  4 +-
 neurodesign/report.py              |  8 ++-
 4 files changed, 90 insertions(+), 5 deletions(-)
 create mode 100644 examples/tutorial.py

diff --git a/.github/workflows/run_examples.yml b/.github/workflows/run_examples.yml
index 801d40d..7652754 100644
--- a/.github/workflows/run_examples.yml
+++ b/.github/workflows/run_examples.yml
@@ -22,7 +22,7 @@ jobs:
         strategy:
             fail-fast: false
             matrix:
-                example: [compare_designs.py, JSS_example.py, optimisation.py]
+                example: [compare_designs.py, JSS_example.py, optimisation.py, tutorial.py]
         steps:
         -   name: Checkout
             uses: actions/checkout@v4
diff --git a/examples/tutorial.py b/examples/tutorial.py
new file mode 100644
index 0000000..7c542ff
--- /dev/null
+++ b/examples/tutorial.py
@@ -0,0 +1,81 @@
+# Adapted from https://lukas-snoek.com/NI-edu/fMRI-introduction/week_3/neurodesign.html
+
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+from rich import print
+
+import neurodesign
+
+# Neurodesigninternally parallelizes some computations using multithreading,
+# which is a massive burden on the CPU. So let's limit the number of threads
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+os.environ["MKL_NUM_THREADS"] = "1"
+
+
+TR = 1.6
+rho = 0.6
+
+n_stimuli = 2
+stim_duration = 1
+duration = 5 * 60
+P = [0.5, 0.5]
+t_pre = 0.1
+t_post = 0.1
+
+ITImodel = "uniform"
+ITImin = 2
+ITImax = 4
+
+C = np.array([[1, -1], [-1, 1]])
+
+#  %%
+exp = neurodesign.Experiment(
+    TR=TR,
+    rho=rho,
+    n_stimuli=n_stimuli,
+    stim_duration=stim_duration,
+    P=P,
+    duration=duration,
+    t_pre=t_pre,
+    t_post=t_post,
+    ITImodel=ITImodel,
+    ITImin=ITImin,
+    ITImax=ITImax,
+    C=C,
+)
+
+#  %%
+weights = [0, 1, 0, 0]  # order: Fe, Fd, Ff, Fc
+outdes = 10
+
+opt = neurodesign.Optimisation(
+    # we have to give our previously created `exp` object to this class as well
+    experiment=exp,
+    weights=weights,
+    preruncycles=10,
+    cycles=1,
+    seed=2,
+    outdes=outdes,
+)
+
+opt.optimise()
+opt.evaluate()
+
+print(f"Onsets: {opt.bestdesign.onsets}")
+print(f"Order: {opt.bestdesign.order}")
+
+Xconv = opt.bestdesign.Xconv
+
+plt.figure(figsize=(15, 5))
+plt.plot(Xconv)
+for ons, cond in zip(opt.bestdesign.onsets, opt.bestdesign.order):
+    c = "tab:blue" if cond == 0 else "tab:orange"
+    plt.plot([ons, ons], [0.35, 0.37], c=c, lw=2)
+
+plt.legend(["Faces", "Houses"])
+plt.grid()
+plt.xlim(0, Xconv.shape[0])
+plt.show()
diff --git a/neurodesign/classes.py b/neurodesign/classes.py
index 8120af2..eb49d1f 100644
--- a/neurodesign/classes.py
+++ b/neurodesign/classes.py
@@ -482,7 +482,7 @@ def __init__(
         if not np.isclose(self.TR % self.resolution, 0):
             self.resolution = _find_new_resolution(self.TR, self.resolution)
             warnings.warn(
-                "the resolution is adjusted to be a multiple of the TR."
+                "the resolution is adjusted to be a multiple of the TR. "
                 f"New resolution: {self.resolution}"
             )
 
@@ -773,7 +773,7 @@ def add_new_designs(self, weights=None, R=None):
         :type experiment: experiment
         :param weights: weights for efficiency calculation.
         :type weights: list of floats, summing to 1
-        :param seed: The seed for ramdom processes.
+        :param seed: The seed for random processes.
         :type seed: integer or None
         """
         # weights
diff --git a/neurodesign/report.py b/neurodesign/report.py
index b1f4717..7e30aba 100644
--- a/neurodesign/report.py
+++ b/neurodesign/report.py
@@ -3,6 +3,7 @@
 import time
 from pathlib import Path
 
+import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib import gridspec
@@ -22,11 +23,12 @@
     Table,
 )
 
-plt.switch_backend("agg")
-
 
 def make_report(population, outfile: str | Path = "NeuroDesign.pdf"):
     """Create a report of a finished design optimisation."""
+    backend = matplotlib.get_backend()
+    plt.switch_backend("agg")
+
     if not isinstance(population.cov, np.ndarray):
         population.evaluate()
 
@@ -172,6 +174,8 @@ def make_report(population, outfile: str | Path = "NeuroDesign.pdf"):
 
     doc.build(story)
 
+    plt.switch_backend(backend)
+
 
 def _form_xo_reader(imgdata):
     (page,) = PdfReader(imgdata).pages