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       = 100       # Number of iterations
problem_id      = 1         # 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  = 10        # 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": "fixed",
     "hs_operator": "differential",
     "hs_variant": "current-to-rand",
     },
]

# 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_1p_1i_2d_100s_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 1. Sphere (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 100 iterations...
... step: 0, best fitness: 81.21773529052734
... step: 10, best fitness: 79.9867172241211
... step: 20, best fitness: 79.5005111694336
... step: 30, best fitness: 79.5005111694336
... step: 40, best fitness: 79.48989868164062
... step: 50, best fitness: 79.4843978881836
... step: 60, best fitness: 79.4843978881836
... step: 70, best fitness: 79.4843978881836
... step: 80, best fitness: 79.48188018798828
... step: 90, best fitness: 79.48102569580078
... step: 99, best fitness: 79.48017120361328
[neuropt:log] Simulation completed. Fetching monitor data... done

(array([ 0.24347829, -1.14769556]), array([79.4801712]))

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

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

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

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

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: 79.4802, f*: 79.4800, error: 1.7120e-04
norm2(g - x*): 1.3030e-02
-1.9109 <= v1 <= 2.6796
-2.1922 <= v2 <= 2.8139

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))

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

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(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))

for i, ax in enumerate(axs.flatten()):

    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")