Example 1: A Simple Neuroptimiser with a problem from IOH

This example demonstrates how to use the Neuroptimiser library to solve an optimisation problem from the IOH framework. The problem is defined by its ID and instance, and the Neuroptimiser is configured with a set of parameters for the agents.

1. Setup

Import necessary libraries and set up the environment for plotting.

# Import necessary libraries
import os
import shutil
from copy import deepcopy
import random

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import ListedColormap

from neuroptimiser import NeurOptimiser
from ioh import get_problem

# Set up the plotting style and parameters
width = 6.5
height = width * 0.618
res_dpi = 333

seed = 42
random.seed(seed)
np.random.seed(seed)

root_figures = "./figures/"

plt.style.context('paper')
sns.set_style("ticks")
plt.rcParams['font.family']         = 'serif'
plt.rcParams['axes.labelsize']      = 11
plt.rcParams['xtick.labelsize']     = 11
plt.rcParams['ytick.labelsize']     = 11
plt.rcParams['legend.fontsize']     = 11
plt.rcParams['axes.grid']           = False
plt.rcParams['axes.spines.top']     = True
plt.rcParams['axes.spines.right']   = True
plt.rcParams['axes.edgecolor']      = 'black'
plt.rcParams['xtick.color']         = 'black'
plt.rcParams['ytick.color']         = 'black'
plt.rcParams['axes.labelcolor']     = 'black'

if shutil.which("pdflatex") is not None:
    plt.rcParams['text.usetex']     = True
    plt.rcParams['font.serif']      = ['Computer Modern Roman']
else:
    plt.rcParams['text.usetex']     = False

def save_this(fig, _in=""):
    if is_saving:
        os.makedirs(f"{root_figures}/" + _in, exist_ok=True)
        fig.savefig(root_figures + _in + "/" + prefix_filename + ".png", transparent=True, dpi=res_dpi, bbox_inches="tight")
    else:
        plt.show()

2. Parameter definition

In this section, we define the problem parameters, including the number of steps, problem ID, instance, dimensions. In addition, we set the custom parameters for the Neuroptimiser, which will replace the default parameters of the spiking core.

# Define the problem parameters
num_steps       = 1000       # Number of iterations
problem_id      = 3         # Problem ID from the IOH framework
problem_ins     = 1         # Problem instance
num_dimensions  = 2         # Number of dimensions for the problem
num_agents      = 30        # Number of agents/units in the Neuroptimiser
num_neighbours  = 15        # Number of neighbours for each agent/unit

# Define the Neuroptimiser parameters
neuropt_name    = "neuropt_Ex1-LinRand"

config_params = dict(
    num_iterations  = num_steps,
    num_agents      = num_agents,
    spiking_core    = "TwoDimSpikingCore",
    num_neighbours  = num_neighbours,
    neuron_topology = "2dr",
    unit_topology   = "random",
)

core_params = [   # Specify the custom parameters for the spiking core
    {"name": "izhikevich",
     "coeffs": "random",
     "spk_cond": "l1",
     "hs_operator": "differential",
     "hs_variant": "current-to-rand",
     "sel_mode": "greedy",
     },
]

# Additional parameters
is_saving       = True      # Whether to save the figures or not
prefix_filename = (f"{neuropt_name}_"
                   f"{problem_id}p_"
                   f"{problem_ins}i_"
                   f"{num_dimensions}d_"
                   f"{num_steps}s_"
                   f"{num_agents}u")
plt.close()
if is_saving:
    print(f"prefix_filename: '{prefix_filename}'")

prefix_filename: 'neuropt_Ex1-LinRand_3p_1i_2d_1000s_30u'

3. Problem setup and optimisation

We first retrieve the problem from the IOH framework.

# Get the problem from the IOH framework
problem = get_problem(fid=problem_id,
                      instance=problem_ins,
                      dimension=num_dimensions,
                      )
problem.reset()
print(problem)

<RealSingleObjectiveProblem 3. Rastrigin (iid=1 dim=2)>

Then, we instantiate the Neuroptimiser with the defined parameters and solve the problem. The debug_mode is set to True to enable detailed logging of the optimisation process.

# Instantiation
optimiser = NeurOptimiser(config_params, core_params)

# Solve the problem
optimiser.solve(problem, debug_mode=True)

