多目的最適化に特化したPythonライブラリpymooの使い方まとめ

$ pip install -U pymoo==0.6.0

import numpy as np
from pymoo.util.misc import stack
from pymoo.core.problem import Problem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.operators.sampling.lhs import LHS
from pymoo.operators.crossover.sbx import SBX
from pymoo.operators.mutation.pm import PolynomialMutation
from pymoo.termination import get_termination
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

# 問題設定（自作関数）
class MyProblem(Problem):
    def __init__(self):
        super().__init__(n_var=2,
                         n_obj=2,
                         n_constr=2,
                         xl=np.array([-2,-2]),
                         xu=np.array([2,2]),
                        )
    def _evaluate(self, X, out, *args, **kwargs):
        f1 = X[:,0]**2 + X[:,1]**2
        f2 = (X[:,0]-1)**2 + X[:,1]**2
        g1 = 2*(X[:, 0]-0.1) * (X[:, 0]-0.9) / 0.18
        g2 = - 20*(X[:, 0]-0.4) * (X[:, 0]-0.6) / 4.8
        out["F"] = np.column_stack([f1, f2])
        out["G"] = np.column_stack([g1, g2])

# パレート解の設定
# パレート解が既知の場合に使用（テスト用）
def func_pf(flatten=True, **kwargs):
    f1_a = np.linspace(0.1**2, 0.4**2, 100)
    f2_a = (np.sqrt(f1_a) - 1)**2
    f1_b = np.linspace(0.6**2, 0.9**2, 100)
    f2_b = (np.sqrt(f1_b) - 1)**2
    a, b = np.column_stack([f1_a, f2_a]), np.column_stack([f1_b, f2_b])
    return stack(a, b, flatten=flatten)
def func_ps(flatten=True, **kwargs):
    x1_a = np.linspace(0.1, 0.4, 50)
    x1_b = np.linspace(0.6, 0.9, 50)
    x2 = np.zeros(50)
    a, b = np.column_stack([x1_a, x2]), np.column_stack([x1_b, x2])
    return stack(a,b, flatten=flatten)
    
class MyTestProblem(MyProblem):
    def _calc_pareto_front(self, *args, **kwargs):
        return func_pf(**kwargs)
    def _calc_pareto_set(self, *args, **kwargs):
        return func_ps(**kwargs)

if __name__ == "__main__":
    # 問題の定義
    problem = MyTestProblem()
    
    # アルゴリズムの初期化（NSGA-IIを使用）
    algorithm = NSGA2(
        pop_size=40,
        n_offsprings=10,
        sampling=LHS(),
        crossover=SBX(prob=0.9, eta=15),
        mutation=PolynomialMutation(eta=20),
        eliminate_duplicates=True
    )
    
    # 終了条件（40世代）
    termination = get_termination("n_gen", 40)
    
    # 最適化の実行
    res = minimize(problem,
                   algorithm,
                   termination,
                   seed=1,
                   save_history=True,
                   verbose=True)
    
    # 結果の可視化
    ps = problem.pareto_set(use_cache=False, flatten=False)
    pf = problem.pareto_front(use_cache=False, flatten=False)
    
    # 入力変数空間
    # plot = Scatter(title = "Design Space", axis_labels="x")
    # plot.add(res.X, s=30, facecolors='none', edgecolors='r')
    # if ps is not None:
    #     plot.add(ps, plot_type="line", color="black", alpha=0.7)
    # plot.do()
    # plot.apply(lambda ax: ax.set_xlim(-0.5, 1.5))
    # plot.apply(lambda ax: ax.set_ylim(-2, 2))
    # plot.show()
    
    # 目的関数空間
    plot = Scatter(title = "Objective Space")
    plot.add(res.F)
    if pf is not None:
        plot.add(pf, plot_type="line", color="black", alpha=0.7)
    plot.show()

import numpy as np
from pymoo.core.problem import Problem

class MyProblem(Problem):
    def __init__(self):
        super().__init__(n_var=2, # 入力変数の次元
                         n_obj=2, # 目的関数の数
                         n_constr=2, # 制約条件の数
                         xl=np.array([-2,-2]), # 入力変数の下限
                         xu=np.array([2,2]), # 入力変数の上限
                        )
    def _evaluate(self, X, out, *args, **kwargs):
        # 目的関数
        f1 = X[:,0]**2 + X[:,1]**2
        f2 = (X[:,0]-1)**2 + X[:,1]**2
        # 制約条件
        g1 = 2*(X[:, 0]-0.1) * (X[:, 0]-0.9) / 0.18
        g2 = - 20*(X[:, 0]-0.4) * (X[:, 0]-0.6) / 4.8

        out["F"] = np.column_stack([f1, f2])
        out["G"] = np.column_stack([g1, g2])

problem = MyProblem()

import numpy as np
from pymoo.core.problem import ElementwiseProblem

class MyProblem(ElementwiseProblem):
    def __init__(self):
        super().__init__(n_var=2, # 入力変数の次元
                         n_obj=2, # 目的関数の数
                         n_constr=2, # 制約条件の数
                         xl=np.array([-2,-2]), # 入力変数の下限
                         xu=np.array([2,2]), # 入力変数の上限
                         )
    def _evaluate(self, x, out, *args, **kwargs):
        # 目的関数
        f1 = x[0]**2 + x[1]**2
        f2 = (x[0]-1)**2 + x[1]**2
        # 制約条件
        g1 = 2*(x[0]-0.1) * (x[0]-0.9) / 0.18
        g2 = - 20*(x[0]-0.4) * (x[0]-0.6) / 4.8

        out["F"] = [f1, f2]
        out["G"] = [g1, g2]

