{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hyperparameter optimization\n", "\n", "[![Download JupyterNotebook](https://img.shields.io/badge/Download-Notebook-orange?style=for-the-badge&logo=Jupyter)](https://raw.githubusercontent.com/ANNarchy/ANNarchy.github.io/master/notebooks/BayesianOptimization.ipynb) [![Download JupyterNotebook](https://img.shields.io/badge/Open_in-Colab-blue?style=for-the-badge&logo=Jupyter)](https://colab.research.google.com/github/ANNarchy/ANNarchy.github.io/blob/master/notebooks/BayesianOptimization.ipynb)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:15:58.101138Z", "iopub.status.busy": "2026-01-05T13:15:58.100803Z", "iopub.status.idle": "2026-01-05T13:15:58.104746Z", "shell.execute_reply": "2026-01-05T13:15:58.104094Z" } }, "outputs": [], "source": [ "#!pip install ANNarchy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Most of the work in computational neuroscience is to guess the values of parameters which are not constrained by the biology. The most basic approach is to simply try out different values, run the simulation, reason about why the results are not what you want, change some parameters, run again, etc. It is very easy to get lost in this process and it requires a great deal of intuition about how the model works.\n", "\n", "If you are able to define an objective function for your model (a single number that tells how well your model performs), you can use search algorithms to find this hyperparameters automatically, at the cost of running your model multiple times.\n", "\n", "Let's take the example of a rate-coded model depending on two hyperparameters `a` and `b`, where is the objective is to have a minimal activity after 1 s of simulation (dummy example):\n", "\n", "```python\n", "net = MyNetwork()\n", "net.compile()\n", "\n", "def run(a, b):\n", " net.pop.a = a\n", " net.pop.b = b\n", " \n", " net.simulate(1000.)\n", " \n", " return (net.pop.r)**2\n", "```\n", "\n", "**Grid search** would iterate over all possible values of the parameters to perform the search:\n", "\n", "```python\n", "min_loss = 1000.\n", "for a in np.linspace(0.0, 1.0, 100):\n", " for b in np.linspace(0.0, 1.0, 100):\n", " loss = run(a, b)\n", " if loss < min_loss:\n", " min_loss = loss\n", " a_best = a ; b_best = b\n", "```\n", "\n", "If you try 100 values for each parameters, you need 10000 simulations to find your parameters. The number of simulations explodes with the number of free parameters. Moreover, you cannot stop the search before the end, as you could miss the interesting region.\n", "\n", "**Random search** samples blindly values for the hyperparameters:\n", "\n", "```python\n", "min_loss = 1000.\n", "for _ in range(1000):\n", " a = np.random.uniform(0.0, 1.0)\n", " b = np.random.uniform(0.0, 1.0)\n", " loss = run(a, b)\n", " if loss < min_loss:\n", " min_loss = loss\n", " a_best = a ; b_best = b\n", "```\n", "\n", "If you are lucky, you may find a good solution quite early in the search, so you can stop it when the loss is below a desired threshold. The main drawback is that the search may spend a lot of time in uninteresting regions: it does not learn anything between two samples.\n", "\n", "An often much more efficient search method is **Bayesian optimization** (also called sequential model-based optimization - SMBO). It is a form of random search that updates beliefs on the hyperparameters. In short, if some parameter values do not lead to good values of the objective function in early samples, they will not be used in later samples. The search becomes more and more focused on the interesting regions of the hyperparameter space. \n", "\n", "As always with Python, there are many libraries for that, including:\n", "\n", "* `hyperopt` \n", "* `optuna` \n", "* `talos` (for keras models) \n", "\n", "This notebook demonstrates how to use `hyperopt` to find some hyperparameters of the COBA models already included in the ANNarchy examples:\n", "\n", "\n", "\n", "Additionally, we will use the `tensorboard` extension to visualize the dependency between the parameters and the objective function." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:15:58.107036Z", "iopub.status.busy": "2026-01-05T13:15:58.106811Z", "iopub.status.idle": "2026-01-05T13:15:58.766968Z", "shell.execute_reply": "2026-01-05T13:15:58.766262Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ANNarchy 5.0 (5.0.0) on linux (posix).\n" ] } ], "source": [ "import numpy as np\n", "\n", "import ANNarchy as ann\n", "from ANNarchy.extensions.tensorboard import Logger" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:15:58.799591Z", "iopub.status.busy": "2026-01-05T13:15:58.799092Z", "iopub.status.idle": "2026-01-05T13:16:08.761790Z", "shell.execute_reply": "2026-01-05T13:16:08.761130Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Compiling network 1... " ] }, { "name": "stdout", "output_type": "stream", "text": [ "OK \n" ] } ], "source": [ "COBA = ann.Neuron(\n", " parameters = dict(\n", " El = -60.0 ,\n", " Vr = -60.0 ,\n", " Erev_exc = 0.0 ,\n", " Erev_inh = -80.0 ,\n", " Vt = -50.0 ,\n", " tau = 20.0 ,\n", " tau_exc = 5.0 ,\n", " tau_inh = 10.0 ,\n", " I = 20.0 ,\n", " ),\n", " equations = [\n", " 'tau * dv/dt = (El - v) + g_exc * (Erev_exc - v) + g_inh * (Erev_inh - v ) + I',\n", " 'tau_exc * dg_exc/dt = - g_exc',\n", " 'tau_inh * dg_inh/dt = - g_inh',\n", " ],\n", " spike = \"v > Vt\",\n", " reset = \"v = Vr\",\n", " refractory = 5.0\n", ")\n", "\n", "class COBANetwork (ann.Network):\n", " \n", " def __init__(self, w_exc:float = 0.6, w_inh:float = 6.7):\n", "\n", " self.P = self.create(geometry=4000, neuron=COBA)\n", " self.P.v = ann.Normal(-55.0, 5.0)\n", " self.P.g_exc = ann.Normal(4.0, 1.5)\n", " self.P.g_inh = ann.Normal(20.0, 12.0)\n", "\n", " self.Ce = self.connect(pre=self.P[:3200], post=self.P, target='exc')\n", " self.Ce.fixed_probability(weights=w_exc, probability=0.02)\n", "\n", " self.Ci = self.connect(pre=self.P[3200:], post=self.P, target='inh')\n", " self.Ci.fixed_probability(weights=w_inh, probability=0.02)\n", "\n", " self.m = self.monitor(self.P, ['spike'])\n", "\n", "net = COBANetwork(w_exc=0.6, w_inh=6.7, dt=0.1)\n", "net.compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the default parameters, the COBA network fires at around 20 Hz:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:16:08.764108Z", "iopub.status.busy": "2026-01-05T13:16:08.763856Z", "iopub.status.idle": "2026-01-05T13:16:08.949987Z", "shell.execute_reply": "2026-01-05T13:16:08.949317Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Firing rate: 21.899\n" ] } ], "source": [ "net.simulate(1000.0)\n", "\n", "data = net.m.get('spike')\n", "fr = net.m.mean_fr(data)\n", "print(\"Firing rate:\", fr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's suppose we now want the network to fire at 30 Hz. Which parameters should we change to obtain that value?\n", "\n", "Many parameters might influence the firing rate of the network (if not all). Here, we make the assumption that the weight values for the excitatory connections (0.6) and inhibitory ones (6.7) are the most critical ones.\n", "\n", "Let's start by importing `hyperopt` (after installing it with `pip install hyperopt`):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:16:08.951883Z", "iopub.status.busy": "2026-01-05T13:16:08.951657Z", "iopub.status.idle": "2026-01-05T13:16:09.595216Z", "shell.execute_reply": "2026-01-05T13:16:09.594509Z" } }, "outputs": [], "source": [ "from hyperopt import fmin, tpe, hp, STATUS_OK" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define a `trial()` method taking values for the two hyperparameters as inputs. It starts by resetting the network, sets the excitatory and inhibitory weights to the desired value, simulates for one second, computes the mean firing rate of the population, logs the parameters and finally returns the objective function: the squared error between the recorded firing rate and 30 Hz." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:16:09.598009Z", "iopub.status.busy": "2026-01-05T13:16:09.597563Z", "iopub.status.idle": "2026-01-05T13:16:09.607350Z", "shell.execute_reply": "2026-01-05T13:16:09.606648Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Logging in runs/Jan05_14-16-09_stand-pc\n" ] } ], "source": [ "logger = Logger()\n", "\n", "def trial(args):\n", " \n", " # Retrieve the parameters\n", " w_exc = args[0]\n", " w_inh = args[1]\n", " \n", " # Create the network as a copy of net with different parameters (avoids recompilation)\n", " new_net = net.copy(w_exc, w_inh)\n", " \n", " # Simulate 1 second\n", " new_net.simulate(1000.0)\n", "\n", " # Retrieve the spike recordings\n", " spikes = new_net.m.get('spike')\n", "\n", " # Compute the population firing rate\n", " fr = new_net.m.mean_fr(spikes)\n", " \n", " # Compute a quadratic loss around 30 Hz\n", " loss = 0.001 * (fr - 30.0)**2 \n", " \n", " # Log the parameters\n", " logger.add_parameters({'w_exc': w_exc, 'w_inh': w_inh},\n", " {'loss': loss, 'firing_rate': fr})\n", " \n", " # Delete the network to avoid memory leak\n", " new_net.clear()\n", " \n", " return {\n", " 'loss': loss,\n", " 'status': STATUS_OK,\n", " # -- store other results like this\n", " 'fr': fr,\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check that the default parameters indeed lead to a firing rate of 20 Hz:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:16:09.610068Z", "iopub.status.busy": "2026-01-05T13:16:09.609779Z", "iopub.status.idle": "2026-01-05T13:16:10.227193Z", "shell.execute_reply": "2026-01-05T13:16:10.226451Z" } }, "outputs": [ { "data": { "text/plain": [ "{'loss': 0.05009454506250008, 'status': 'ok', 'fr': 22.922249999999995}" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trial([0.6, 6.7])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use `hyperopt` to find the hyperparameters making the network fire at 30 Hz.\n", "\n", "The `fmin()` function takes:\n", "\n", "* `fn`: the objective function for a set of parameters.\n", "* `space`: the search space for the hyperparameters (the prior). \n", "* `algo`: which algorithm to use, either tpe.suggest or random.suggest\n", "* `max_evals`: number of samples (simulations) to make.\n", "\n", "Here, we will sample the excitatory weights between 0.1 and 1, the inhibitory ones between 1 and 10. Of course, the smaller the range, the better. Refer to the doc of hyperopt for other sampling priors." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2026-01-05T13:16:10.229927Z", "iopub.status.busy": "2026-01-05T13:16:10.229613Z", "iopub.status.idle": "2026-01-05T13:17:17.679549Z", "shell.execute_reply": "2026-01-05T13:17:17.678896Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\r", " 0%| | 0/100 [00:00