[neuropt:log] Debug mode is enabled. Monitoring will be activated.
[neuropt:log] Parameters are set up.
[neuropt:log] Initial positions and topologies are set up.
[neuropt:log] Tensor contraction layer, neighbourhood manager, and high-level selection unit are created.
[neuropt:log] Population of nheuristic units is created.
[neuropt:log] Connections between nheuristic units and auxiliary processes are established.
[neuropt:log] Monitors are set up.
[neuropt:log] Starting simulation with 1000 iterations...
... step: 0, best fitness: -457.0638427734375
... step: 100, best fitness: -461.8068542480469
... step: 200, best fitness: -461.9642639160156
... step: 300, best fitness: -462.08636474609375
... step: 400, best fitness: -462.0899963378906
... step: 500, best fitness: -462.0899963378906
... step: 600, best fitness: -462.0899963378906
... step: 700, best fitness: -462.0899963378906
... step: 800, best fitness: -462.0899963378906
... step: 900, best fitness: -462.0899963378906
... step: 999, best fitness: -462.0899963378906
[neuropt:log] Simulation completed. Fetching monitor data... done

(array([-2.34097064,  2.29994103]), -462.0899963378906)

# Show the overall configuration parameters of the optimiser
print(optimiser.config_params)

{'num_iterations': 1000, 'num_agents': 30, 'spiking_core': 'TwoDimSpikingCore', 'num_neighbours': 15, 'neuron_topology': '2dr', 'unit_topology': 'random', 'search_space': array([[-5.,  5.],
       [-5.,  5.]]), 'seed': 69, 'num_dimensions': 2, 'function': <function AbstractSolver._rescale_problem.<locals>.scaled_problem at 0x120b35c60>, 'core_params': {}}

# Show the core parameters used in the optimisation
pd.DataFrame(optimiser.core_params)

	name	coeffs	spk_cond	hs_operator	hs_variant	sel_mode	spk_alpha	seed	noise_std	hs_params	...	thr_k	spk_weights	approx	thr_alpha	thr_mode	thr_max	alpha	thr_min	spk_q_ord	init_position
0	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	69	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.250919762305275, 0.9014286128198323]
1	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	70	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.4639878836228102, 0.1973169683940732]
2	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	71	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.687962719115127, -0.6880109593275947]
3	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	72	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.8838327756636011, 0.7323522915498704]
4	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	73	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.2022300234864176, 0.416145155592091]
5	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	74	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.9588310114083951, 0.9398197043239886]
6	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	75	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.6648852816008435, -0.5753217786434477]
7	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	76	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.6363500655857988, -0.6331909802931324]
8	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	77	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.39151551408092455, 0.04951286326447568]
9	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	78	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.13610996271576847, -0.4175417196039162]
10	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	79	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.22370578944475894, -0.7210122786959163]
11	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	80	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.4157107029295637, -0.2672763134126166]
12	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	81	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.08786003156592814, 0.5703519227860272]
13	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	82	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.6006524356832805, 0.02846887682722321]
14	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	83	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.18482913772408494, -0.9070991745600046]
15	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	84	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.21508970380287673, -0.6589517526254169]
16	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	85	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.869896814029441, 0.8977710745066665]
17	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	86	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.9312640661491187, 0.6167946962329223]
18	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	87	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.39077246165325863, -0.8046557719872323]
19	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	88	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.3684660530243138, -0.1196950125207974]
20	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	89	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.7559235303104423, -0.00964617977745963]
21	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	90	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.9312229577695632, 0.8186408041575641]
22	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	91	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.48244003679996617, 0.32504456870796394]
23	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	92	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.3765778478211781, 0.040136042355621626]
24	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	93	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.0934205586865593, -0.6302910889489459]
25	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	94	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.9391692555291171, 0.5502656467222291]
26	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	95	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.8789978831283782, 0.7896547008552977]
27	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	96	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[0.19579995762217028, 0.8437484700462337]
28	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	97	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.823014995896161, -0.6080342751617096]
29	izhikevich	random	l1	differential	current-to-rand	greedy	0.25	98	0.1	{}	...	0.05	[0.5, 0.5]	rk4	1.0	diff_pg	1.0	1.0	0.000001	2	[-0.9095454221789239, -0.3493393384734713]

30 rows × 24 columns

4. Results processing and visualisation

In this section, we process the results of the optimisation and visualise them using various plots. The results include the fitness values, agent positions, and phase portraits.

# Recover the results from the optimiser
fp              = optimiser.results["fp"]
fg              = optimiser.results["fg"]
positions       = np.array(optimiser.results["p"])
best_position   = np.array(optimiser.results["g"])
v1              = np.array(optimiser.results["v1"])
v2              = np.array(optimiser.results["v2"])

# Calculate the absolute error in fitness values
efp             = np.abs(np.array(fp) - problem.optimum.y)
efg             = np.abs(np.array(fg) - problem.optimum.y)

# Convert the spikes to integer type
spikes          = np.array(optimiser.results["s"]).astype(int)