problem = MyProblem()

import numpy as np
from pymoo.problems.functional import FunctionalProblem
# 目的関数
objs = [
    lambda x: x[0]**2 + x[1]**2,
    lambda x: (x[0]-1)**2 + x[1]**2
]
# 制約条件
constr_ieq = [
    lambda x: 2*(x[0]-0.1) * (x[0]-0.9) / 0.18,
    lambda x: - 20*(x[0]-0.4) * (x[0]-0.6) / 4.8
]
functional_problem = FunctionalProblem(2,
                                       objs,
                                       constr_ieq=constr_ieq,
                                       xl=np.array([-2,-2]),
                                       xu=np.array([2,2]))

from pymoo.util.misc import stack

def func_pf(flatten=True, **kwargs):
        f1_a = np.linspace(0.1**2, 0.4**2, 100)
        f2_a = (np.sqrt(f1_a) - 1)**2

        f1_b = np.linspace(0.6**2, 0.9**2, 100)
        f2_b = (np.sqrt(f1_b) - 1)**2

        a, b = np.column_stack([f1_a, f2_a]), np.column_stack([f1_b, f2_b])
        return stack(a, b, flatten=flatten)

def func_ps(flatten=True, **kwargs):
        x1_a = np.linspace(0.1, 0.4, 50)
        x1_b = np.linspace(0.6, 0.9, 50)
        x2 = np.zeros(50)

        a, b = np.column_stack([x1_a, x2]), np.column_stack([x1_b, x2])
        return stack(a,b, flatten=flatten)

import numpy as np
from pymoo.util.misc import stack
from pymoo.core.problem import Problem

class MyTestProblem(MyProblem):
    def __init__(self):
        '''省略'''
    def _evaluate(self, X, out, *args, **kwargs):
        '''省略'''
    def _calc_pareto_front(self, *args, **kwargs):
        return func_pf(**kwargs)

    def _calc_pareto_set(self, *args, **kwargs):
        return func_ps(**kwargs)

problem = MyTestProblem()

from pymoo.problems.functional import FunctionalProblem


functional_test_problem = FunctionalProblem(2,
                                            objs,
                                            constr_ieq=constr_ieq,
                                            xl=-2,
                                            xu=2,
                                            func_pf=func_pf,
                                            func_ps=func_ps
                                            )

from pymoo.core.problem import ElementwiseProblem

class MyProblem(ElementwiseProblem):
    def __init__(self, **kwargs):
        super().__init__(n_var=10, n_obj=1, n_constr=0, xl=-5, xu=5, **kwargs)
    def _evaluate(self, x, out, *args, **kwargs):
         out["F"] = (x ** 2).sum()

from multiprocessing.pool import ThreadPool
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.optimize import minimize

n_threads = 8 # スレッド数
pool = ThreadPool(n_threads) # 初期化
problem = MyProblem(parallelization = ('starmap', pool.starmap)) # 問題設定
res = minimize(problem, GA()) # 計算の実行
print(f'Threads: {res.exec_time} sec')
pool.close()

import multiprocessing
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.optimize import minimize

n_process = 8 # プロセス数
pool = multiprocessing.Pool(n_process) # 初期化
problem = MyProblem(parallelization = ('starmap', pool.starmap)) # 問題設定
res = minimize(problem, GA()) # 計算の実行
print(f'Processes: {res.exec_time} sec')
pool.close()

from dask.distributed import Client
import numpy as np
from pymoo.core.problem import ElementwiseProblem
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.optimize import minimize

client = Client()

# 評価関数の定義
def fun(x):
    return {
        "F": np.sum(x ** 2)
    }

# 問題の定義（評価関数は使用しない）
class MyProblem(ElementwiseProblem):
    def __init__(self, *args, **kwargs):
        super().__init__(n_var=10, n_obj=1, n_constr=0, xl=-5, xu=5, *args, **kwargs)

# 問題を作成し，並列化をdaskに設定
problem = MyProblem(parallelization=("dask", client, fun))

# 実行
res = minimize(problem, GA())
print(f'Dask: {res.exec_time} sec')

client.close()

from multiprocessing.pool import ThreadPool
from pymoo.core.problem import Problem
import numpy as np
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.optimize import minimize

class MyProblem(Problem):
    def __init__(self, **kwargs):
        super().__init__(n_var=10, n_obj=1, n_constr=0, xl=-5, xu=5, **kwargs)
        self.pool = ThreadPool(8)

    def _evaluate(self, X, out, *args, **kwargs):
        # 評価関数の定義
        def my_eval(x):
            return (x ** 2).sum()
        # 変数の処理
        params = [[X[k]] for k in range(len(X))]
        # 関数値を並列処理で計算
        F = self.pool.starmap(my_eval, params)
        # 関数値を返す
        out["F"] = np.array(F)

problem = MyProblem()
res = minimize(MyProblem(), GA())
print(f'Threads: {res.exec_time} sec')
problem.pool.close()

import numpy as np
from dask.distributed import Client
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.core.problem import Problem
from pymoo.optimize import minimize


class MyProblem(Problem):

    def __init__(self, client, *args, **kwargs):
        super().__init__(n_var=10, n_obj=1, n_constr=0, xl=-5, xu=5, *args, **kwargs)
        self.client = client

    def _evaluate(self, X, out, *args, **kwargs):
        def fun(x):
            return np.sum(x ** 2)

        jobs = [self.client.submit(fun, x) for x in X]
        out["F"] = np.row_stack([job.result() for job in jobs])
        
