Merge pull request #71 from neurolib-dev/dev/exploration_tools

Exploration and Evolution: improvements

Author: caglorithm

README.md

@@ -30,13 +30,15 @@
 
 # What is neurolib?
 
-`neurolib` allows you to build, simulate, and optimize your own state-of-the-art whole-brain models. To simulate the neural activity of each brain area, the main implementation provides an advanced neural mass mean-field model of spiking adaptive exponential integrate-and-fire neurons (AdEx) called `aln`. Each brain area is represented by two populations of excitatory and inhibitory neurons. An extensive analysis and validation of the `aln` model can be found in our [paper](https://arxiv.org/abs/1906.00676) and its associated [github page](https://github.com/caglarcakan/stimulus_neural_populations).
+Please read the [gentle introduction](https://caglorithm.github.io/notebooks/neurolib-intro/) to `neurolib` for an overview of the basic functionality and some background information on the science behind whole-brain simulations.
+
+`neurolib` allows you to build, simulate, and optimize your own state-of-the-art whole-brain models. To simulate the neural activity of each brain area, the main implementation provides an advanced neural mass mean-field model of spiking adaptive exponential integrate-and-fire neurons (AdEx) called `aln`. Each brain area is represented by two populations of excitatory and inhibitory neurons. An extensive analysis and validation of the `aln` model can be found in our [paper](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007822) and its associated [github page](https://github.com/caglarcakan/stimulus_neural_populations).
 
 `neurolib` provides a simulation and optimization framework which allows you to easily implement your own neural mass model, simulate fMRI BOLD activity, analyse the results and fit your model to empirical data.
 
 Please reference the following paper if you use `neurolib` for your own research:
 
-**Reference:** Cakan, C., Obermayer, K. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation. PLOS Computational Biology ([ArXiv](https://arxiv.org/abs/1906.00676)).
+**Reference:** Cakan, C., Obermayer, K. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation. PLOS Computational Biology ([Link](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007822)).
 
 The figure below shows a schematic of how a brain network can be constructed:
 
@@ -145,7 +147,7 @@ aln.params['duration'] = 5*60*1000 # in ms, simulates for 5 minutes
 
 aln.run(bold=True)
 ```
-This can take several minutes to compute, since we are simulating 82 nodes for 5 minutes realtime. Note that we specified `bold=True` which simulates the BOLD model in parallel to the neuronal model. The resulting firing rates and BOLD functional connectivity looks like this:
+This can take several minutes to compute, since we are simulating 80 brain regions for 5 minutes realtime. Note that we specified `bold=True` which simulates the BOLD model in parallel to the neuronal model. The resulting firing rates and BOLD functional connectivity looks like this:
 <p align="center">
   <img src="https://github.com/neurolib-dev/neurolib/raw/master/resources/gw_simulated.png">
 </p>
@@ -200,13 +202,21 @@ We can plot the results to get something close to a bifurcation diagram!
 
 A detailed example is available as a [IPython Notebook](examples/example-2-evolutionary-optimization-minimal.ipynb). 
 
-`neurolib` also implements evolutionary parameter optimization, which works particularly well with brain networks. In an evolutionary algorithm, each simulation is represented as an individual and the parameters of the simulation, for example coupling strengths or noise level values, are represented as the genes of each individual. An individual is a part of a population. In each generation, individuals are evaluated and ranked according to a fitness criterion. For whole-brain network simulations, this could be the fit of the simulated activity to empirical data. Then, individuals with a high fitness value are `selected` as parents and `mate` to create offspring. These offspring undergo random `mutations` of their genes. After all offspring are evaluated, the best individuals of the population are selected to transition into the next generation. This process goes on for a given amount generations until a stopping criterion is reached. This could be a predefined maximum number of generations or when a large enough population with high fitness values is found.
+`neurolib` also implements evolutionary parameter optimization, which works particularly well with brain networks. In an evolutionary algorithm, each simulation is represented as an individual and the parameters of the simulation, for example coupling strengths or noise level values, are represented as the genes of each individual. An individual is a part of a population. In each generation, individuals are evaluated and ranked according to a fitness criterion. For whole-brain network simulations, this could be the fit of the simulated activity to empirical data. Then, individuals with a high fitness value are `selected` as parents and `mate` to create offspring. These offspring undergo random `mutations` of their genes. After all offspring are evaluated, the best individuals of the population are selected to transition into the next generation. This process goes on for a given amount generations until a stopping criterion is reached. This could be a predefined maximum number of generations or when a large enough population with high fitness values is found. 
+
+An example genealogy tree is shown below. You can see the evolution starting at the top and individuals reproducing generation by generation. The color indicates the fitness.
+
+<p align="center">
+  <img src="https://github.com/neurolib-dev/neurolib/raw/master/resources/evolution_tree.png", width="600">
+</p>
+
+`neurolib` makes it very easy to set up your own evolutionary optimization and everything else is handled under the hood. You can chose between two implemented evolutionary algorithms: `adaptive` is a gaussian mutation and rank selection algorithm with adaptive step size that ensures convergence (a schematic is shown in the image below). `nsga2` is an implementation of the popular multi-objective optimization algorithm by Deb et al. 2002. 
 
 <p align="center">
   <img src="https://github.com/neurolib-dev/neurolib/raw/master/resources/evolutionary-algorithm.png", width="600">
 </p>
 
-`neurolib` makes it very easy to set up your own evolutionary optimization and everything else is handled under the hood. Of course, if you like, you can dig deeper, define your own selection, mutation and mating operators. In the following demonstration, we will simply evaluate the fitness of each individual as the distance to the unit circle. After a couple of generations of mating, mutating and selecting, only individuals who are close to the circle should survive:
+Of course, if you like, you can dig deeper, define your own selection, mutation and mating operators. In the following demonstration, we will simply evaluate the fitness of each individual as the distance to the unit circle. After a couple of generations of mating, mutating and selecting, only individuals who are close to the circle should survive:
 
 ```python
 from neurolib.utils.parameterSpace import ParameterSpace

codecov.yml

@@ -12,5 +12,6 @@ coverage:
 	- "neurolib/models/thalamus/timeIntegration.py" 
 	- "neurolib/models/bold/timeIntegration.py" 
 	- "neurolib/utils/devUtils.py"
-  status:
-    patch: off
+	- "neurolib/utils/atlases.py"
+  status:	
+    patch: off

examples/example-1.2-brain-network-exploration.ipynb

@@ -546,7 +546,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "search.dfResults = eu.processExplorationResults(search.results, search.dfResults, model=model, ds=ds, bold_transient=10000)"
+    "eu.processExplorationResults(search, model=model, ds=ds, bold_transient=10000)"
    ]
   },
   {
@@ -842,7 +842,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.4"
   }
  },
  "nbformat": 4,

examples/example-2-evolutionary-optimization-minimal.ipynb

@@ -1,5 +1,16 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Evolutionary optimization of a whole-brain model\n",
+    "\n",
+    "This notebook provides an example for the use of the evolutionary optimization framework built-in to the library. Under the hood, the implementation of the evolutionary algorithm is powered by `deap` and `pypet` cares about the parallelization and storage of the simulation data for us.\n",
+    "\n",
+    "We want to optimize a whole-brain network that should produce simulated BOLD activity (fMRI data) that is similar to the empirical dataset. We measure the fitness of each simulation by computing the `func.matrix_correlation` of the functional connectivity `func.fc(model.BOLD.BOLD)` to the empirical data `ds.FCs`. The ones that are closest to the empirical data get a higher fitness and have a higher chance of reproducing and survival. "
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -44,28 +55,42 @@
     "plt.rcParams['image.cmap'] = 'plasma'"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We create a brain network model using the empirical dataset `ds`:"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
-    "aln = ALNModel(Cmat = ds.Cmat, Dmat = ds.Dmat) # simulates the whole-brain model in 10s chunks by default if bold == True\n",
+    "model = ALNModel(Cmat = ds.Cmat, Dmat = ds.Dmat) # simulates the whole-brain model in 10s chunks by default if bold == True\n",
     "# Resting state fits\n",
-    "aln.params['mue_ext_mean'] = 1.57\n",
-    "aln.params['mui_ext_mean'] = 1.6\n",
-    "aln.params['sigma_ou'] = 0.09\n",
-    "aln.params['b'] = 5.0\n",
-    "aln.params['signalV'] = 2\n",
-    "aln.params['dt'] = 0.2\n",
-    "aln.params['duration'] = 0.2 * 60 * 1000 #ms\n",
+    "model.params['mue_ext_mean'] = 1.57\n",
+    "model.params['mui_ext_mean'] = 1.6\n",
+    "model.params['sigma_ou'] = 0.09\n",
+    "model.params['b'] = 5.0\n",
+    "model.params['signalV'] = 2\n",
+    "model.params['dt'] = 0.2\n",
+    "model.params['duration'] = 0.2 * 60 * 1000 #ms\n",
     "# testing: aln.params['duration'] = 0.2 * 60 * 1000 #ms\n",
     "# real: aln.params['duration'] = 1.0 * 60 * 1000 #ms"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our evaluation function will do the following: first it will simulate the model for a short time to see whether there is any sufficient activity. This speeds up the evolution considerably, since large regions of the state space show almost no neuronal activity. Only then do we simulate the model for the full duration and compute the fitness using the empirical dataset."
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -79,26 +104,16 @@
     "    \n",
     "    # Stage 1 : simulate for a few seconds to see if there is any activity\n",
     "    # ---------------------------------------\n",
-    "    model.params['dt'] = 0.1\n",
     "    model.params['duration'] = 3*1000.\n",
     "    model.run()\n",
     "    \n",
     "    # check if stage 1 was successful\n",
-    "    if np.max(aln.rates_exc[:, aln.t > 500]) > 300 or np.max(aln.rates_exc[:, aln.t > 500]) < 10:\n",
-    "        return invalid_result, {}\n",
-    "    \n",
-    "    # Stage 2: simulate BOLD for a few seconds to see if it moves\n",
-    "    # ---------------------------------------\n",
-    "    model.params['dt'] = 0.2\n",
-    "    model.params['duration'] = 20*1000.\n",
-    "    model.run(bold = True)\n",
-    "    \n",
-    "    if np.std(aln.BOLD.BOLD[:, 5:10]) < 0.001:\n",
+    "    if np.max(model.output[:, model.t > 500]) > 160 or np.max(model.output[:, model.t > 500]) < 10:\n",
     "        return invalid_result, {}\n",
+    "\n",
     "    \n",
-    "    # Stage 3: full and final simulation\n",
+    "    # Stage 2: full and final simulation\n",
     "    # ---------------------------------------\n",
-    "    model.params['dt'] = 0.2\n",
     "    model.params['duration'] = defaultDuration\n",
     "    model.run(chunkwise=True, bold = True)\n",
     "    \n",
@@ -118,6 +133,13 @@
     "    return fitness_tuple, {}"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We specify the parameter space that we want to search."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 14,
@@ -128,15 +150,29 @@
     "                      [[0.0, 3.0], [0.0, 3.0], [0.0, 100.0], [0.0, 0.3], [0.0, 500.0], [0.0, 400.0]])"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that we chose `algorithm='nsga2'` when we create the `Evolution`. This will use the multi-objective optimization algorithm by Deb et al. 2002. Although we have only one objective here (namely the FC fit), we could in principle add more objectives, like the `FCD` matrix fit or other objectives. For this, we would have to add these values to the fitness in the evaluation function above and add more `weights` in the definition of the `Evolution`. We can use positive weights for that objective to be maximized and negative ones for minimization. Please refer to the [DEAP documentation](https://deap.readthedocs.io/) for more information."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "evolution = Evolution(evaluateSimulation, pars, weightList = [1.0] * len(ds.BOLDs), model = aln, POP_INIT_SIZE=4, POP_SIZE = 4, NGEN=2, filename=\"example-2.2.hdf\")\n",
-    "#testing: evolution = Evolution(evaluateSimulation, pars, weightList = [1.0] * len(ds.BOLDs), model = aln, POP_INIT_SIZE=4, POP_SIZE = 4, NGEN=2)\n",
-    "# real: evolution = Evolution(evaluateSimulation, pars, weightList = [1.0] * len(ds.BOLDs), model = aln, POP_INIT_SIZE=1600, POP_SIZE = 160, NGEN=100)"
+    "evolution = Evolution(evaluateSimulation, pars, algorithm = 'nsga2', weightList = [1.0] * len(ds.BOLDs), model = model, POP_INIT_SIZE=4, POP_SIZE = 4, NGEN=2, filename=\"example-2.2.hdf\")\n",
+    "#testing: evolution = Evolution(evaluateSimulation, pars, algorithm = 'nsga2', weightList = [1.0] * len(ds.BOLDs), model = model, POP_INIT_SIZE=4, POP_SIZE = 4, NGEN=2)\n",
+    "# real: evolution = Evolution(evaluateSimulation, pars, algorithm = 'nsga2', weightList = [1.0] * len(ds.BOLDs), model = model, POP_INIT_SIZE=1600, POP_SIZE = 160, NGEN=100)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "That's it, we can run the evolution now."
    ]
   },
   {
@@ -148,13 +184,27 @@
     "evolution.run(verbose = False)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We could now save the full evolution object for later analysis using `evolution.saveEvolution()`."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "# Analysis"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `info()` method gives us a useful overview of the evolution, like a summary of the evolution parameters, the statistics of the population and also scatterplots of the individuals in our search space."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 17,
@@ -297,18 +347,6 @@
     "plt.xlabel(\"Generation #\")\n",
     "plt.ylabel(\"Score\")"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 26,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# if we were lazy, we could just dump the whole evolution object into a file:\n",
-    "# import dill\n",
-    "# fname = os.path.join(\"data/\", \"evolution-\" + evolution.trajectoryName + \".dill\")\n",
-    "# dill.dump(evolution, open(\"~evolution.dill\", \"wb\"))"
-   ]
   }
  ],
  "metadata": {
@@ -327,7 +365,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

examples/example-2.1-evolutionary-optimization-aln.ipynb

@@ -31,7 +31,7 @@ def simulateBOLD(Z, dt, voxelCounts, X=None, F=None, Q=None, V=None):
     if "voxelCounts" not in globals():
         voxelCounts = np.ones((N,))
 
-    # Balloon-Windkessel model parameters (Deco 2013, Friston 2003):
+    # Balloon-Windkessel model parameters (Deco 2013, Friston 2000):
     # Note: the distribution of each Balloon-Windkessel models parameters are given per voxel
     # Since we usually average the empirical fMRI of each voxel for a given area, the standard
     # deviation of the gaussian distribution should be divided by the number of voxels in each area

examples/example-2.2-evolution-brain-network-aln-resting-state-fit.ipynb

@@ -92,7 +92,7 @@ def simulateBold(self, t, variables, append=False):
                         self.setOutput("BOLD.BOLD", BOLD)
                     else:
                         logging.warn(
-                            f"Will not simulate BOLD if output {bold_input.shape[1]} not at least of duration {self.boldModel.samplingRate_NDt*self.params['dt']}"
+                            f"Will not simulate BOLD if output {bold_input.shape[1]*self.params['dt']} not at least of duration {self.boldModel.samplingRate_NDt*self.params['dt']}"
                         )
         else:
             logging.warn("BOLD model not initialized, not simulating BOLD. Use `run(bold=True)`")

neurolib/models/bold/timeIntegration.py

@@ -13,6 +13,14 @@
 #     return ParametersInterval(*(individual[: len(paramInterval)]))._asdict().copy()
 
 
+def randomParameters(paramInterval):
+    """
+    Generate a sequence of random parameters from a ParamsInterval using a uniform distribution.
+    Format: [mean_par1, mean_par2, ...]
+    """
+    params = [np.random.uniform(*pI) for pI in paramInterval]
+    return params
+
 def randomParametersAdaptive(paramInterval):
     """
     Generate a sequence of random parameters from a ParamsInterval using a uniform distribution.
@@ -28,7 +36,7 @@ def randomParametersAdaptive(paramInterval):
     return params
 
 
-def mutateUntilValid(pop, paramInterval, toolbox, maxTries=0):
+def mutateUntilValid(pop, paramInterval, toolbox, MUTATE_P={}, maxTries=100):
     """Checks the validity of new individuals' parameter. If they are invalid 
     (for example if they are out of the predefined paramter space bounds), 
     mutate the individual, until valid.
@@ -42,7 +50,7 @@ def mutateUntilValid(pop, paramInterval, toolbox, maxTries=0):
     for i, ind in enumerate(pop):
 
         ind_bak = copy.copy(ind)
-        toolbox.mutate(pop[i])
+        toolbox.mutate(pop[i], **MUTATE_P)
 
         nMutations = 0
         while not checkParamValidity(pop[i], paramInterval) and nMutations < maxTries:
@@ -329,4 +337,3 @@ def gaussianAdaptiveMutation_nStepSizes(individual, gamma_gl=None, gamma=None):
     individual[:] = newParams + newSigmas
 
     return (individual,)
-