Merge pull request #21 from neurolib-dev/fix/better_documentation

Fix/better documentation
neurolib-dev · Feb 3, 2020 · cf3593e · cf3593e
2 parents 30ce50e + b0b9d85
commit cf3593e
Show file tree

Hide file tree

Showing 5 changed files with 202 additions and 39 deletions.
diff --git a/README.md b/README.md
@@ -3,15 +3,30 @@
 # Neurolib
 *Easy whole-brain neural mass modeling* 👩‍🔬💻🧠
 
-Neurolib allows you to easily create your own state-of-the-art whole-brain models. The main implementation is a neural mass firing rate model of spiking adaptive exponential integrate-and-fire neurons (AdEx) called `aln` which consists of two populations of excitatory and a inhibitory neurons. An extensive analysis of the model can be found in the paper *Biophysically grounded mean-field models of neural populations under electrical stimulation*, Cakan et al. 2020 ([ArXiv](https://arxiv.org/abs/1906.00676)), and its associated [github page](https://github.com/caglarcakan/stimulus_neural_populations).
+`Neurolib` allows you to easily create your own state-of-the-art whole-brain models. The main implementation is a neural mass firing rate model of spiking adaptive exponential integrate-and-fire neurons (AdEx) called `aln` which consists of two populations of excitatory and inhibitory neurons. 
+
+An extensive analysis of this model can be found in our paper and its associated [github page](https://github.com/caglarcakan/stimulus_neural_populations).
+
+**Reference:** Cakan, C., Obermayer, K. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation ([ArXiv](https://arxiv.org/abs/1906.00676)).
 
 <p align="center">
   <img src="resources/pipeline.png" width="700">
 </p>
 
-### Whole-brain modeling
+<p style="text-align: center;"> 
+<br>
+<b> Quick examples:
+<a href="#single-node">Single node simulation</a> 
+<a href="#whole-brain-network">Whole-brain network</a>
+<a href="#parameter-exploration">Parameter exploration</a>
+<a href="#evolutionary-optimization">Evolutionary optimization</a>
+</b>
+<br>
+</p>
+
+## Whole-brain modeling
 
-In combination with structural brain data, for example from diffusion tensor imaging (DTI) [tractography](https://en.wikipedia.org/wiki/Tractography), and resting state [BOLD](https://en.wikipedia.org/wiki/Blood-oxygen-level-dependent_imaging) data from magnetic resonance imaging (rs-fMRI), a network model of a whole brain can be created. Structural connectivity matrices from DTI tractography define 1) the connection strengths between areas, represented for example by the number of axonal fibers between each two brain areas and 2) the signal transmission delays measured from the length of the axonal fibers. 
+In combination with structural brain data, for example from diffusion tensor imaging (DTI) [tractography](https://en.wikipedia.org/wiki/Tractography), and functional resting state [BOLD](https://en.wikipedia.org/wiki/Blood-oxygen-level-dependent_imaging) time series data from magnetic resonance imaging (rs-fMRI), a network model of a whole brain can be created. Structural connectivity matrices from DTI tractography define 1) the connection strengths between areas, represented for example by the number of axonal fibers between each two brain areas and 2) the signal transmission delays measured from the length of the axonal fibers. 
 
 The resulting whole-brain model consists of interconnected brain areas, with each brain area having their internal neural dynamics. The neural activity is used to simulate BOLD activity using the Balloon-Windkessel model. The resulting simulated [resting state functional connectivity](https://en.wikipedia.org/wiki/Resting_state_fMRI#Functional) can then be used to fit the model to empirical functional brain data. 
 
@@ -21,7 +36,7 @@ Below is an animation in which the neural activity from such a model is plotted
   <img src="resources/brain_slow_waves_small.gif">
 </p>
 
-## Installation
+# Installation
 The easiest way to get going is to install the pypi release of `neurolib` using
 
 ```
@@ -33,12 +48,16 @@ Alternatively, you can also clone this repository and install all dependencies w
 git clone https://github.com/neurolib-dev/neurolib.git
 cd neurolib/
 pip install -r requirements.txt
+pip install .
 ```
 
-## Usage
-Example iPython notebooks on how to use the library can be found in the `./examples/` directory. A basic overview is given here. 
+# Usage
+Example [IPython Notebooks](examples/) on how to use the library can be found in the `./examples/` directory, don't forget to check them out! A basic overview of the functionality that `Neurolib` provides is given here. 
+
+## Single node
+
+A detailed example is available as a [IPython Notebook](examples/example-0-aln-minimal.ipynb). 
 
-### Single node
 To create a single `aln` model with the default parameters, simply run
 
 ```python
@@ -61,7 +80,9 @@ plt.plot(aln.t, aln.rates_exc.T)
   <img src="resources/single_timeseries.png">
 </p>
 
-### Whole-brain network
+## Whole-brain network
+
+A detailed example is available as a [IPython Notebook](examples/example-0-aln-minimal.ipynb). 
 
 To simulate a whole-brain network model, first we need to load a DTI and a resting-state fMRI dataset (an example dataset called `gw` is provided in the `neurolib/data/datasets/` directory).
 
@@ -89,7 +110,9 @@ This can take several minutes to compute, since we are simulating 90 nodes for 5
   <img src="resources/gw_simulated.png">
 </p>
 
-The quality of the fit of this simulation to the functional empirical data can now be computed per subject or for the whole group on average:
+The quality of the fit of this simulation can be computed by correlating the simulated functional connectivity matrix above to the empirical resting-state functional connectivity for each subject. This gives us an estimate of how well the model reproduces inter-areal BOLD correlations. As a rule of thumb, a value above 0.5 is considered good. 
+
+We can compute this by using the builtin functions `func.fc` to calculate a functional connectivity from `N` (`N` = number of regions) time series and and `fund.matrix_correlation` to compare this to the empirical data.
 
 ```python
 scores = []
@@ -113,4 +136,78 @@ Subject 9: 0.54
 Subject 10: 0.49
 Mean simulated FC to empirical FC correlation: 0.57
 ```
+## Parameter exploration
+A detailed example is available as a [IPython Notebook](examples/example-1-aln-parameter-exploration.ipynb). 
+
+Whenever you work with a model, it is of great importance to know what kind of dynamics it exhibits given a certain set of parameter. For this, it is useful to get an overview of the state space of a given model. For example in the case of `aln`, the dynamics depends a lot on the mean inputs to the excitatory and the inhibitory population. `Neurolib` makes it very easy to quickly explore parameter spaces of a given model:
+
+```python
+# create model
+aln = ALNModel()
+# define the parameter space to explore
+parameters = {'mue_ext_mean' : np.linspace(0, 3, 11).tolist(),
+              'mui_ext_mean' : np.linspace(0, 3, 11).tolist()}       
+
+# define exploration              
+search = BoxSearch(aln, parameters)
+search.initializeExploration()
+search.run()                
+```
+That's it!. You can now use the builtin functionality to read the simulation results from disk and do your analysis:
+
+```python
+search.loadResults()
+
+# calculate maximum firing rate for each parameter
+for i in search.dfResults.index:
+    search.dfResults.loc[i, 'max_r'] = np.max(search.results[i]['rates_exc'][:, -int(1000/aln.params['dt']):])
+```
+We can plot the results to get something close to a bifurcation diagram!
+
+<p align="center">
+  <img src="resources/exploration_aln.png">
+</p>
+
+## Evolutionary optimization
+
+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. An individual is a part of a population. In each generation, individuals are evaluated according to a fitness criterion. Afterwards, all individuals of a population which have a high fitness value are able to mate and create offspring. This offspring undergoes some random changes, which is often called mutation. Then all offsprings are evaluated and the best part of the population is selected. This process goes on for a given amount generations until a good population with high-fitness individuals is found.
+
+`Neurolib` makes it very easy to set up your own evolutionary optimization and everything else is handled under the hood. In this example, we will simply calculate the example of each individual as the distance to the unit circle. 
+
+```python
+from neurolib.utils.parameterSpace import ParameterSpace
+from neurolib.optimize.evolution import Evolution
+
+def optimize_me(traj):
+    ind = evolution.getIndividualFromTraj(traj)
+
+    # let's make a circle
+    fitness_result = abs((ind.x**2 + ind.y**2) - 1)
+
+    # gather results
+    fitness_tuple = (fitness_result ,)
+    result_dict = {"result" : [fitness_result]}
+
+    return fitness_tuple, result_dict
+
+# we define a parameter space and its boundaries
+pars = ParameterSpace(['x', 'y'], [[-5.0, 5.0], [-5.0, 5.0]])
+
+# initialize the evolution and go
+evolution = Evolution(optimize_me, pars, weightList = [-1.0], POP_INIT_SIZE= 100, POP_SIZE = 50, NGEN=10)
+evolution.run()    
+```
+
+That's all! Now you can check the results!
+
+```python
+evolution.loadResults()
+evolution.info(plot=True)
+```
+This will give you a summary of the last generation and plot a distribution of the individuals (and their parameters). As you can see in the parameter space cross sections below, all remaining individuals lie on a circle.
 
+<p align="center">
+  <img src="resources/evolution_minimal.png">
+</p>
diff --git a/neurolib/models/aln/model.py b/neurolib/models/aln/model.py
@@ -20,7 +20,16 @@ class ALNModel(Model):
     modelOutputNames = ["rates_exc", "rates_inh"]
 
     def __init__(
-        self, params=None, Cmat=[], Dmat=[], lookupTableFileName=None, seed=None, simulateChunkwise=False, chunkSize=10000, simulateBOLD=False, saveAllActivity=False,
+        self,
+        params=None,
+        Cmat=[],
+        Dmat=[],
+        lookupTableFileName=None,
+        seed=None,
+        simulateChunkwise=False,
+        chunkSize=10000,
+        simulateBOLD=False,
+        saveAllActivity=False,
     ):
         """
         :param params: parameter dictionary of the model
@@ -46,29 +55,52 @@ def __init__(
         self.simulateBOLD = simulateBOLD  # BOLD
         if simulateBOLD:
             self.simulateChunkwise = True  # Override this setting if BOLD is simulated!
-        self.chunkSize = chunkSize  # Size of integration chunks in chunkwise integration in case of simulateBOLD == True
-        self.saveAllActivity = saveAllActivity  # Save data of all chunks? Can be very memory demanding if simulations are long or large
+        self.chunkSize = (
+            chunkSize  # Size of integration chunks in chunkwise integration in case of simulateBOLD == True
+        )
+        self.saveAllActivity = (
+            saveAllActivity  # Save data of all chunks? Can be very memory demanding if simulations are long or large
+        )
 
         # load default parameters if none were given
         if params == None:
-            self.params = dp.loadDefaultParams(Cmat=self.Cmat, Dmat=self.Dmat, lookupTableFileName=self.lookupTableFileName, seed=self.seed,)
+            self.params = dp.loadDefaultParams(
+                Cmat=self.Cmat, Dmat=self.Dmat, lookupTableFileName=self.lookupTableFileName, seed=self.seed,
+            )
         else:
             self.params = params
 
     def run(self):
         """
-        Runs an aLN mean-field model simulation
+        Runs an aLN mean-field model simulation.
         """
 
         if self.simulateChunkwise:
-            t_BOLD, BOLD, return_tuple = cw.chunkwiseTimeIntAndBOLD(self.params, self.chunkSize, self.simulateBOLD, self.saveAllActivity)
+            t_BOLD, BOLD, return_tuple = cw.chunkwiseTimeIntAndBOLD(
+                self.params, self.chunkSize, self.simulateBOLD, self.saveAllActivity
+            )
             rates_exc, rates_inh, t, mufe, mufi, IA, seem, seim, siem, siim, seev, seiv, siev, siiv = return_tuple
             self.t_BOLD = t_BOLD
             self.BOLD = BOLD
             Model.setOutput(self, "BOLD.t", t_BOLD)
             Model.setOutput(self, "BOLD.BOLD", BOLD)
         else:
-            rates_exc, rates_inh, t, mufe, mufi, IA, seem, seim, siem, siim, seev, seiv, siev, siiv = ti.timeIntegration(self.params)
+            (
+                rates_exc,
+                rates_inh,
+                t,
+                mufe,
+                mufi,
+                IA,
+                seem,
+                seim,
+                siem,
+                siim,
+                seev,
+                seiv,
+                siev,
+                siiv,
+            ) = ti.timeIntegration(self.params)
 
         # convert output from kHz to Hz
         rates_exc = rates_exc * 1000.0

diff --git a/neurolib/optimize/evolution/evolution.py b/neurolib/optimize/evolution/evolution.py
@@ -19,11 +19,23 @@
 
 
 class Evolution:
+    """Evolutionary parameter optimization.
+    """
+
     def __init__(
-        self, evalFunction, parameterSpace, weightList, model=None, hdf_filename="evolution.hdf", ncores=None, POP_INIT_SIZE=100, POP_SIZE=20, NGEN=10, CXPB=1 - 0.96,
+        self,
+        evalFunction,
+        parameterSpace,
+        weightList,
+        model=None,
+        hdf_filename="evolution.hdf",
+        ncores=None,
+        POP_INIT_SIZE=100,
+        POP_SIZE=20,
+        NGEN=10,
+        CXPB=1 - 0.96,
     ):
-        """Evolutionary parameter optimization
-        
+        """
         :param model: Model to run
         :type model: Model
         :param evalFunction: Evaluation function of a run that provides a fitness vector and simulation outputs
@@ -49,7 +61,14 @@ def __init__(
 
         # initialize pypet environment
         # env = pp.Environment(trajectory=trajectoryName, filename=trajectoryFileName)
-        env = pp.Environment(trajectory=trajectoryName, filename=trajectoryFileName, use_pool=False, multiproc=True, ncores=ncores, complevel=9,)
+        env = pp.Environment(
+            trajectory=trajectoryName,
+            filename=trajectoryFileName,
+            use_pool=False,
+            multiproc=True,
+            ncores=ncores,
+            complevel=9,
+        )
 
         # Get the trajectory from the environment
         traj = env.traj
@@ -74,10 +93,13 @@ def __init__(
 
         self.initialPopulationSimulated = False
 
-        # ------------- settings
+        # settings
         self.verbose = False
 
-        # ------------- define parameters
+        # environment parameters
+        self.evaluationCounter = 0
+
+        # simulation parameters
         self.ParametersInterval = parameterSpace.named_tuple_constructor
         self.paramInterval = parameterSpace.named_tuple
 
@@ -92,7 +114,7 @@ def __init__(
             self.traj, self.paramInterval, self.ParametersInterval, self.POP_SIZE, self.CXPB, self.NGEN, self.model,
         )
 
-        # ------------- initialize population
+        # initialize population
         self.last_id = 0
         self.pop = self.toolbox.population(n=self.POP_INIT_SIZE)
         self.pop = self.tagPopulation(self.pop)
@@ -164,7 +186,10 @@ def initDEAP(self, toolbox, pypetEnvironment, paramInterval, evalFunction, weigh
         # need to create a lambda funciton because du.generateRandomParams wants an argument but
         # toolbox.register cannot pass an argument to it.
         toolbox.register(
-            "individual", deap.tools.initIterate, deap.creator.Individual, lambda: du.generate_random_pars_adapt(paramInterval),
+            "individual",
+            deap.tools.initIterate,
+            deap.creator.Individual,
+            lambda: du.generate_random_pars_adapt(paramInterval),
         )
         toolbox.register("population", deap.tools.initRepeat, list, toolbox.individual)
 
@@ -196,21 +221,28 @@ def evalPopulationUsingPypet(self, traj, toolbox, pop, gIdx):
         # This is for only having the cartesian product
         # between ``generation x (ind_idx AND individual)``, so that every individual has just one
         # unique index within a generation.
-        traj.f_expand(pp.cartesian_product({"generation": [gIdx], "id": [x.id for x in pop], "individual": [list(x) for x in pop],}, [("id", "individual"), "generation"],))  # the current generation  # unique id of each individual
-        # SIMULUATE INDIVIDUALS
+        traj.f_expand(
+            pp.cartesian_product(
+                {"generation": [gIdx], "id": [x.id for x in pop], "individual": [list(x) for x in pop],},
+                [("id", "individual"), "generation"],
+            )
+        )  # the current generation  # unique id of each individual
 
-        results = toolbox.map(toolbox.evaluate)
+        # increment the evaluationCounter
+        self.evaluationCounter += len(pop)
 
-        assert len(results) > 0, "No results returned from simulations."
+        # run simulations for one generation
+        evolutionResult = toolbox.map(toolbox.evaluate)
 
-        for idx, result in enumerate(results):
-            runIdx, packed_result = result
-            # this is the return from the evaluation function
-            fitnesses_result, outputs = packed_result
+        assert len(evolutionResult) > 0, "No results returned from simulations."
 
+        for idx, result in enumerate(evolutionResult):
+            runIndex, packedReturnFromEvalFunction = result
+            # this is the return from the evaluation function
+            fitnessesResult, outputs = packedReturnFromEvalFunction
             # store outputs of simulations in population
             pop[idx].outputs = outputs
-            pop[idx].fitness.values = fitnesses_result
+            pop[idx].fitness.values = fitnessesResult
             # mean fitness value
             pop[idx].fitness.score = np.nansum(pop[idx].fitness.wvalues) / (len(pop[idx].fitness.wvalues))
         return pop
@@ -308,7 +340,9 @@ def runEvolution(self):
 
             if self.verbose:
                 eu.printParamDist(self.pop, self.paramInterval, self.gIdx)
-                eu.printPopFitnessStats(self.pop, self.paramInterval, self.gIdx, draw_scattermatrix=True, save_plots="evo")
+                eu.printPopFitnessStats(
+                    self.pop, self.paramInterval, self.gIdx, draw_scattermatrix=True, save_plots="evo"
+                )
 
             # save all simulation data to pypet
             try:
@@ -345,20 +379,20 @@ def loadResults(self, filename=None, trajectoryName=None):
         """
         if filename == None:
             filename = self.HDF_FILE
-        trajLoaded = pu.loadPypetTrajectory(filename, trajectoryName)
-        return trajLoaded
+        self.traj = pu.loadPypetTrajectory(filename, trajectoryName)
 
     def getScoresDuringEvolution(self, traj=None, drop_first=True, reverse=False):
         if traj == None:
             traj = self.traj
 
         generation_names = list(traj.results.evolution.f_to_dict(nested=True).keys())
 
-        if drop_first:
-            # drop first (initial) generation 0
-            generation_names = generation_names[1:]
         if reverse:
             generation_names = generation_names[::-1]
+        if drop_first:
+            # drop first (initial) generation 0
+            # generation_names = generation_names[1:]
+            generation_names.remove("gen_000000")
 
         npop = len(traj.results.evolution[generation_names[0]].scores)
 

diff --git a/resources/evolution_minimal.png b/resources/evolution_minimal.png
diff --git a/resources/exploration_aln.png b/resources/exploration_aln.png