client = Client(processes=False)
problem = MyProblem(client)
res = minimize(problem, GA())
print(f'Dask: {res.exec_time} sec')
client.close()

from pymoo.core.problem import ElementwiseProblem
import numpy as np
import matplotlib.pyplot as plt

class ConstrainedProblem(ElementwiseProblem):

    def __init__(self, **kwargs):
        super().__init__(n_var=2, n_obj=1, n_ieq_constr=1, n_eq_constr=0, xl=0, xu=2, **kwargs)

    def _evaluate(self, x, out, *args, **kwargs):
        out["F"] = x[0] ** 2 + x[1] ** 2
        out["G"] = 1.0 - (x[0] + x[1])

X1, X2 = np.meshgrid(np.linspace(0, 2, 500), np.linspace(0, 2, 500))

F = X1**2 + X2**2
plt.rc('font', family='serif')

levels = 5 * np.linspace(0, 1, 10)
plt.figure(figsize=(7, 5))
CS = plt.contour(X1, X2, F, levels, colors='black', alpha=0.5)
CS.collections[0].set_label("$f(x)$")

X = np.linspace(0, 1, 500)
plt.plot(X, 1-X, linewidth=2.0, color="green", linestyle='dotted', label="g(x)")

plt.scatter([0.5], [0.5], marker="*", color="red", s=200, label="Optimum (0.5, 0.5)")

plt.xlim(0, 2)
plt.ylim(0, 2)
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")

plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.12),
          ncol=4, fancybox=True, shadow=False)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from pymoo.core.problem import ElementwiseProblem

class ConstrainedProblemWithEquality(ElementwiseProblem):

    def __init__(self, **kwargs):
        super().__init__(n_var=2, n_obj=1, n_ieq_constr=1, n_eq_constr=1, xl=0, xu=1, **kwargs)

    def _evaluate(self, x, out, *args, **kwargs):
        out["F"] = x[0] + x[1]
        out["G"] = 1.0 - (x[0] + x[1])
        out["H"] = 3 * x[0] - x[1]

X1, X2 = np.meshgrid(np.linspace(0, 2, 500), np.linspace(0, 2, 500))

F = X1**2 + X2**2
plt.rc('font', family='serif')

levels = 5 * np.linspace(0, 1, 10)
plt.figure(figsize=(7, 5))
CS = plt.contour(X1, X2, F, levels, colors='black', alpha=0.5)
CS.collections[0].set_label("$f(x)$")

X = np.linspace(0, 1, 500)
plt.plot(X, 1-X, linewidth=2.0, color="green", linestyle='dotted', label="g(x)")
plt.plot(X, 3*X, linewidth=2.0, color="blue", linestyle='dotted', label="h(x)")

plt.scatter([0.25], [0.75], marker="*", color="red", s=200, label="Optimum (0.5, 0.5)")

plt.xlim(0, 2)
plt.ylim(0, 2)
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")

plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.12),
          ncol=4, fancybox=True, shadow=False)

plt.tight_layout()
plt.show()

from pymoo.algorithms.soo.nonconvex.de import DE
from pymoo.constraints.as_penalty import ConstraintsAsPenalty
from pymoo.optimize import minimize
from pymoo.core.evaluator import Evaluator
from pymoo.core.individual import Individual

problem = ConstrainedProblem()

algorithm = DE()

res = minimize(ConstraintsAsPenalty(problem, penalty=100.0),
               algorithm,
               seed=1,
               verbose=False)

res = Evaluator().eval(problem, Individual(X=res.X))

print("Best solution found: \nX = %s\nF = %s\nCV = %s" % (res.X, res.F, res.CV))

from pymoo.core.repair import Repair

class MyRepair(Repair):

    def _do(self, problem, X, **kwargs):
        X[:, 0] = 1/3 * X[:, 1]
        return X

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.optimize import minimize

algorithm = GA(repair=MyRepair())

res = minimize(ConstrainedProblemWithEquality(),
               algorithm,
               ('n_gen', 20),
               seed=1,
               verbose=True)

print("Best solution found: \nX = %s\nF = %s\nCV = %s" % (res.X, res.F, res.CV))

from pymoo.problems import get_problem
problem = get_problem("zdt1") # zdt1の場合

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.factory import get_problem
from pymoo.optimize import minimize

problem = get_problem("g1")
algorithm = GA(
    pop_size=100,
    eliminate_duplicates=True)
res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)
print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

from pymoo.core.problem import Problem
from pymoo.operators.sampling.rnd import FloatRandomSampling
from pymoo.util import plotting

problem = Problem(n_var=2, xl=0, xu=1)

sampling = FloatRandomSampling()

X = sampling(problem, 200).get("X")
plotting.plot(X, no_fill=True)

from pymoo.operators.sampling.lhs import LHS

sampling = LHS()

X = sampling(problem, 200).get("X")
plotting.plot(X, no_fill=True)

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.operators.selection.rnd import RandomSelection

selection = RandomSelection()
algorithm = GA(
    pop_size=100,
    selection=selection,
    eliminate_duplicates=True)

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.operators.selection.tournament import TournamentSelection
from pymoo.optimize import minimize
from pymoo.problems import get_problem


