Split objective (#8)

* First implementation of split

* Add tests

* Update test

* Add objectives to init

* Add evaluation function to objectives

* Add objective class support

optimizer/__init__.py

@@ -4,3 +4,4 @@
 from .file_manager import FileManager
 from .optimizer import Optimizer
 from .mopso import MOPSO
+from .objective import Objective, ElementWiseObjective, BatchObjective

optimizer/mopso.py

@@ -103,17 +103,6 @@ def set_state(self, velocity, position, best_position, fitness, best_fitness):
         self.fitness = fitness
         self.best_fitness = best_fitness
 
-    def evaluate_fitness(self, objective_functions):
-        """
-        Evaluate the fitness of the particle based on the provided objective functions.
-        Calls the `set_fitness` method to update the particle's fitness and best position.
-
-        Parameters:
-            objective_functions (list): List of objective functions used for fitness evaluation.
-        """
-        fitness = [obj_func(self.position) for obj_func in objective_functions]
-        self.set_fitness(fitness)
-
     def update_best(self):
         """
         Update particle's fitness and best position
@@ -201,26 +190,21 @@ class MOPSO(Optimizer):
             Calculate the crowding distance for particles in the Pareto front.
     """
 
-    def __init__(self, objective_functions,
+    def __init__(self,
+                 objective,
                  lower_bounds, upper_bounds, num_particles=50,
                  inertia_weight=0.5, cognitive_coefficient=1, social_coefficient=1,
                  num_iterations=100,
                  optimization_mode='individual',
                  max_iter_no_improv=None,
-                 num_objectives=None,
                  incremental_pareto=False):
-        self.objective_functions = objective_functions
+        self.objective = objective
         if FileManager.loading_enabled:
             try:
                 self.load_checkpoint(num_additional_iterations=num_iterations)
                 return
             except FileNotFoundError as e:
                 print("Checkpoint not found. Fallback to standard construction.")
-
-        if num_objectives is None:
-            self.num_objectives = len(self.objective_functions)
-        else:
-            self.num_objectives = num_objectives
         self.num_particles = num_particles
         self.num_params = len(lower_bounds)
         self.lower_bounds = lower_bounds
@@ -231,7 +215,7 @@ def __init__(self, objective_functions,
         self.num_iterations = num_iterations
         self.max_iter_no_improv = max_iter_no_improv
         self.optimization_mode = optimization_mode
-        self.particles = [Particle(lower_bounds, upper_bounds, num_objectives, num_particles)
+        self.particles = [Particle(lower_bounds, upper_bounds, objective.num_objectives, num_particles)
                           for _ in range(num_particles)]
         self.iteration = 0
         self.incremental_pareto = incremental_pareto
@@ -248,7 +232,6 @@ def save_attributes(self):
         pso_attributes = {
             'lower_bounds': self.lower_bounds,
             'upper_bounds': self.upper_bounds,
-            'num_objectives': self.num_objectives,
             'num_particles': self.num_particles,
             'num_params': self.num_params,
             'inertia_weight': self.inertia_weight,
@@ -305,7 +288,6 @@ def load_checkpoint(self, num_additional_iterations):
         # restore pso attributes
         self.lower_bounds = pso_attributes['lower_bounds']
         self.upper_bounds = pso_attributes['upper_bounds']
-        self.num_objectives = pso_attributes['num_objectives']
         self.num_particles = pso_attributes['num_particles']
         self.num_params = pso_attributes['num_params']
         self.inertia_weight = pso_attributes['inertia_weight']
@@ -321,7 +303,7 @@ def load_checkpoint(self, num_additional_iterations):
         self.particles = []
         for i in range(self.num_particles):
             particle = Particle(self.lower_bounds, self.upper_bounds,
-                                num_objectives=self.num_objectives,
+                                num_objectives=self.objective.num_objectives,
                                 num_particles=self.num_particles)
             particle.set_state(
                 position=np.array(
@@ -330,7 +312,7 @@ def load_checkpoint(self, num_additional_iterations):
                     individual_states[i][self.num_params:2*self.num_params], dtype=float),
                 best_position=np.array(
                     individual_states[i][2*self.num_params:3*self.num_params], dtype=float),
-                fitness=[np.inf] * self.num_objectives,
+                fitness=[np.inf] * self.objective.num_objectives,
                 best_fitness=np.array(
                     individual_states[i][3*self.num_params:], dtype=float)
             )
@@ -340,7 +322,7 @@ def load_checkpoint(self, num_additional_iterations):
         self.pareto_front = []
         for i in range(len(pareto_front)):
             particle = Particle(self.lower_bounds, self.upper_bounds,
-                                num_objectives=self.num_objectives,
+                                num_objectives=self.objective.num_objectives,
                                 num_particles=self.num_particles)
             particle.set_state(position=pareto_front[i][:self.num_params],
                                fitness=pareto_front[i][self.num_params:],
@@ -364,16 +346,9 @@ def optimize(self):
             list: List of Particle objects representing the Pareto front of non-dominated solutions.
         """
         for _ in range(self.num_iterations):
-            if self.optimization_mode == 'global':
-                optimization_output = [objective_function([particle.position for
-                                                           particle in self.particles])
-                                       for objective_function in self.objective_functions]
-            for p_id, particle in enumerate(self.particles):
-                if self.optimization_mode == 'individual':
-                    particle.evaluate_fitness(self.objective_functions)
-                if self.optimization_mode == 'global':
-                    particle.set_fitness([output[p_id]
-                                         for output in optimization_output])
+            optimization_output = self.objective.evaluate([particle.position for particle in self.particles])
+            [particle.set_fitness(optimization_output[:,p_id]) for p_id, particle in enumerate(self.particles)]
+
             FileManager.save_csv([np.concatenate([particle.position, np.ravel(
                                  particle.fitness)]) for particle in self.particles],
                                  'history/iteration' + str(self.iteration) + '.csv')

optimizer/objective.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+class Objective():
+    def __init__(self, objective_functions, num_objectives = None) -> None:
+        self.objective_functions = objective_functions
+        if num_objectives is None:
+            self.num_objectives = len(self.objective_functions)
+        else:
+            self.num_objectives = num_objectives
+        pass
+
+    def evaluate(self, items):
+        return np.array([objective_function([item for item in items]) for objective_function in self.objective_functions])
+
+    def type(self):
+        return self.__class__.__name__
+
+class ElementWiseObjective(Objective):
+    def __init__(self, objective_functions, num_objectives=None) -> None:
+        super().__init__(objective_functions, num_objectives)
+
+    def evaluate(self, items):
+        return np.array([[obj_func(item) for item in items] for obj_func in self.objective_functions])
+
+class BatchObjective(Objective):
+    def __init__(self, objective_functions, num_objectives=None) -> None:
+        super().__init__(objective_functions, num_objectives)

tests/schaffer.py

@@ -24,8 +24,10 @@ def f(params):
 optimizer.FileManager.working_dir="tmp/schaffer/"
 optimizer.FileManager.loading_enabled = True
 
-pso = optimizer.MOPSO(objective_functions=[f],lower_bounds=lb, upper_bounds=ub, 
-            num_objectives=2, num_particles=num_agents, num_iterations=num_iterations, 
+objective = optimizer.Objective([f])
+
+pso = optimizer.MOPSO(objective=objective,lower_bounds=lb, upper_bounds=ub, 
+            num_particles=num_agents, num_iterations=num_iterations, 
             inertia_weight=0.5, cognitive_coefficient=1, social_coefficient=1, 
             max_iter_no_improv=None, optimization_mode='global')
 

tests/zdt1.py

@@ -0,0 +1,51 @@
+import optimizer
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+
+
+num_agents = 100
+num_iterations = 200
+num_params = 30
+
+lb = [0] * num_params
+ub = [1] * num_params
+
+def zdt1_objective1(x):
+    return x[0]
+
+def zdt1_objective2(x):
+    f1 = x[0]
+    g = 1 + 9.0 / (len(x)-1) * sum(x[1:])
+    h = 1.0 - np.sqrt(f1 / g)
+    f2 = g * h
+    return f2
+
+optimizer.FileManager.working_dir="tmp/zdt1/"
+optimizer.FileManager.loading_enabled = False
+optimizer.FileManager.saving_enabled = False
+
+objective = optimizer.ElementWiseObjective([zdt1_objective1, zdt1_objective2])
+
+pso = optimizer.MOPSO(objective=objective,lower_bounds=lb, upper_bounds=ub, 
+            num_particles=num_agents, num_iterations=num_iterations, 
+            inertia_weight=0.6, cognitive_coefficient=1, social_coefficient=2, 
+            max_iter_no_improv=None)
+
+# run the optimization algorithm
+pso.optimize()
+
+fig, ax = plt.subplots()
+
+pareto_front = pso.get_current_pareto_front()
+n_pareto_points = len(pareto_front)
+pareto_x = [particle.fitness[0] for particle in pareto_front]
+pareto_y = [particle.fitness[1] for particle in pareto_front]
+real_x = (np.linspace(0, 1, n_pareto_points))
+real_y = 1-np.sqrt(real_x)
+plt.scatter(real_x, real_y, s=5, c='red')
+plt.scatter(pareto_x, pareto_y, s=5)
+
+plt.savefig('tmp/pf.png')
+

tests/zdt2.py

@@ -0,0 +1,52 @@
+import optimizer
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+
+
+num_agents = 100
+num_iterations = 300
+num_params = 30
+
+lb = [0] * num_params
+ub = [1] * num_params
+
+def zdt2_objective1(x):
+    return x[0]
+
+def zdt2_objective2(x):
+    f1 = x[0]
+    g = 1.0 + 9.0 * sum(x[1:]) / (len(x) - 1)
+    h = 1.0 - np.power((f1 *1.0 / g),2)
+    f2 = g * h
+    return f2
+
+optimizer.FileManager.working_dir="tmp/zdt2/"
+optimizer.FileManager.loading_enabled = False
+optimizer.FileManager.saving_enabled = True
+
+objective = optimizer.ElementWiseObjective([zdt2_objective1, zdt2_objective2])
+
+pso = optimizer.MOPSO(objective=objective,lower_bounds=lb, upper_bounds=ub, 
+            num_particles=num_agents, num_iterations=num_iterations, 
+            inertia_weight=0.4, cognitive_coefficient=0, social_coefficient=2, 
+            max_iter_no_improv=None)
+
+# run the optimization algorithm
+pso.optimize()
+
+fig, ax = plt.subplots()
+
+pareto_front = pso.get_current_pareto_front()
+n_pareto_points = len(pareto_front)
+pareto_x = [particle.fitness[0] for particle in pareto_front]
+pareto_y = [particle.fitness[1] for particle in pareto_front]
+
+real_x = (np.linspace(0, 1, n_pareto_points))
+real_y = 1 - np.power(real_x, 2)
+plt.scatter(real_x, real_y, s=5, c='red')
+plt.scatter(pareto_x, pareto_y, s=5)
+
+plt.savefig('tmp/pf.png')
+

tests/zdt3.py

@@ -0,0 +1,66 @@
+import optimizer
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+
+
+num_agents = 200
+num_iterations = 300
+num_params = 30
+
+lb = [0] * num_params
+ub = [1] * num_params
+
+def zdt3_objective1(x):
+    return x[0]
+
+def zdt3_objective2(x):
+    f1 = x[0]
+    g = 1.0 + 9.0 * sum(x[1:]) / (len(x) - 1)
+    h = (1.0 - np.power(f1 * 1.0 / g, 0.5) - (f1 * 1.0 / g) * np.sin(10 * np.pi * f1))
+    f2 = g * h
+    return f2
+
+optimizer.FileManager.working_dir="tmp/zdt3/"
+optimizer.FileManager.loading_enabled = False
+optimizer.FileManager.saving_enabled = True
+
+objective = optimizer.ElementWiseObjective([zdt3_objective1, zdt3_objective2])
+
+pso = optimizer.MOPSO(objective=objective,lower_bounds=lb, upper_bounds=ub, 
+            num_particles=num_agents, num_iterations=num_iterations, 
+            inertia_weight=0.4, cognitive_coefficient=1, social_coefficient=2, 
+            max_iter_no_improv=None)
+
+# run the optimization algorithm
+pso.optimize()
+
+fig, ax = plt.subplots()
+
+pareto_front = pso.get_current_pareto_front()
+n_pareto_points = len(pareto_front)
+pareto_x = [particle.fitness[0] for particle in pareto_front]
+pareto_y = [particle.fitness[1] for particle in pareto_front]
+
+regions = [[0, 0.0830015349],
+                   [0.182228780, 0.2577623634],
+                   [0.4093136748, 0.4538821041],
+                   [0.6183967944, 0.6525117038],
+                   [0.8233317983, 0.8518328654]]
+
+pf = []
+
+for r in regions:
+    x1 = np.linspace(r[0], r[1], int(n_pareto_points / len(regions)))
+    x2 = 1 - np.sqrt(x1) - x1 * np.sin(10 * np.pi * x1)
+    pf.append([x1, x2])
+
+real_x = np.concatenate([x for x, _ in pf])
+real_y = np.concatenate([y for _, y in pf])
+
+plt.scatter(real_x, real_y, s=5, c='red')
+plt.scatter(pareto_x, pareto_y, s=5)
+
+plt.savefig('tmp/pf.png')
+