SMA Implementation
This notebook implements the slime mould algorithm as described in Slime Mould Algorithm. This implementation references the author’s MatLab code to resolve ambiguities.
This notebook uses the MIT license. Parts of the notebook transcribe (from MatLab to Python) work by Heidari under the MIT license. Heidari is also an author of the base paper.
[11]:
from typing import Literal as L, Callable
import numpy as np
import operator
import random
from numpy.random import default_rng
Set random seeds. Note that these seeds are also used in later cells, for example by default_rng.
[12]:
GLOBAL_SEED: int = 24601
np.random.seed(GLOBAL_SEED)
random.seed(GLOBAL_SEED)
Declare dimension of the problem: the size D of each solution and the number N of solutions in the population. Also declare related type aliases.
[13]:
#: Number of solutions in the population.
N: int
#: Size of each solution; or dimension of the problem.
D: int
type Location = np.ndarray[tuple[L['D']],
np.dtype[np.float64]]
"""Type of a solution in R^D. Location of "one"
slime mould in a D-dimensional space.
Alias for D-vectors.
"""
type Population = np.ndarray[tuple[L['N'], L['D']],
np.dtype[np.float64]]
"""Type of the population in R^{N\\times D}. Each ith
item is a Location.
"""
#: Fitness of one solution. Semantic alias for Float64
type Fitness = np.float64
type Fitnesses = np.ndarray[tuple[L['N']],
np.dtype[np.float64]]
"""Type of the vector of fitnesses. Contains one item
for each of the N solutions.
Alias for D-vectors.
"""
"""A hyperparameter for the algorithm. The canonical code
uses this value.
"""
MAX_ITERATION: int
Define the population initialiser. Parameters UB and LB limit the search space.
[14]:
def initialise_population(N: int,
dim: int,
UB: Location,
LB: Location)\
-> Population:
"""Given parameters, initialise a Population from
parameters.
"""
if np.any(UB<LB):
raise ValueError(
"For each index i, UB[i]"
"must not be less than LB[i].")
Xs = default_rng(GLOBAL_SEED).random((N, dim))
print(Xs*(UB - LB) + LB)
return Xs*(UB - LB) + LB
UB: Location
"""Upper bounds for items in a solution.
Xs[i][j] must not exceed UB[j].
"""
LB: Location
"""Lower bounds for items in a solution.
Xs[i][j] must not fall below LB[j].
""";
Initialise the fitness. Also define whether the task is to maximise or minimise the objective function.
[15]:
from dataclasses import dataclass
@dataclass
class SMAReport:
#: Solutions in each generation. 0th is earliest.
solutions: list[Population]
#: Fitnesses in each generation. 0th is earliest.
fitnesses: list[Fitnesses]
#: Best solution ever.
best_solution: Location
#: Best fitness ever.
best_fitness: Fitness
#: Best solutions in each generation. 0th is earliest.
gen_best_solutions: list[Location]
#: Best fitnesses in each generation. 0th is earliest.
gen_best_fitnesses: list[Fitness]
def sma(objective: Callable[[Location], Fitness],
is_maximising: bool,
max_iteration: int,
pop_ub: Location,
pop_lb: Location,
ub: Location,
lb: Location,
N: int,
D: int,
Z: float,) -> SMAReport:
historical_solutions: list[Population] = []
historical_fitnesses: list[Fitnesses] = []
gen_best_solutions: list[Location] = []
gen_best_fitnesses: list[Fitness] = []
Xs = initialise_population(N, D, pop_ub, pop_lb)
#: The best solution.
best_solution: Location
#: Fitness of the best solution.
best_fitness: Fitness
# Initialise :var:`best_fitness` as the worst
# possible fitness. If algorithm is maximising,
# then best_fitness starts small. Otherwise,
# best_fitness starts large.
if is_maximising:
best_fitness = np.float64(np.finfo(np.float64).min)
else:
best_fitness = np.float64(np.finfo(np.float64).max)
#: N-vector. S[i] is the fitness of the i^th solution.
S: Fitnesses
for t in range(max_iteration):
# Bound each Xs[i][j] to between UB[j] and LB[j]
for i in range(len(Xs)):
# Xs[i,:] is the ith slice
_ub_violation_mask = Xs[i,:] > ub
Xs[i,:] = np.multiply(Xs[i,:], ~_ub_violation_mask) +\
np.multiply(ub, _ub_violation_mask)
_lb_violation_mask = Xs[i,:] < lb
Xs[i,:] = np.multiply(Xs[i,:], ~_lb_violation_mask) +\
np.multiply(lb, _lb_violation_mask)
# Calculate fitnesses. S[i] is fitness(X[i,:])
S = np.apply_along_axis(objective, 1, Xs)
# Algorithm 1: Calculate the W by Eq. (2.5).
# small_dict is an ordered list. It maps
# each index `smell_index: int` to the
# smell_index^{nd} best solution.
# smell_indices: (dim,) of int
# smells: (dim,) of floats
# Implemented according to the canonical implementation.
# The best is at index 0; the worst is at index N.
# As such, sort according to :var:`IS_MAXIMISING`:
smell_dict = np.array(sorted(list(enumerate(S)), # type: ignore
key=lambda x: x[1],
reverse=True if is_maximising\
else False))
weights = np.zeros(shape=(N, D))
bF: np.float64 = smell_dict[0][1]
wF: np.float64 = smell_dict[-1][1]
# Partition smell_dict into two roughly-equal-sized
# matrices containing the better / worse halves.
smell_dict_condition = smell_dict[0:int(N/2)]
smell_dict_others = smell_dict[int(N/2):-1]
# This variable is named S in the canonical implementation.
# However, in the paper S means something different.
# This is name is to avoid ambiguity.
# `eps` is to avoid division by zero if bF=wF.
# The number is so infinitesimally small that we
# could probably ignore the probability that
# bF and wF differ by exactly eps.
S_eps = best_fitness - wF + np.finfo(np.float64).eps # was bF in paper ... maybe try to be more aggressive?
for j in range(0, D):
for smell_index, smell_order in smell_dict_condition:
weights[int(smell_index), j] = 1 + random.random()*np.log10((bF-smell_order)/(S_eps)+1)
for smell_index, smell_order in smell_dict_others:
weights[int(smell_index), j] = 1 - random.random()*np.log10((bF-smell_order)/(S_eps)+1)
if is_maximising:
is_better_than = operator.gt
else:
is_better_than = operator.lt
if is_better_than(bF, best_fitness):
# Recall that the 0th solution is best.
best_fitness = smell_dict[0][1]
best_solution = Xs[int(smell_dict[0][0])]
# Replacing t in original with t+1 here. This is because
# Matlab indices begin as 0.
a = np.arctanh(1-(t+1)/max_iteration) # Code says it's (2.4)
b = 1 - (t+1)/max_iteration
historical_solutions.append(Xs.copy())
historical_fitnesses.append(S.copy())
gen_best_solutions.append(Xs[int(smell_dict[0][0])].copy())
gen_best_fitnesses.append(smell_dict[0][1].copy())
for i in range(N):
rand=random.random()
# The code ... uses rand directly?
if rand < Z: # Code says it's (2.7)
Xs[i,:] = (pop_ub-pop_lb)*random.random() + pop_lb
else:
# AllFitness in code is S here
p = np.tanh(np.abs(S[i] - best_fitness))
vb = np.random.uniform(low=-a, high=a, size=(D,)) # Code says it's 2.3
vc = np.random.uniform(low=-b, high=b, size=(D,))
for j in range(D):
r = random.random()
A = random.randint(0, N-1)
B = random.randint(0, N-1)
if r < p:
Xs[i, j] = best_solution[j] + vb[j] * (weights[i, j]*Xs[A,j]-Xs[B,j])
else:
Xs[i, j] = vc[j] * Xs[i,j]
result: SMAReport = SMAReport(
solutions=historical_solutions,
fitnesses=historical_fitnesses,
best_solution=best_solution,
best_fitness=best_fitness,
gen_best_solutions=gen_best_solutions,
gen_best_fitnesses=gen_best_fitnesses,
)
return result
Parameterise and run the algorithm.
[16]:
from ograph import ofunc, oplot
[ ]:
[17]:
N = 200
D = 2
POP_UB = np.ones(shape=(D,)) * -1
POP_LB = np.ones(shape=(D,)) * -4
UB = np.ones(shape=(D,)) * 7
LB = np.ones(shape=(D,)) * -7
IS_MAXIMISING: bool = True
Z = 0.03
MAX_ITERATION = 200
#: The objective function
OBJECTIVE: Callable[[Location], Fitness] =\
lambda x: ofunc.himmelblau(x[0],
x[1]) * -1\
if IS_MAXIMISING else 1 # type: ignore
report: SMAReport = sma(objective = OBJECTIVE,
is_maximising = IS_MAXIMISING,
max_iteration = MAX_ITERATION,
pop_ub=POP_UB,
pop_lb=POP_LB,
ub=UB,
lb=LB,
N=N,
D=D,
Z=Z)
[[-2.47592429 -1.81997011]
[-2.95540877 -3.33771139]
[-2.74391284 -1.82419344]
[-3.70899605 -1.68464243]
[-3.12115934 -1.37014237]
[-2.53291771 -3.64338561]
[-1.72119328 -1.28885633]
[-2.33623507 -1.08980774]
[-3.26921986 -2.55601543]
[-3.12279309 -2.57700669]
[-3.38698749 -1.15188707]
[-2.05444234 -1.46576723]
[-1.60774563 -3.08985395]
[-1.35949565 -1.44753345]
[-2.03237569 -2.18611263]
[-1.91113253 -3.40816876]
[-3.45566262 -3.36016751]
[-3.61281961 -2.52430802]
[-1.02348297 -2.21315978]
[-1.92142152 -2.66819629]
[-2.39064646 -2.91871284]
[-1.23466717 -2.51284576]
[-2.55051423 -3.798853 ]
[-3.95026681 -3.83633743]
[-3.96621018 -1.98419591]
[-3.50945244 -1.27658987]
[-3.22636507 -2.85088508]
[-2.44282669 -2.4260564 ]
[-1.53628828 -2.69364699]
[-2.36810088 -1.00410217]
[-3.11530769 -2.05087961]
[-2.05258559 -3.36806285]
[-3.53886252 -3.9694065 ]
[-1.3621931 -1.16921042]
[-1.696319 -2.74432217]
[-3.81658648 -2.77322841]
[-3.93657671 -1.2421323 ]
[-1.34358889 -3.06222432]
[-3.76172238 -3.62649093]
[-1.38924747 -2.12576489]
[-1.80712843 -2.31343234]
[-3.40268273 -2.27668321]
[-1.16427419 -2.73559948]
[-1.93870581 -1.73864863]
[-2.37392866 -2.75147314]
[-3.87192287 -2.24350892]
[-2.34202431 -1.41521903]
[-1.03185307 -1.42448362]
[-2.27137564 -1.49339086]
[-2.5378265 -3.74476359]
[-2.10174376 -3.41777756]
[-3.18254033 -1.21881549]
[-1.32006398 -3.64269858]
[-1.91151534 -1.08159944]
[-1.47428545 -3.52296 ]
[-2.79978241 -3.83048691]
[-2.73554389 -2.50570139]
[-2.60446389 -1.32918957]
[-2.44779074 -2.661567 ]
[-1.29374689 -2.51676992]
[-1.47920299 -2.18365621]
[-2.16098851 -1.57941823]
[-2.82107172 -2.55913883]
[-3.59611955 -3.69247776]
[-3.28155639 -2.52486625]
[-3.35494193 -1.91654424]
[-3.16305946 -1.7731639 ]
[-2.22308459 -1.69209623]
[-1.09012889 -1.71263931]
[-3.88693587 -3.09800484]
[-3.3473536 -2.84332726]
[-3.73451444 -3.01135429]
[-3.76125346 -2.45202997]
[-1.35278209 -3.56194624]
[-1.63793391 -3.17910374]
[-2.9649655 -2.84549329]
[-1.11609482 -3.01301184]
[-1.2262371 -2.9566359 ]
[-2.28486943 -2.06914252]
[-3.7267008 -1.32296343]
[-2.56948376 -1.09182448]
[-3.98085922 -1.74046479]
[-2.91341855 -1.56751214]
[-1.33330165 -2.74246792]
[-3.26042757 -3.32370859]
[-3.47616085 -3.90004547]
[-2.08280041 -2.98357382]
[-2.47823356 -2.89552476]
[-1.7924004 -2.40325816]
[-3.31049065 -3.8383048 ]
[-1.53415311 -1.37880931]
[-3.52002152 -2.45316001]
[-3.93736291 -3.62328251]
[-1.24878303 -1.17312098]
[-3.62471116 -1.49570033]
[-2.80068436 -1.13425278]
[-3.2638772 -2.18697822]
[-3.45478549 -3.18581588]
[-3.62457202 -3.98949789]
[-3.68974153 -2.13974449]
[-3.24142622 -1.37405943]
[-2.66465865 -2.32745027]
[-3.64719944 -3.91604168]
[-1.72132094 -3.48397851]
[-2.27790202 -1.19762802]
[-1.95870719 -2.47135787]
[-2.99962105 -2.56463065]
[-3.34058563 -2.95747866]
[-1.07865836 -2.86619631]
[-2.86550952 -3.51947678]
[-3.15997956 -1.31631882]
[-3.81358945 -3.47630911]
[-1.80326271 -2.84731933]
[-3.81049119 -3.3414486 ]
[-2.88653917 -2.26224947]
[-3.75765269 -2.03059638]
[-1.69842732 -1.67950422]
[-1.10621879 -1.91947618]
[-2.59497986 -1.20477833]
[-2.53341649 -2.10843739]
[-2.39596934 -1.01293467]
[-2.05177285 -2.77754899]
[-1.24294954 -3.30699857]
[-3.74531197 -2.48200374]
[-1.81253217 -3.13723949]
[-2.93803093 -1.93859849]
[-3.24906568 -2.76776741]
[-2.18681072 -2.05932403]
[-1.83690268 -1.54028669]
[-2.31175863 -1.71190059]
[-1.18953702 -1.93295709]
[-2.21999599 -1.7823127 ]
[-2.38500489 -1.82205898]
[-1.4252258 -2.58258284]
[-2.42025759 -3.91852256]
[-3.23610232 -2.62886081]
[-3.62572128 -2.27520051]
[-2.32607002 -1.25987613]
[-1.46239442 -2.77781081]
[-1.025142 -3.56547457]
[-1.26883583 -2.30054105]
[-1.23096913 -3.72684842]
[-2.77318501 -2.79151009]
[-2.33445358 -3.19602229]
[-2.96505354 -2.30971645]
[-2.41090855 -2.54689573]
[-1.32452311 -1.11226622]
[-2.53921083 -1.25020504]
[-2.36262288 -3.2611987 ]
[-1.69692747 -1.2936879 ]
[-3.43556928 -1.87757262]
[-3.95422027 -2.40760427]
[-1.10131194 -1.90829358]
[-1.80428306 -3.2513341 ]
[-3.44366843 -2.2996864 ]
[-3.62591078 -3.70612528]
[-2.74077055 -2.77437767]
[-3.66911411 -3.26080062]
[-3.72013603 -3.3772269 ]
[-3.81412233 -2.94955127]
[-2.3324628 -1.47515259]
[-2.22100761 -2.58766684]
[-2.07797463 -3.51482905]
[-1.47731733 -3.6234165 ]
[-3.47019464 -3.63316133]
[-2.90143483 -1.99771576]
[-1.75228291 -1.39501959]
[-3.8114033 -3.68825511]
[-1.20786874 -3.96460392]
[-2.83593277 -2.3252146 ]
[-2.62628278 -1.50187703]
[-3.29921907 -1.23407264]
[-2.32970278 -3.15306616]
[-2.30941998 -1.75723721]
[-3.1245127 -3.63280641]
[-1.05138074 -1.96040215]
[-2.36792012 -3.23299844]
[-1.42315267 -2.35999401]
[-3.47842811 -1.36905547]
[-3.11673082 -2.93499453]
[-3.63518322 -3.97754817]
[-3.3953347 -1.47998179]
[-1.15114295 -3.86829738]
[-2.28704496 -1.44676427]
[-1.16091465 -1.79337725]
[-1.85292888 -1.15332267]
[-2.76692174 -2.68306539]
[-2.12602633 -2.85424294]
[-1.34284003 -3.05060179]
[-1.78918981 -1.7868207 ]
[-1.87361448 -1.09093413]
[-3.10341507 -3.53840054]
[-1.99017561 -2.84100373]
[-2.83066363 -2.27584255]
[-2.12206382 -1.15996083]
[-2.66494694 -1.04418838]
[-1.42596273 -1.30656911]
[-1.45538376 -1.81875416]
[-2.7928368 -1.27597375]
[-3.14304053 -2.39254982]]
Visualise
[18]:
print(f"Best solution is {report.best_solution}, with fitness {report.best_fitness}")
Best solution is [-3.31078708 -3.3121009 ], with fitness -0.032185163502593615
[ ]:
from ograph.swarm import plot_fitnesses
plot_fitnesses(report.gen_best_fitnesses,
best_selector=max)
[20]:
plot_bests(report.solutions, report.gen_best_fitnesses)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 plot_bests(report.solutions, report.gen_best_fitnesses)
NameError: name 'plot_bests' is not defined
Visualise Particles
[ ]:
from ograph.swarm import plot_positions, animate_positions
[ ]:
plot_positions(report.solutions,
objective=lambda x, y: OBJECTIVE((x, y)))
Because the ipython notebook might not correctly render gifs, the animation will be saved to the working directory.
[ ]:
animate_positions(report.solutions,
objective=lambda x, y: OBJECTIVE((x, y)))