def binary_tournament(pop, P, _, **kwargs):
    # The P input defines the tournaments and competitors
    n_tournaments, n_competitors = P.shape

    if n_competitors != 2:
        raise Exception("Only pressure=2 allowed for binary tournament!")

    # the result this function returns
    import numpy as np
    S = np.full(n_tournaments, -1, dtype=np.int)

    # now do all the tournaments
    for i in range(n_tournaments):
        a, b = P[i]

        # if the first individual is better, choose it
        if pop[a].F < pop[b].F:
            S[i] = a

        # otherwise take the other individual
        else:
            S[i] = b

    return S


selection = TournamentSelection(pressure=2, func_comp=binary_tournament)

problem = get_problem("rastrigin")

algorithm = GA(pop_size=100, eliminate_duplicates=True)

res = minimize(problem, algorithm, termination=('n_gen', 100), verbose=False)

print(res.X)

import matplotlib.pyplot as plt
import numpy as np

from pymoo.core.individual import Individual
from pymoo.core.problem import Problem
from pymoo.operators.crossover.sbx import SBX


def show(eta_cross):
    problem = Problem(n_var=1, xl=0.0, xu=1.0)

    a, b = Individual(X=np.array([0.2])), Individual(X=np.array([0.8]))

    parents = [[a, b] for _ in range(5000)]

    off = SBX(prob=1.0, prob_var=1.0, eta=eta_cross).do(problem, parents)
    Xp = off.get("X")

    plt.hist(Xp, range=(0, 1), bins=200, density=True, color="red")
    plt.show()


show(1)

import matplotlib.pyplot as plt
import numpy as np

from pymoo.core.individual import Individual
from pymoo.core.problem import Problem
from pymoo.operators.crossover.pntx import PointCrossover, SinglePointCrossover, TwoPointCrossover

n_var, n_matings = 50, 30

problem = Problem(n_var=n_var, xl=0.0, xu=1.0, var_type=int)

a, b = Individual(X=np.arange(1, n_var + 1)), Individual(X=-np.arange(1, n_var + 1))

parents = [[a, b] for _ in range(n_matings)]

def show(M):
    plt.figure(figsize=(4, 6))
    plt.imshow(M, cmap='Greys', interpolation='nearest')
    plt.xlabel("Variables")
    plt.ylabel("Individuals")
    plt.show()


print("One Point Crossover")
off = SinglePointCrossover(prob=1.0).do(problem, parents)
Xp = off.get("X")
show(Xp[:n_matings] != a.X)

print("Two Point Crossover")
off = TwoPointCrossover(prob=1.0).do(problem, parents)
Xp = off.get("X")
show(Xp[:n_matings] != a.X)

print("K Point Crossover (k=4)")
off = PointCrossover(prob=1.0, n_points=4).do(problem, parents)
Xp = off.get("X")
show(Xp[:n_matings] != a.X)

from pymoo.operators.crossover.expx import ExponentialCrossover

off = ExponentialCrossover(prob=1.0, prob_exp=0.9).do(problem, parents)
Xp = off.get("X")
show((Xp[:n_matings] != a.X))

from pymoo.operators.crossover.ux import UniformCrossover
off = UniformCrossover(prob=1.0).do(problem, parents)
Xp = off.get("X")
show(Xp[:n_matings] != a.X)

from pymoo.operators.crossover.hux import HalfUniformCrossover

n_var, n_matings = 100, 100

problem = Problem(n_var=n_var, xl=0.0, xu=1.0, var_type=int)

a = Individual(X=np.full(n_var, False))
b = Individual(X=np.array([k % 5 == 0 for k in range(n_var)]))

parents = [[a, b] for _ in range(n_matings)]

off = HalfUniformCrossover(prob=1.0).do(problem, parents)
Xp = off.get("X")
show(Xp[:n_matings] != a.X)

diff_a_to_b = (a.X != b.X).sum()
diff_a_to_off = (a.X != Xp[:n_matings]).sum()

print("Difference in bits (a to b): ", diff_a_to_b)
print("Difference in bits (a to off): ", diff_a_to_off)
print("Crossover Rate: ", diff_a_to_off / diff_a_to_b)

import matplotlib.pyplot as plt
import numpy as np

from pymoo.core.population import Population
from pymoo.core.problem import Problem
from pymoo.operators.mutation.pm import PolynomialMutation


def show(eta_mut):
    problem = Problem(n_var=1, xl=0.0, xu=1.0)

    X = np.full((5000, 1), 0.5)
    pop = Population.new(X=X)

    mutation = PolynomialMutation(prob=1.0, eta=eta_mut)

    off = mutation(problem, pop)
    Xp = off.get("X")

    plt.hist(Xp, range=(0.0, 1.0), bins=200, density=True, color="red")
    plt.show()


show(30)

import matplotlib.pyplot as plt
import numpy as np

from pymoo.core.population import Population
from pymoo.core.problem import Problem
from pymoo.operators.mutation.bitflip import BitflipMutation

n_var, n_matings = 100, 50

problem = Problem(n_var=n_var, vtype=bool)

X = np.full((100, 100), False)
pop = Population.new(X=X)

mutation = BitflipMutation(prob=0.5, prob_var=0.3)
Xp = mutation(problem, pop).get("X")

plt.figure(figsize=(4, 4))
plt.imshow(X != Xp, cmap='Greys', interpolation='nearest')
plt.show()

from pymoo.algorithms.soo.nonconvex.de import DE
from pymoo.factory import get_problem
from pymoo.operators.sampling.lhs import LHS
from pymoo.optimize import minimize