# Print some minimal information about the results
print(f"fg: {fg[-1][0]:.4f}, f*: {problem.optimum.y:.4f}, error: {efg[-1][0]:.4e}")
print(f"norm2(g - x*): {np.linalg.norm(best_position[-1] - problem.optimum.x):.4e}")
print(f"{v1.min():.4f} <= v1 <= {v1.max():.4f}")
print(f"{v2.min():.4f} <= v2 <= {v2.max():.4f}")

fg: -462.0900, f*: -462.0900, error: 3.6621e-06
norm2(g - x*): 1.8054e-04
-6.9711 <= v1 <= 8.0792
-6.1717 <= v2 <= 6.8563

Error convergence

This plot shows the convergence of the absolute error in fitness values over the iterations of the optimisation process.

fig, ax = plt.subplots(figsize=(width*0.9, height/1.8))

plt.plot(efp, color="silver", alpha=0.5)
plt.plot(np.max(efp, axis=1), '--', color="red", label=r"Max.")
plt.plot(np.average(efp, axis=1), '--', color="black", label=r"Mean")
plt.plot(np.median(efp, axis=1), '--', color="blue", label=r"Median")
plt.plot(efg, '--', color="green", label=r"Min.")

plt.xlabel(r"Step, $t$")
plt.ylabel(r"Abs. Error, $\varepsilon_f$")

lgd = plt.legend(ncol=2, loc="lower left")

plt.xscale("log")
plt.yscale("log")

ax.patch.set_alpha(0)
fig.tight_layout()

save_this(fig, _in="fitness")

Position evolution in 2D

This plot shows the evolution of the unit positions in the 2D space over the iterations of the optimisation process.

fig, ax = plt.subplots(figsize=(width/2, width/2))

cmap = plt.get_cmap('viridis', num_agents)
color = cmap(np.linspace(0, 1, num_agents))

for agent, c in enumerate(color):
    plt.plot(positions[:, agent, 0], positions[:,agent, 1], "--o",
             color=c, alpha=0.9, markersize=2, linewidth=1,
             label=f"Agent {agent}")

plt.plot(problem.optimum.x[0], problem.optimum.x[1], "*",
         color="black", markeredgecolor="cyan", markersize=12, label="Optimum")

plt.plot(best_position[:, 0], best_position[:, 1], "--*",
         color="red", markersize=3, label="Best position")

# plt.legend()
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")

ax.patch.set_alpha(0)
fig.tight_layout()

save_this(fig, _in="positions_2d")

Position evolution in 3D

This plot shows the evolution of the unit positions in the 3D space over the iterations of the optimisation process.

fig = plt.figure(figsize=(width/2, width/2))

ax = fig.add_subplot(111, projection='3d')
ax.set_proj_type('ortho')

cmap = plt.get_cmap('viridis', num_agents)
color = cmap(np.linspace(0, 1, num_agents))

steps = np.arange(optimiser._num_iterations) + 1

for agent, c in enumerate(color):
    ax.plot3D(positions[:, agent, 0], positions[:, agent, 1], steps,
              "--o", color=c, alpha=0.9, markersize=2, linewidth=1,
              label=f"Agent {agent}")

ax.plot3D(best_position[:, 0], best_position[:, 1], steps,
          "-", color="red", markersize=2,  linewidth=1.5,
          label="Best position")
ax.plot3D(problem.optimum.x[0], problem.optimum.x[1], steps,
          ":", color="black", markersize=2,  linewidth=1.5,
          label="Optimum")

for axis in [ax.xaxis, ax.yaxis, ax.zaxis]:
    axis.pane.set_edgecolor('black')
    axis.pane.set_linewidth(1.0)

# ax.viewfig, _init(elev=35, azim=135)
ax.view_init(elev=30, azim=100)
# ax.legend()
ax.set_xlabel(r"$x_1$", labelpad=1)
ax.set_ylabel(r"$x_2$", labelpad=1)
ax.set_zlabel(r"$t$", labelpad=0)
ax.set_box_aspect([1, 1, 0.8])

ax.set_zlim(1, num_steps)

ax.grid(False)
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False

ax.patch.set_alpha(0)
fig.subplots_adjust(left=0.1, right=0.95, top=0.95, bottom=0.1)
# fig.tight_layout()

save_this(fig, _in="positions_3d")

Phase portrait

This section visualises the phase portrait of the optimisation process in 2D. The phase portrait shows the trajectory of each agent in the 2D space over the iterations.

fig, axs = plt.subplots(nrows=np.ceil(num_dimensions / 2).astype(int),
                        ncols=2, figsize=(width, width/2))

steps = np.arange(optimiser._num_iterations) + 1
cmap = plt.get_cmap('Spectral', num_agents)
color = cmap(np.linspace(0, 1, num_agents))

