From baf6e40e152470b3997ca1a3af474b5478a4c570 Mon Sep 17 00:00:00 2001
From: coolneng <akasroua@gmail.com>
Date: Tue, 22 Jun 2021 09:58:34 +0200
Subject: [PATCH] Implement simulated annealing

---
 src/simulated_annealing.py | 111 +++++++++++++++++++++++++++++++++++++
 1 file changed, 111 insertions(+)
 create mode 100644 src/simulated_annealing.py

diff --git a/src/simulated_annealing.py b/src/simulated_annealing.py
new file mode 100644
index 0000000..cfd1926
--- /dev/null
+++ b/src/simulated_annealing.py
@@ -0,0 +1,111 @@
+from numpy.random import choice, randint, random
+from pandas import DataFrame
+from itertools import combinations
+from math import exp, log
+
+
+def get_row_distance(source, destination, data):
+    row = data.query(
+        """(source == @source and destination == @destination) or \
+        (source == @destination and destination == @source)"""
+    )
+    return row["distance"].values[0]
+
+
+def compute_distance(element, solution, data):
+    accumulator = 0
+    distinct_elements = solution.query(f"point != {element}")
+    for _, item in distinct_elements.iterrows():
+        accumulator += get_row_distance(
+            source=element, destination=item.point, data=data
+        )
+    return accumulator
+
+
+def generate_first_solution(n, m, data):
+    solution = DataFrame(columns=["point", "distance", "fitness"])
+    solution["point"] = choice(n, size=m, replace=False)
+    solution["distance"] = solution["point"].apply(
+        func=compute_distance, solution=solution, data=data
+    )
+    solution = evaluate_solution(solution, data)
+    return solution
+
+
+def evaluate_solution(solution, data):
+    fitness = 0
+    comb = combinations(solution.index, r=2)
+    for index in list(comb):
+        elements = solution.loc[index, :]
+        fitness += get_row_distance(
+            source=elements["point"].head(n=1).values[0],
+            destination=elements["point"].tail(n=1).values[0],
+            data=data,
+        )
+    solution["fitness"] = fitness
+    return solution
+
+
+def element_in_dataframe(solution, element):
+    duplicates = solution.query(f"point == {element}")
+    return not duplicates.empty
+
+
+def generate_neighbour(previous, n, data):
+    solution = previous.copy()
+    random_index = randint(len(solution.point.values))
+    random_element = randint(n)
+    while element_in_dataframe(solution=solution, element=random_element):
+        random_element = randint(n)
+    solution["point"].loc[random_index] = random_element
+    solution["distance"].loc[random_index] = compute_distance(
+        element=solution["point"].loc[random_index], solution=solution, data=data
+    )
+    solution = evaluate_solution(solution, data)
+    return solution
+
+
+def get_temperatures(initial_fitness, mu=0.3, phi=0.3):
+    T0 = (mu * initial_fitness) / -(log(phi))
+    Tf = 1e-3
+    return T0, Tf
+
+
+def get_metaparameters(n, max_iterations):
+    max_neighbours = 10 * n
+    max_successes = 0.1 * max_neighbours
+    M = max_iterations / max_neighbours
+    return max_neighbours, max_successes, M
+
+
+def beta(T0, Tf, M):
+    return (T0 - Tf) / (M * T0 * Tf)
+
+
+def cool_down(T, T0, Tf, M):
+    return T / (1 + beta(T0, Tf, M) * T)
+
+
+def simulated_annealing(n, m, data, max_iterations=10000):
+    initial_solution = generate_first_solution(n, m, data)
+    T0, Tf = get_temperatures(initial_fitness=initial_solution.fitness.values[0])
+    max_neighbours, max_successes, M = get_metaparameters(n, max_iterations)
+    current_solution = initial_solution
+    best_solution = initial_solution
+    T = T0
+    while T > T0:
+        num_successes = 0
+        for _ in range(max_neighbours):
+            neighbour = generate_neighbour(current_solution, n, data)
+            delta = neighbour.fitness.values[0] - current_solution.fitness.values[0]
+            if delta > 0 or random() < exp((-delta) / T):
+                current_sol = neighbour
+                num_successes += 1
+                if current_solution.fitness.values[0] > best_solution.fitness.values[0]:
+                    best_solution = current_sol
+            if num_successes >= max_successes:
+                break
+        if num_successes == 0:
+            break
+        T = cool_down(T, T0, Tf, M)
+    return best_solution