problem = get_problem("ackley", n_var=10)
algorithm = DE(
    pop_size=100,
    sampling=LHS(),
    variant="DE/rand/1/bin",
    CR=0.5,
    F=0.3,
    dither="vector",
    jitter=False
)
res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)
print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

from pymoo.algorithms.soo.nonconvex.pso import PSO
from pymoo.factory import get_problem
from pymoo.operators.sampling.lhs import LHS
from pymoo.optimize import minimize

problem = get_problem("ackley", n_var=10)
algorithm = PSO(
    pop_size=25,
    sampling=LHS(),
    adaptive = True,
    w = 0.9,
    c1 = 2.0,
    c2 = 2.0,
    initial_velocity = "random",
    max_velocity_ratio = 0.20,
    pertube_best = True,
)
res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)
print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.factory import get_problem
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

problem = get_problem("zdt1")

algorithm = NSGA2(pop_size=100)

res = minimize(problem,
               algorithm,
               ('n_gen', 200),
               seed=1,
               verbose=False)

plot = Scatter()
plot.add(problem.pareto_front(), plot_type="line", color="black", alpha=0.7)
plot.add(res.F, color="red")
plot.show()

import numpy as np

from pymoo.algorithms.moo.rnsga2 import RNSGA2
from pymoo.factory import get_problem
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

problem = get_problem("zdt1", n_var=30)
pf = problem.pareto_front()

ref_points = np.array([[0.5, 0.2], [0.1, 0.6]])

algorithm = RNSGA2(
    ref_points=ref_points,
    pop_size=40,
    epsilon=0.01,
    normalization='front',
    extreme_points_as_reference_points=False,
    weights=np.array([0.5, 0.5]))

res = minimize(problem,
               algorithm,
               save_history=True,
               termination=('n_gen', 250),
               seed=1,
               pf=pf,
               disp=False)

Scatter().add(pf, label="pf").add(res.F, label="F").add(ref_points, label="ref_points").show()

from pymoo.algorithms.moo.nsga3 import NSGA3
from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

ref_dirs = get_reference_directions("das-dennis", 3, n_partitions=12)

algorithm = NSGA3(pop_size=92,
                  ref_dirs=ref_dirs)

res = minimize(get_problem("dtlz1"),
               algorithm,
               seed=1,
               termination=('n_gen', 600))

Scatter().add(res.F).show()

from pymoo.visualization.scatter import Scatter
from pymoo.util.ref_dirs import get_reference_directions
ref_dirs = get_reference_directions(
    "das-dennis", 
    3, # 次元
    n_partitions=12 # 分割数
)
Scatter().add(ref_dirs).show()

from pymoo.visualization.scatter import Scatter
from pymoo.util.ref_dirs import get_reference_directions

ref_dirs = get_reference_directions(
    "energy", 
    3, # 次元
    90, # 点の数
    seed=1, # 乱数シード
)
Scatter().add(ref_dirs).show()

from pymoo.visualization.scatter import Scatter
from pymoo.util.ref_dirs import get_reference_directions

ref_dirs = get_reference_directions(
    "multi-layer",
    get_reference_directions("das-dennis", 3, n_partitions=12, scaling=1.0),
    get_reference_directions("das-dennis", 3, n_partitions=12, scaling=0.5)
)
Scatter().add(ref_dirs).show()

from pymoo.visualization.pcp import PCP
from pymoo.util.ref_dirs import get_reference_directions

partitions = [3, 2, 1, 1]
ref_dirs = get_reference_directions("layer-energy", 8, partitions)
PCP().add(ref_dirs).show()

import numpy as np
from pymoo.algorithms.moo.unsga3 import UNSGA3
from pymoo.problems import get_problem
from pymoo.optimize import minimize

problem = get_problem("ackley", n_var=30)
ref_dirs = np.array([[1.0]])
algorithm = UNSGA3(ref_dirs, pop_size=100)
res = minimize(problem,
               algorithm,
               termination=('n_gen', 150),
               save_history=True,
               seed=1)

print("UNSGA3: Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

import numpy as np

from pymoo.algorithms.moo.rnsga3 import RNSGA3
from pymoo.factory import get_problem
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

problem = get_problem("zdt1")
pf = problem.pareto_front()

# Define reference points
ref_points = np.array([[0.3, 0.4], [0.8, 0.5]])

# Get Algorithm
algorithm = RNSGA3(
    ref_points=ref_points,
    pop_per_ref_point=50,
    mu=0.1)

res = minimize(problem,
               algorithm=algorithm,
               termination=('n_gen', 300),
               pf=pf,
               seed=1,
               verbose=False)

reference_directions = res.algorithm.survival.ref_dirs

plot = Scatter()
plot.add(pf, label="pf")
plot.add(res.F, label="F")
plot.add(ref_points, label="ref_points")
plot.add(reference_directions, label="ref_dirs")
plot.show()

import numpy as np
from pymoo.core.problem import ElementwiseProblem
from pymoo.core.duplicate import ElementwiseDuplicateElimination
from pymoo.algorithms.soo.nonconvex.brkga import BRKGA
from pymoo.optimize import minimize

class MyProblem(ElementwiseProblem):
    def __init__(self, my_list):
        self.correct = np.argsort(my_list)
        super().__init__(n_var=len(my_list), n_obj=1, n_ieq_constr=0, xl=0, xu=1)

    def _evaluate(self, x, out, *args, **kwargs):
        pheno = np.argsort(x)
        out["F"] = - float((self.correct == pheno).sum())
        out["pheno"] = pheno
        out["hash"] = hash(str(pheno))
        
class MyElementwiseDuplicateElimination(ElementwiseDuplicateElimination):
    def is_equal(self, a, b):
        return a.get("hash") == b.get("hash")
    
np.random.seed(2)
list_to_sort = np.random.random(20)
problem = MyProblem(list_to_sort)
print("Sorted by", np.argsort(list_to_sort))

algorithm = BRKGA(
    n_elites=200,
    n_offsprings=700,
    n_mutants=100,
    bias=0.7,
    eliminate_duplicates=MyElementwiseDuplicateElimination())

res = minimize(problem,
               algorithm,
               ("n_gen", 75),
               seed=1,
               verbose=False)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))