axes = axs.flatten() if isinstance(axs, np.ndarray) else [axs]
for i, ax in enumerate(axes[:num_dimensions]):

    for agent, c in enumerate(color):
        ax.plot(v1[agent, :, i], v2[agent, :, i], "-o",
                color=c, alpha=0.5, markersize=1, linewidth=1,
                label=f"Agent {agent}")
        ax.plot(v1[agent, 0, i], v2[agent, 0, i], "-s",
                color=c, alpha=0.5, markersize=1, linewidth=1,
                label=f"Agent {agent}")

    ax.set_xlabel(r"$v_{}$".format("{1," + str(i+1) + "}"))
    ax.set_ylabel(r"$v_{}$".format("{2," + str(i+1) + "}"))

    ax.set_xlim(v1.min(), v1.max())
    ax.set_ylim(v2.min(), v2.max())

    ax.patch.set_alpha(0)

fig.tight_layout()

save_this(fig, _in="portrait_2d")

Phase portrait in 3D

This section visualises the phase portrait of the optimisation process in 3D. The phase portrait shows the trajectory of each agent in the 3D space over the iterations.

fig = plt.figure(figsize=(width, width/2))

num_rows = np.ceil(num_dimensions / 2).astype(int)
num_cols = 2

axs = [
    fig.add_subplot(
        num_rows, num_cols, i + 1, projection='3d'
    ) for i in range(num_dimensions)
]

steps = np.arange(optimiser._num_iterations) + 1
cmap = plt.get_cmap('Spectral', num_agents)
color = cmap(np.linspace(0, 1, num_agents))

for i, ax in enumerate(axs):
    ax.set_proj_type('ortho')
    ax.set_box_aspect([1, 1, 0.8])
    ax.view_init(elev=35, azim=110)

    for agent, c in enumerate(color):
        ax.plot3D(v1[agent, :, i], v2[agent, :, i], steps, "-",
                  color=c, alpha=0.9,
                  markersize=1, linewidth=1,
                  label=f"Agent {agent}")

        ax.plot3D(v1[agent, 0, i], v2[agent, 0, i], 0, "-s",
                  color=c, alpha=0.9,
                  markersize=1, linewidth=1,
                  label=f"Agent {agent}")

    ax.set_xlabel(r"$v_{}$".format("{" + str(i+1) + ",1}"))
    ax.set_ylabel(r"$v_{}$".format("{" + str(i+1) + ",2}"))
    ax.set_zlabel(r"$t$", labelpad=0.1)

    ax.set_xlim(v1.min(), v1.max())
    ax.set_ylim(v2.min(), v2.max())

    for axis in [ax.xaxis, ax.yaxis, ax.zaxis]:
        axis.pane.set_edgecolor('black')
        axis.pane.set_linewidth(1.0)

    ax.grid(False)
    ax.xaxis.pane.fill = False
    ax.yaxis.pane.fill = False
    ax.zaxis.pane.fill = False
    ax.patch.set_alpha(0)

fig.tight_layout()

save_this(fig, _in="portrait_3d")

Spike activity heatmap

This section visualises the spike activity of the agents over the iterations of the optimisation process. The heatmap shows the summed spike counts across all dimensions for each agent at each step.

fig, ax = plt.subplots(figsize=(width*0.9, height * 0.6))

spikes_sum = np.sum(spikes, axis=2)
cmap = plt.get_cmap("YlGnBu_r", 3)
im = ax.imshow(spikes_sum.T, aspect='auto', origin='lower',
               cmap=cmap, vmin=0, vmax=2)

ax.set_xlabel(r"Step, $t$")
ax.set_ylabel(r"Agent index, $i$")

cbar = fig.colorbar(im, ax=ax, pad=0.02, ticks=[0, 1, 2])
cbar.set_label(r"Spike count, $s_{1}+s_{2}$")

ax.patch.set_alpha(0)
fig.tight_layout()

save_this(fig, _in="spikes_hm")

window = max(5, num_steps // 50)  # choose a small window
# mean activity per step (across agents and dims)
rate = spikes.mean(axis=(1, 2)).astype(float)

# simple moving average
kernel = np.ones(window) / window
rate_ma = np.convolve(rate, kernel, mode='same')

fig, ax = plt.subplots(figsize=(width*0.9, height*0.4))
ax.plot(np.arange(num_steps), rate, alpha=0.4, label="Instantaneous")
ax.plot(np.arange(num_steps), rate_ma, linewidth=1.5, label=f"MA (w={window})")

ax.set_xlabel(r"Step, $t$")
ax.set_ylabel(r"Spike rate")
ax.legend(loc="upper right", ncol=2, frameon=False)

ax.patch.set_alpha(0)
fig.tight_layout()
save_this(fig, _in="spikes_rate_ma")