MOEAD

.gitignore

@@ -10,8 +10,6 @@ __pycache__/
 *.so
 
 
-pymoo/experimental/*
-
 
 pymoo/cython/*html
 

pymoo/algorithms/moead.py

@@ -1,30 +1,34 @@
 import numpy as np
 from scipy.spatial.distance import cdist
 
+from pymoo.cython.decomposition import cython_pbi
 from pymoo.algorithms.genetic_algorithm import GeneticAlgorithm
-from pymoo.operators.crossover.real_differental_evolution_crossover import DifferentialEvolutionCrossover
+from pymoo.operators.crossover.real_simulated_binary_crossover import SimulatedBinaryCrossover
 from pymoo.operators.default_operators import set_if_none
+from pymoo.operators.mutation.real_polynomial_mutation import PolynomialMutation
 from pymoo.operators.sampling.real_random_sampling import RealRandomSampling
 from pymoo.rand import random
-from pymoo.util.decomposition import tchebi
+from pymoo.util.decomposition import decompose
 from pymoo.util.display import disp_multi_objective
-from pymop.util import get_uniform_weights
 
 
 class MOEAD(GeneticAlgorithm):
     def __init__(self,
-                 pop_size=100,
+                 ref_dirs,
                  n_neighbors=15,
-                 prob=0.5,
-                 weight=0.8,
-                 method="binomial",
+                 decomposition="cython_pbi",
+                 crossover="sbx",
+                 prob_neighbor_mating=0.7,
                  **kwargs):
 
         self.n_neighbors = n_neighbors
+        self.prob_neighbor_mating = prob_neighbor_mating
+        self.decomposition = decomposition
 
-        set_if_none(kwargs, 'pop_size', pop_size)
+        set_if_none(kwargs, 'pop_size', len(ref_dirs))
         set_if_none(kwargs, 'sampling', RealRandomSampling())
-        set_if_none(kwargs, 'crossover', DifferentialEvolutionCrossover(prob=prob, weight=weight, method=method))
+        set_if_none(kwargs, 'crossover', SimulatedBinaryCrossover(prob_cross=0.9, eta_cross=20))
+        set_if_none(kwargs, 'mutation', PolynomialMutation(eta_mut=15))
         set_if_none(kwargs, 'selection', None)
         set_if_none(kwargs, 'mutation', None)
         set_if_none(kwargs, 'survival', None)
@@ -34,52 +38,67 @@ def __init__(self,
         self.func_display_attrs = disp_multi_objective
 
         # initialized when problem is known
-        self.weights = None
-        self.neighbours = None
-        self.ideal_point = None
-
-    def _initialize(self):
-
-        # weights to be used for decomposition
-        self.weights = get_uniform_weights(self.pop_size, self.problem.n_obj)
+        self.ref_dirs = ref_dirs
 
         # neighbours includes the entry by itself intentionally for the survival method
-        self.neighbours = np.argsort(cdist(self.weights, self.weights), axis=1, kind='quicksort')[:, :self.n_neighbors + 1]
+        self.neighbors = np.argsort(cdist(self.ref_dirs, self.ref_dirs), axis=1, kind='quicksort')[:, :self.n_neighbors]
 
-        # set the initial ideal point
+        self.ideal_point = None
+
+    def _initialize(self):
         pop = super()._initialize()
         self.ideal_point = np.min(pop.F, axis=0)
-
         return pop
 
     def _next(self, pop):
 
-        # iterate for each member of the population
-        for i in range(self.pop_size):
+        # iterate for each member of the population in random order
+        for i in random.perm(self.pop_size):
 
-            # all neighbors shuffled (excluding the individual itself)
-            neighbors = self.neighbours[i][1:][random.perm(self.n_neighbors - 1)]
-            parents = np.concatenate([[i], neighbors[:self.crossover.n_parents - 1]])
+            if random.random() < self.prob_neighbor_mating:
+                parents = self.neighbors[i, :][random.perm(self.n_neighbors)][:self.crossover.n_parents]
+            else:
+                parents = np.arange(self.pop_size)[random.perm(self.pop_size)][:self.crossover.n_parents]
 
             # do recombination and create an offspring
-            X = self.crossover.do(self.problem, pop.X[None, parents, :], X=pop.X[[i], :])
+            X = self.crossover.do(self.problem, pop.X[parents, None, :])
+            X = self.mutation.do(self.problem, X)
+            X = X[0, :]
 
             # evaluate the offspring
             F, _ = self.evaluator.eval(self.problem, X)
 
             # update the ideal point
-            self.ideal_point = np.min(np.concatenate([self.ideal_point[None, :], F], axis=0), axis=0)
+            self.ideal_point = np.min(np.vstack([self.ideal_point, F]), axis=0)
 
-            # for each offspring that was created
-            for k in range(self.crossover.n_children):
-                # the weights of each neighbor
-                weights_of_neighbors = self.weights[self.neighbours[i], :]
+            batch = False
+
+            if batch:
+
+                # all neighbors of this individual and corresponding weights
+                N = self.neighbors[i, :]
+                w = self.ref_dirs[N, :]
 
                 # calculate the decomposed values for each neighbour
-                FV = tchebi(pop.F[self.neighbours[i]], weights_of_neighbors, self.ideal_point)
-                off_FV = tchebi(F[[k], :], weights_of_neighbors, self.ideal_point)
+                FV = decompose(pop.F[N, :], w, method=self.decomposition, ideal_point=self.ideal_point)[0]
+                off_FV = decompose(F[None, :], w, method=self.decomposition, ideal_point=self.ideal_point)[0]
 
                 # get the absolute index in F where offspring is better than the current F (decomposed space)
-                off_is_better = self.neighbours[i][np.where(off_FV < FV)[0]]
-                pop.F[off_is_better, :] = F[k, :]
-                pop.X[off_is_better, :] = X[k, :]
+                neighbors_to_update = N[off_FV < FV]
+                pop.F[neighbors_to_update, :] = F
+                pop.X[neighbors_to_update, :] = X
+
+            else:
+
+                for neighbor in self.neighbors[i][random.perm(self.n_neighbors)]:
+
+                    # the weight vector of this neighbor
+                    w = self.ref_dirs[neighbor, :]
+
+                    FV = decompose(pop.F[[neighbor], :], w, method=self.decomposition, ideal_point=self.ideal_point)
+                    off_FV = decompose(F[None, :], w, method=self.decomposition, ideal_point=self.ideal_point)
+
+                    # replace if better
+                    if np.all(off_FV < FV):
+                        pop.X[neighbor, :], pop.F[neighbor, :] = X, F
+

pymoo/algorithms/nsga3.py

@@ -10,7 +10,7 @@
 from pymoo.rand import random
 from pymoo.util.display import disp_multi_objective
 from pymoo.util.dominator import compare
-from pymoo.util.reference_directions import get_uniform_weights, get_ref_dirs_from_section, get_multi_layer_ref_dirs
+from pymoo.util.reference_directions import get_ref_dirs_from_section, get_multi_layer_ref_dirs
 
 
 class NSGA3(GeneticAlgorithm):

pymoo/cython/decomposition.pyx

@@ -0,0 +1,47 @@
+# distutils: language = c++
+# cython: boundscheck=False, wraparound=False, cdivision=True
+
+
+from my_math cimport c_norm
+from libcpp.vector cimport vector
+
+
+def cython_pbi(double[:,:] F, double[:] weights, double[:] ideal_point, double theta=0.01):
+    return c_pbi(F, weights, ideal_point, theta)
+
+
+
+cdef extern from "math.h":
+    double sqrt(double m)
+    double pow(double base, double exponent)
+
+
+cdef vector[double] c_pbi(double[:,:] F, double[:] weights, double[:] ideal_point, double theta=0.01):
+    cdef:
+        double d1, d2, f_max, norm
+        int i, j, n_dim
+        vector[double] pbi
+
+    n_points = F.shape[0]
+    n_obj = F.shape[1]
+
+    pbi = vector[double](n_points)
+    norm = c_norm(weights)
+
+    for i in range(n_points):
+
+        d1 = 0
+        for j in range(n_obj):
+            d1 += pow(((F[i,j] - ideal_point[j]) * weights[j]) / norm, 2.0)
+        d1 = sqrt(d1)
+
+        d2 = 0
+        for j in range(n_obj):
+            d2 += pow(F[i,j] - (ideal_point[j] - d1 * weights[j]), 2.0)
+        d2 = sqrt(d2)
+
+        pbi[i] = d1 + theta * d2
+
+    return pbi
+
+

pymoo/cython/my_math.pxd

@@ -0,0 +1,6 @@
+# distutils: language = c++
+# cython: boundscheck=False, wraparound=False, cdivision=True
+
+
+cdef double c_norm(double[:] v)
+cdef double[:,:] c_calc_perpendicular_distance(double[:,:] P, double[:,:] L)

pymoo/cython/my_math.pyx

@@ -0,0 +1,68 @@
+# distutils: language = c++
+# cython: boundscheck=False, wraparound=False, cdivision=True
+
+
+import numpy as np
+from libcpp.vector cimport vector
+
+
+
+def cython_calc_perpendicular_distance(double[:,:] P, double[:,:] L):
+    return np.array(c_calc_perpendicular_distance(P, L))
+
+
+
+cdef extern from "math.h":
+    double sqrt(double m)
+    double pow(double base, double exponent)
+
+
+cdef double c_norm(double[:] v):
+    cdef:
+        double val
+        int i
+    val = 0
+    for i in range(v.shape[0]):
+        val += pow(v[i], 2)
+    val = sqrt(val)
+    return val
+
+
+cdef double[:,:] c_calc_perpendicular_distance(double[:,:] P, double[:,:] L):
+    cdef :
+        int s_L, s_P, n_dim, i, j, k
+        double[:,:] M
+        vector[double] N
+        double norm, dot, perp_dist, norm_scalar_proj
+
+    s_L = L.shape[0]
+    s_P = P.shape[0]
+    n_dim = L.shape[1]
+
+    M = np.zeros((s_P, s_L), dtype=np.float64)
+
+    for i in range(s_L):
+
+        norm = c_norm(L[i, :])
+
+        N = vector[double](n_dim)
+        for k in range(n_dim):
+            N[k] = L[i, k] / norm
+
+        for j in range(s_P):
+
+            dot = 0
+            for k in range(n_dim):
+                dot += L[i, k] * P[j, k]
+            norm_scalar_proj = dot / norm
+
+            perp_dist = 0
+            for k in range(n_dim):
+                perp_dist += pow(norm_scalar_proj * N[k] - P[j, k], 2)
+            perp_dist = sqrt(perp_dist)
+
+            M[j, i] = perp_dist
+
+    return M
+
+