print("Solution", res.opt.get("pheno")[0])
print("Optimum ", np.argsort(list_to_sort))

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))
print("Solution", res.opt.get("pheno")[0])
print("Optimum ", np.argsort(list_to_sort))

from pymoo.algorithms.soo.nonconvex.nelder import NelderMead
from pymoo.problems import get_problem
from pymoo.optimize import minimize

problem = get_problem("sphere")

algorithm = NelderMead()

res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

from pymoo.algorithms.soo.nonconvex.pattern import PatternSearch
from pymoo.factory import get_problem
from pymoo.optimize import minimize

problem = get_problem("ackley", n_var=30)

algorithm = PatternSearch()

res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

from pymoo.algorithms.soo.nonconvex.cmaes import CMAES
from pymoo.problems import get_problem
from pymoo.optimize import minimize

problem = get_problem("sphere")

algorithm = CMAES()

res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)

print(f"Best solution found: \nX = {res.X}\nF = {res.F}\nCV= {res.CV}")

res = minimize(problem,
               algorithm,
               ('n_iter', 10),
               seed=1,
               verbose=True)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

res = minimize(problem,
               algorithm,
               ('n_evals', 50),
               seed=1,
               verbose=True)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

problem = get_problem("rastrigin")

algorithm = CMAES(restarts=10, restart_from_best=True)

res = minimize(problem,
               algorithm,
               ('n_evals', 2500),
               seed=1,
               verbose=False)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

import numpy as np
from pymoo.util.normalization import denormalize

# define an intitial point for the search
np.random.seed(1)
x0 = denormalize(np.random.random(problem.n_var), problem.xl, problem.xu)

algorithm = CMAES(x0=x0,
                 sigma=0.5,
                 restarts=2,
                 maxfevals=np.inf,
                 tolfun=1e-6,
                 tolx=1e-6,
                 restart_from_best=True,
                 bipop=True)

res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))

from pymoo.algorithms.moo.moead import MOEAD
from pymoo.optimize import minimize
from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions
from pymoo.visualization.scatter import Scatter

problem = get_problem("dtlz2")

ref_dirs = get_reference_directions("uniform", 3, n_partitions=12)

algorithm = MOEAD(
    ref_dirs,
    n_neighbors=15,
    prob_neighbor_mating=0.7,
)

res = minimize(problem,
               algorithm,
               ('n_gen', 200),
               seed=1,
               verbose=False)

Scatter().add(res.F).show()

from pymoo.algorithms.moo.ctaea import CTAEA
from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

problem = get_problem("c1dtlz1", None, 3, k=5)

ref_dirs = get_reference_directions("das-dennis", 3, n_partitions=12)

# create the algorithm object
algorithm = CTAEA(ref_dirs=ref_dirs)

# execute the optimization
res = minimize(problem,
               algorithm,
               ('n_gen', 600),
               seed=1,
               verbose=False
               )

sc = Scatter(legend=False, angle=(45, 30))
sc.add(problem.pareto_front(ref_dirs), plot_type='surface', alpha=0.2, label="PF", color="blue")
sc.add(res.F, facecolor="none", edgecolor="red")
sc.show()

from pymoo.termination.default import DefaultMultiObjectiveTermination

termination = DefaultMultiObjectiveTermination(
    xtol=1e-8,
    cvtol=1e-6,
    ftol=0.0025,
    period=30,
    n_max_gen=1000,
    n_max_evals=100000
)

from pymoo.termination.default import DefaultSingleObjectiveTermination

termination = DefaultSingleObjectiveTermination(
    xtol=1e-8,
    cvtol=1e-6,
    ftol=1e-6,
    period=20,
    n_max_gen=1000,
    n_max_evals=100000
)

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.problems import get_problem
from pymoo.termination import get_termination
from pymoo.optimize import minimize

problem = get_problem("zdt3")
algorithm = NSGA2(pop_size=100)
termination = get_termination("n_eval", 300)

res = minimize(problem,
               algorithm,
               termination,
               pf=problem.pareto_front(),
               seed=1,
               verbose=True)

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.problems import get_problem
from pymoo.optimize import minimize
from pymoo.termination.xtol import DesignSpaceTermination
from pymoo.termination.robust import RobustTermination

problem = get_problem("zdt3")
algorithm = NSGA2(pop_size=100)
termination = RobustTermination(DesignSpaceTermination(tol=0.01), period=20)

res = minimize(problem,
               algorithm,
               termination,
               pf=problem.pareto_front(),
               seed=1,
               verbose=False)

print(res.algorithm.n_gen)

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.problems import get_problem
from pymoo.optimize import minimize
from pymoo.termination.ftol import MultiObjectiveSpaceTermination
from pymoo.visualization.scatter import Scatter

