generalized_opt_routine

.github/workflows/unit_tests.yml

@@ -43,7 +43,7 @@ jobs:
       - name: Install HNN-core
         shell: bash -el {0}
         run: |
-          pip install --verbose '.[gui]'
+          pip install --verbose '.[gui,opt]'
       - name: Lint with flake8
         shell: bash -el {0}
         run: |
@@ -55,4 +55,4 @@ jobs:
       - name: Upload code coverage
         shell: bash -el {0}
         run: |
-          bash <(curl -s https://codecov.io/bash) -f ./coverage.xml
+          bash <(curl -s https://codecov.io/bash) -f ./coverage.xml

.github/workflows/windows_unit_tests.yml

@@ -40,7 +40,7 @@ jobs:
           python -m pip install --upgrade pip
           pip install flake8 pytest pytest-cov
           pip install psutil joblib
-          pip install --verbose .[gui]
+          pip install --verbose .[gui,opt]
       - name: Test with pytest
         shell: cmd
         run: |

README.rst

@@ -28,6 +28,11 @@ GUI
 * ipywidgets (<=7.7.1)
 * voila (<=0.3.6)
 
+Optimization
+~~~~~~~~~~~~
+
+* scikit-learn
+
 Parallel processing
 ~~~~~~~~~~~~~~~~~~~
 
@@ -61,6 +66,13 @@ To check if everything worked fine, you can do::
 
 and it should not give any error messages.
 
+**Installing optimization dependencies**
+
+If you are using bayesian optimization, then scikit-learn is required. Install 
+hnn-core with scikit-learn using the following command::
+
+    $ pip install hnn_core[opt]
+
 **GUI installation**
 
 To install the GUI dependencies along with ``hnn-core``, a simple tweak to the above command is needed::

doc/api.rst

@@ -40,7 +40,7 @@ Optimization (:py:mod:`hnn_core.optimization`):
 .. autosummary::
    :toctree: generated/
 
-   optimize_evoked
+   Optimizer
 
 Dipole (:py:mod:`hnn_core.dipole`):
 -----------------------------------

doc/whats_new.rst

@@ -38,6 +38,10 @@ API
 - :func:`~hnn_core.Dipole.write` and :func:`~hnn_core.Dipole.read_dipoles`
   now support hdf5 format for read/write Dipole object, by
   `Rajat Partani`_ in :gh:`648`
+
+- Add ability to optimize parameters associated with evoked drives and plot 
+  convergence. User can constrain parameter ranges and specify solver,
+  by `Carolina Fernandez Pujol`_ in :gh:`652`
 
 .. _0.3:
 

examples/howto/optimize_evoked.py

@@ -7,21 +7,21 @@
 of the model simulation to match an experimental dipole waveform.
 """
 
-# Authors: Blake Caldwell <[email protected]>
+# Authors: Carolina Fernandez <[email protected]>
+#          Nick Tolley <[email protected]>
+#          Ryan Thorpe <[email protected]>
 #          Mainak Jas <[email protected]>
 
 import os.path as op
 
-import numpy as np
 import matplotlib.pyplot as plt
 
 ###############################################################################
 # Let us import hnn_core
 
 import hnn_core
-from hnn_core import (MPIBackend, jones_2009_model, read_params,
-                      simulate_dipole, read_dipole)
-
+from hnn_core import (MPIBackend, jones_2009_model, simulate_dipole,
+                      read_dipole)
 
 hnn_core_root = op.join(op.dirname(hnn_core.__file__))
 
@@ -34,6 +34,7 @@
 # This is a different experiment than the one to which the base parameters were
 # tuned. So, the initial RMSE will be large, giving the optimization procedure
 # a lot to work with.
+
 from urllib.request import urlretrieve
 
 data_url = ('https://raw.githubusercontent.com/jonescompneurolab/hnn/master/'
@@ -42,82 +43,172 @@
 exp_dpl = read_dipole('S1_SupraT.txt')
 
 ###############################################################################
-# Read the base parameters from a file
-params_fname = op.join(hnn_core_root, 'param', 'default.json')
-params = read_params(params_fname)
+# Let’s then simulate the dipole with some initial parameters.
 
-###############################################################################
-# Let's first simulate the dipole with some initial parameters. The parameter
-# definitions also contain the drives. Even though we could add drives
-# explicitly through our API
-# (see :ref:`sphx_glr_auto_examples_workflows_plot_simulate_evoked.py`),
-# for conciseness,
-# we add them automatically from the parameter files
-
-scale_factor = 3000.
-smooth_window_len = 30.
 tstop = exp_dpl.times[-1]
-net = jones_2009_model(params=params, add_drives_from_params=True)
-with MPIBackend(n_procs=n_procs):
-    print("Running simulation with initial parameters")
-    initial_dpl = simulate_dipole(net, tstop=tstop, n_trials=1)[0]
-    initial_dpl = initial_dpl.scale(scale_factor).smooth(smooth_window_len)
+scale_factor = 3000
+smooth_window_len = 30
+
+net_init = jones_2009_model()
+
+# Proximal 1
+weights_ampa_p1 = {'L2_basket': 0.2913, 'L2_pyramidal': 0.9337,
+                   'L5_basket': 0.1951, 'L5_pyramidal': 0.3602}
+weights_nmda_p1 = {'L2_basket': 0.9240, 'L2_pyramidal': 0.0845,
+                   'L5_basket': 0.5849, 'L5_pyramidal': 0.65105}
+synaptic_delays_p = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
+                     'L5_basket': 1., 'L5_pyramidal': 1.}
+net_init.add_evoked_drive('evprox1',
+                          mu=5.6813,
+                          sigma=20.3969,
+                          numspikes=1,
+                          location='proximal',
+                          weights_ampa=weights_ampa_p1,
+                          weights_nmda=weights_nmda_p1,
+                          synaptic_delays=synaptic_delays_p)
+
+# Distal
+weights_ampa_d1 = {'L2_basket': 0.8037, 'L2_pyramidal': 0.5738,
+                   'L5_pyramidal': 0.3626}
+weights_nmda_d1 = {'L2_basket': 0.2492, 'L2_pyramidal': 0.6183,
+                   'L5_pyramidal': 0.1121}
+synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
+                      'L5_pyramidal': 0.1}
+net_init.add_evoked_drive('evdist1',
+                          mu=58.6539,
+                          sigma=5.5810,
+                          numspikes=1,
+                          location='distal',
+                          weights_ampa=weights_ampa_d1,
+                          weights_nmda=weights_nmda_d1,
+                          synaptic_delays=synaptic_delays_d1)
+
+# Proximal 2
+weights_ampa_p2 = {'L2_basket': 0.01, 'L2_pyramidal': 0.01, 'L5_basket': 0.01,
+                   'L5_pyramidal': 0.01}
+weights_nmda_p2 = {'L2_basket': 0.01, 'L2_pyramidal': 0.01, 'L5_basket': 0.01,
+                   'L5_pyramidal': 0.01}
+net_init.add_evoked_drive('evprox2',
+                          mu=80,
+                          sigma=1,
+                          numspikes=1,
+                          location='proximal',
+                          weights_ampa=weights_ampa_p2,
+                          weights_nmda=weights_nmda_p2,
+                          synaptic_delays=synaptic_delays_p)
+
+with MPIBackend(n_procs=n_procs, mpi_cmd='mpiexec'):
+    init_dpl = simulate_dipole(net_init, tstop=tstop, n_trials=1)[0]
+init_dpl.scale(scale_factor)
+init_dpl.smooth(smooth_window_len)
 
 ###############################################################################
 # Now we start the optimization!
-
-from hnn_core.optimization import optimize_evoked
-
-with MPIBackend(n_procs=n_procs):
-    net_opt = optimize_evoked(net, tstop=tstop, n_trials=1,
-                              target_dpl=exp_dpl, initial_dpl=initial_dpl,
-                              scale_factor=scale_factor,
-                              smooth_window_len=smooth_window_len)
-
-###############################################################################
-# Now, let's simulate the dipole with the optimized drive parameters.
-with MPIBackend(n_procs=n_procs):
-    best_dpl = simulate_dipole(net_opt, tstop=tstop, n_trials=1)[0]
-    best_dpl = best_dpl.scale(scale_factor).smooth(smooth_window_len)
+#
+# First, we define a function that will tell the optimization routine how to
+# modify the network drive parameters. The function will take in the Network
+# object with no attached drives, and a dictionary of the paramters we wish to
+# optimize.
+
+
+def set_params(net, params):
+
+    # Proximal 1
+    net.add_evoked_drive('evprox1',
+                         mu=5.6813,
+                         sigma=20.3969,
+                         numspikes=1,
+                         location='proximal',
+                         weights_ampa=weights_ampa_p1,
+                         weights_nmda=weights_nmda_p1,
+                         synaptic_delays=synaptic_delays_p)
+
+    # Distal
+    net.add_evoked_drive('evdist1',
+                         mu=58.6539,
+                         sigma=5.5810,
+                         numspikes=1,
+                         location='distal',
+                         weights_ampa=weights_ampa_d1,
+                         weights_nmda=weights_nmda_d1,
+                         synaptic_delays=synaptic_delays_d1)
+
+    # Proximal 2
+    weights_ampa_p2 = {'L2_basket':
+                       params['evprox2_ampa_L2_basket'],
+                       'L2_pyramidal':
+                       params['evprox2_ampa_L2_pyramidal'],
+                       'L5_basket':
+                       params['evprox2_ampa_L5_basket'],
+                       'L5_pyramidal':
+                       params['evprox2_ampa_L5_pyramidal']}
+    weights_nmda_p2 = {'L2_basket':
+                       params['evprox2_nmda_L2_basket'],
+                       'L2_pyramidal':
+                       params['evprox2_nmda_L2_pyramidal'],
+                       'L5_basket':
+                       params['evprox2_nmda_L5_basket'],
+                       'L5_pyramidal':
+                       params['evprox2_nmda_L5_pyramidal']}
+    net.add_evoked_drive('evprox2',
+                         mu=params['evprox2_mu'],
+                         sigma=params['evprox2_sigma'],
+                         numspikes=1,
+                         location='proximal',
+                         weights_ampa=weights_ampa_p2,
+                         weights_nmda=weights_nmda_p2,
+                         synaptic_delays=synaptic_delays_p)
 
 ###############################################################################
-# Then, we can plot the pre- and post-optimization simulations alongside the
-# experimental data. 
+# Then, we define the constraints.
+#
+# The constraints must be a dictionary of tuples where the first value in each
+# tuple is the lower bound and the second value is the upper bound for the
+# corresponding parameter.
+#
+# The following synaptic weight parameter ranges (units of micro-siemens)
+# were chosen so as to keep the model in physiologically realistic regimes.
 
-fig, axes = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
 
-exp_dpl.plot(ax=axes[0], layer='agg', show=False, color='tab:blue')
-initial_dpl.plot(ax=axes[0], layer='agg', show=False, color='tab:orange')
-best_dpl.plot(ax=axes[0], layer='agg', show=False, color='tab:green')
-axes[0].legend(['experimental', 'initial', 'optimized'])
-net_opt.cell_response.plot_spikes_hist(ax=axes[1])
+constraints = dict({'evprox2_ampa_L2_basket': (0.01, 1.),
+                    'evprox2_ampa_L2_pyramidal': (0.01, 1.),
+                    'evprox2_ampa_L5_basket': (0.01, 1.),
+                    'evprox2_ampa_L5_pyramidal': (0.01, 1.),
+                    'evprox2_nmda_L2_basket': (0.01, 1.),
+                    'evprox2_nmda_L2_pyramidal': (0.01, 1.),
+                    'evprox2_nmda_L5_basket': (0.01, 1.),
+                    'evprox2_nmda_L5_pyramidal': (0.01, 1.),
+                    'evprox2_mu': (100., 120.),
+                    'evprox2_sigma': (2., 30.)})
 
 ###############################################################################
-# Finally, let's explore which parameters were changed to cause the 
-# improved dipole fit. 
-# As an example, we will look at the new dynamics of the first proximal drive, 
-# as well as the synaptic weight of its layer 2/3 pyramidal AMPA receptors.
+# Now we define and fit the optimizer.
 
-from hnn_core.network import pick_connection
+from hnn_core.optimization import Optimizer
 
-dynamics_opt = net_opt.external_drives['evdist1']['dynamics']
+net = jones_2009_model()
+optim = Optimizer(net, tstop=tstop, constraints=constraints,
+                  set_params=set_params, scale_factor=scale_factor,
+                  smooth_window_len=smooth_window_len, max_iter=40)
+with MPIBackend(n_procs=n_procs, mpi_cmd='mpiexec'):
+    optim.fit(exp_dpl)
 
-conn_indices = pick_connection(net=net_opt, src_gids='evprox1', target_gids='L2_pyramidal',
-    loc='proximal', receptor='ampa')
-conn_idx = conn_indices[0]
-weight_opt = net_opt.connectivity[conn_idx]
+###############################################################################
+# Finally, we can plot the experimental data alongside the post-optimization
+# simulation dipole as well as the convergence plot.
 
-print("Optimized dynamic properties: ", dynamics_opt)
-print("\nOptimized weight: ", weight_opt)
+with MPIBackend(n_procs=n_procs, mpi_cmd='mpiexec'):
+    opt_dpl = simulate_dipole(optim.net_, tstop=tstop, n_trials=1)[0]
+opt_dpl.scale(scale_factor)
+opt_dpl.smooth(smooth_window_len)
 
-###############################################################################
-# Let's compare to the initial dynamic properties and weight.
-dynamics_initial = net.external_drives['evdist1']['dynamics']
+fig, axes = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
 
-conn_indices = pick_connection(net=net, src_gids='evprox1', target_gids='L2_pyramidal',
-    loc='proximal', receptor='ampa')
-conn_idx = conn_indices[0]
-weight_initial = net.connectivity[conn_idx]
+# plot original
+exp_dpl.plot(ax=axes[0], layer='agg', show=False, color='tab:blue')
+init_dpl.plot(ax=axes[0], layer='agg', show=False, color='tab:orange')
+opt_dpl.plot(ax=axes[0], layer='agg', show=False, color='tab:green')
+axes[0].legend(['experimental', 'initial', 'optimized'])
+optim.net_.cell_response.plot_spikes_hist(ax=axes[1])
 
-print("Initial dynamic properties: ", dynamics_initial)
-print("\nInitial weight: ", weight_initial)
+fig1 = optim.plot_convergence()