problem = get_problem("zdt3")

algorithm = NSGA2(pop_size=100)

termination = RobustTermination(
    MultiObjectiveSpaceTermination(tol=0.005, n_skip=5), period=20)


res = minimize(problem,
               algorithm,
               termination,
               pf=True,
               seed=1,
               verbose=False)

print("Generations", res.algorithm.n_gen)
plot = Scatter(title="ZDT3")
plot.add(problem.pareto_front(use_cache=False, flatten=False), plot_type="line", color="black")
plot.add(res.F, facecolor="none", edgecolor="red", alpha=0.8, s=20)
plot.show()

from pymoo.optimize import minimize

res = minimize(problem,
               algorithm,
               termination,
               seed=1,
               save_history=True,
               verbose=True)

import numpy as np

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.problems import get_problem
from pymoo.optimize import minimize
from pymoo.util.display.column import Column
from pymoo.util.display.output import Output

class MyOutput(Output):
    def __init__(self):
        super().__init__()
        self.x_mean = Column("x_mean", width=13)
        self.x_std = Column("x_std", width=13)
        self.columns += [self.x_mean, self.x_std]

    def update(self, algorithm):
        super().update(algorithm)
        self.x_mean.set(np.mean(algorithm.pop.get("X")))
        self.x_std.set(np.std(algorithm.pop.get("X")))


problem = get_problem("zdt2")

algorithm = NSGA2(pop_size=100)

res = minimize(problem,
               algorithm,
               ('n_gen', 10),
               seed=1,
               output=MyOutput(),
               verbose=True)

import matplotlib.pyplot as plt
import numpy as np

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.problems import get_problem
from pymoo.core.callback import Callback
from pymoo.optimize import minimize


class MyCallback(Callback):

    def __init__(self) -> None:
        super().__init__()
        self.data["best"] = []

    def notify(self, algorithm):
        self.data["best"].append(algorithm.pop.get("F").min())


problem = get_problem("sphere")

algorithm = GA(pop_size=100)

res = minimize(problem,
               algorithm,
               ('n_gen', 10),
               seed=1,
               callback=MyCallback(),
               save_history=True,
               verbose=True)

val = res.algorithm.callback.data["best"]
plt.plot(np.arange(len(val)), val)
plt.show()

val = [e.pop.get("F").min() for e in res.history]
plt.plot(np.arange(len(val)), val)
plt.show()

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.problems import get_problem
from pymoo.optimize import minimize

problem = get_problem("sphere")
algorithm = GA(pop_size=5)
res = minimize(problem,
               algorithm,
               ('n_gen', 30),
               seed=1)

res.X, res.F, res.G, res.CV

res.algorithm, res.pop

res.pop.get("X"), res.pop.get("F")

import copy

# アルゴリズムをコピー
obj = copy.deepcopy(algorithm)

# 問題設定
obj.setup(problem, termination=termination, seed=1)

# 終了基準が満たされるまで実行
while obj.has_next():

    # アルゴリズムの反復
    obj.next()

    # 結果のプリント
    print(f"gen: {obj.n_gen} n_nds: {len(obj.opt)} constr: {obj.opt.get('CV').min()} ideal: {obj.opt.get('F').min(axis=0)}")

# 結果を取得
result = obj.result()

from pymoo.visualization.scatter import Scatter
from pymoo.problems import get_problem
F = get_problem("zdt3").pareto_front(use_cache=False, flatten=True)
plot = Scatter(legend=True)
plot.add(F, s=30, facecolors='none', edgecolors='r', label='red point')
plot.add(F, plot_type="line", color="black", linewidth=2, label='black line')
plot.show()

from pymoo.visualization.scatter import Scatter
from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions

ref_dirs = get_reference_directions("uniform", 3, n_partitions=12)
F = get_problem("dtlz1").pareto_front(ref_dirs)

plot = Scatter()
plot.add(F)
plot.show()

import numpy as np
from pymoo.visualization.scatter import Scatter

F = np.random.random((30, 4))

plot = Scatter(tight_layout=True)
plot.add(F, s=10)
plot.add(F[10], s=30, color="red")
plot.show()

from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions
from pymoo.visualization.pcp import PCP

# データセットの作成，DTLZ1のパレートフロントをスケーリングしたもの
ref_dirs = get_reference_directions("das-dennis", 6, n_partitions=5) * [2, 4, 8, 16, 32, 64]
F = get_problem("dtlz1").pareto_front(ref_dirs)

PCP().add(F).show()

from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions
from pymoo.visualization.pcp import PCP

# データセットの作成，DTLZ1のパレートフロントをスケーリングしたもの
ref_dirs = get_reference_directions("das-dennis", 6, n_partitions=5) * [2, 4, 8, 16, 32, 64]
F = get_problem("dtlz1").pareto_front(ref_dirs)

plot = PCP()
plot.set_axis_style(color="grey", alpha=0.5)
plot.add(F, color="grey", alpha=0.3)
plot.add(F[50], linewidth=5, color="red")
plot.add(F[75], linewidth=5, color="blue")
plot.show()

import numpy as np
from pymoo.visualization.heatmap import Heatmap

np.random.seed(1234)
F = np.random.random((4, 6))
Heatmap().add(F).show()

import numpy as np
from pymoo.visualization.heatmap import Heatmap

np.random.seed(1234)
F = np.random.random((4, 6))
plot = Heatmap(title=("Optimization", {'pad': 15}),
               cmap="Oranges_r",
               solution_labels=["Solution A", "Solution B", "Solution C", "Solution D"],
               labels=["profit", "cost", "sustainability", "environment", "satisfaction", "time"])
plot.add(F)
plot.show()

import numpy as np
from pymoo.visualization.petal import Petal

np.random.seed(1234)
F = np.random.random((1, 6))
print(F)
Petal(bounds=[0, 1]).add(F).show()

plot = Petal(bounds=[0, 1],
             cmap="tab20",
             labels=["profit", "cost", "sustainability", "environment", "satisfaction", "time"],
             title=("Solution A", {'pad': 20}))
plot.add(F)
plot.show()

F = np.random.random((6, 6))
plot = Petal(bounds=[0, 1], title=["Solution %s" % t for t in ["A", "B", "C", "D", "E", "F"]])
plot.add(F[:3])
plot.add(F[3:])
plot.show()

import numpy as np
from pymoo.visualization.radar import Radar

np.random.seed(3)

ideal_point = np.array([0.15, 0.1, 0.2, 0.1, 0.1])
nadir_point = np.array([0.85, 0.9, 0.95, 0.9, 0.85])

F = np.random.random((1, 5)) * (nadir_point - ideal_point) + ideal_point
print(F)

plot = Radar(bounds=[ideal_point, nadir_point], normalize_each_objective=False)
plot.add(F)
plot.show()

plot = Radar(bounds=[ideal_point, nadir_point], normalize_each_objective=True)
plot.add(F)
plot.show()

F = np.random.random((6, 5)) * (nadir_point - ideal_point) + ideal_point

plot = Radar(bounds=[ideal_point, nadir_point],
             axis_style={"color": 'blue'},
             point_style={"color": 'red', 's': 30})
plot.add(F[:3], color="red", alpha=0.8)
plot.add(F[3:], color="green", alpha=0.8)
plot.show()

from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions
from pymoo.visualization.radviz import Radviz

ref_dirs = get_reference_directions("uniform", 6, n_partitions=5)
F = get_problem("dtlz1").pareto_front(ref_dirs)
Radviz().add(F).show()

plot = Radviz(title="Optimization",
              legend=(True, {'loc': "upper left", 'bbox_to_anchor': (-0.1, 1.08, 0, 0)}),
              labels=["profit", "cost", "sustainability", "environment", "satisfaction", "time"],
              endpoint_style={"s": 70, "color": "green"})
plot.set_axis_style(color="black", alpha=1.0)
plot.add(F, color="grey", s=20)
plot.add(F[65], color="red", s=70, label="Solution A")
plot.add(F[72], color="blue", s=70, label="Solution B")
plot.show()

from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions
from pymoo.visualization.radviz import Radviz

ref_dirs = get_reference_directions("uniform", 6, n_partitions=5)
F = get_problem("dtlz1").pareto_front(ref_dirs)
Radviz().add(F).show()

plot = Radviz(title="Optimization",
              legend=(True, {'loc': "upper left", 'bbox_to_anchor': (-0.1, 1.08, 0, 0)}),
              labels=["profit", "cost", "sustainability", "environment", "satisfaction", "time"],
              endpoint_style={"s": 70, "color": "green"})
plot.set_axis_style(color="black", alpha=1.0)
plot.add(F, color="grey", s=20)
plot.add(F[65], color="red", s=70, label="Solution A")
plot.add(F[72], color="blue", s=70, label="Solution B")
plot.show()

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.problems import get_problem
from pymoo.optimize import minimize
import numpy as np
import matplotlib.pyplot as plt

problem = get_problem("ackley")
algorithm = GA()

res = minimize(problem,
               algorithm,
               termination=('n_gen', 10),
               seed=1,
               save_history=True)

n_evals = np.array([e.evaluator.n_eval for e in res.history])
opt = np.array([e.opt[0].F for e in res.history])

plt.title("Convergence")
plt.plot(n_evals, opt, "--")
plt.xlabel("Number of evaluations")
plt.ylabel("Fitness")
plt.show()

res.history[5].opt[0].X

import numpy as np
from pymoo.problems import get_problem
from pymoo.visualization.scatter import Scatter

# zdt1のパレートフロント
pf = get_problem("zdt1").pareto_front()

# アルゴリズムにより得られた結果
A = pf[::10] * 1.1

# 結果の描画
Scatter(legend=True).add(pf, label="Pareto-front").add(A, label="Result").show()

from pymoo.indicators.gd import GD

ind = GD(pf)
print("GD", ind(A))

from pymoo.indicators.gd_plus import GDPlus

ind = GDPlus(pf)
print("GD+", ind(A))

from pymoo.indicators.igd import IGD

ind = IGD(pf)
print("IGD", ind(A))

from pymoo.indicators.igd_plus import IGDPlus

ind = IGDPlus(pf)
print("IGD+", ind(A))

from pymoo.indicators.hv import HV

ref_point = np.array([1.2, 1.2])

ind = HV(ref_point=ref_point)
print("HV", ind(A))

import numpy as np
from pymoo.problems import get_problem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.optimize import minimize
from pymoo.util.running_metric import RunningMetricAnimation

problem = get_problem("zdt1")
algorithm = NSGA2(pop_size=50)

res = minimize(problem,
               algorithm,
               ('n_gen', 40),
               seed=1,
               verbose=False,
              save_history=True)

running = RunningMetricAnimation(delta_gen=10,
                        n_plots=4,
                        key_press=False,
                        do_show=True)

for algorithm in res.history:
    running.update(algorithm)