constraint_optimization_problem.py

import math
import re

import numpy as np


FEASIBILITY_EPS = 1e-4


CEC2017_NAMES = {
    1: "C01",
    2: "C02",
    3: "C03",
    4: "C04",
    5: "C05",
    6: "C06",
    7: "C07",
    8: "C08",
    9: "C09",
    10: "C10",
    11: "C11",
    12: "C12",
    13: "C13",
    14: "C14",
    15: "C15",
    16: "C16",
    17: "C17",
    18: "C18",
    19: "C19",
    20: "C20",
    21: "C21",
    22: "C22",
    23: "C23",
    24: "C24",
    25: "C25",
    26: "C26",
    27: "C27",
    28: "C28",
}


def _make_rng(rng=None, rng_seed=None):
    if isinstance(rng, np.random.Generator):
        return rng
    if isinstance(rng, np.random.RandomState):
        return np.random.default_rng(int(rng.randint(0, 2**32 - 1)))
    if rng_seed is None:
        rng_seed = rng if isinstance(rng, (int, np.integer)) else 42
    return np.random.default_rng(int(rng_seed))


def _as_array(x, dtype=np.float64):
    return np.asarray(x, dtype=dtype)


def _as_vector(x, dim, name, dtype=np.float64):
    x = _as_array(x, dtype=dtype).reshape(-1)
    if x.size == 1:
        return np.full(dim, x.item(), dtype=dtype)
    if x.size != dim:
        raise ValueError(f"{name} must have exactly {dim} elements, got {x.size}.")
    return x


def _ensure_2d(x, dtype=np.float64):
    x = _as_array(x, dtype=dtype)
    return x.reshape(1, -1) if x.ndim == 1 else x


def rotate_gen(dim, rng=None, rng_seed=None, dtype=np.float64):
    rng = _make_rng(rng=rng, rng_seed=rng_seed)
    a = rng.normal(loc=0.0, scale=1.0, size=(dim, dim)).astype(dtype)
    q, _ = np.linalg.qr(a)
    if np.linalg.det(q) < 0:
        q[:, 0] = -q[:, 0]
    return q


def _positive_violation(x, epsilon):
    x = np.array(x, copy=True)
    x[x <= 0.0] = 0.0
    return x


def _equal_violation(x, epsilon):
    x = np.abs(x)
    x[x <= epsilon] = 0.0
    return x


def _append_zero_column(x):
    return np.concatenate((x[:, 1:], np.zeros((x.shape[0], 1), dtype=x.dtype)), axis=1)


def _rosenbrock_chain(z):
    return np.sum(100.0 * (z[:, :-1] ** 2 - z[:, 1:]) ** 2 + (z[:, :-1] - 1.0) ** 2, axis=1)


def _discontinuous_transform(y):
    return np.where(np.abs(y) < 0.5, y, 0.5 * np.round(2.0 * y))


def _rotated_cumulative_sum(x):
    return np.cumsum(x, axis=1)


def _safe_nan(x):
    return np.nan_to_num(x, nan=1e9, posinf=1e9, neginf=-1e9)


class Basic_Problem:
    def __init__(self, lb, ub, name=None):
        self.lb = _as_vector(lb, self.dim, "lb")
        self.ub = _as_vector(ub, self.dim, "ub")
        self.name = name or self.__class__.__name__

    def eval(self, x):
        x = self.transform_solution(x)
        if x.ndim == 1:
            return self.func(x.reshape(1, -1))
        if x.ndim == 2:
            return self.func(x)
        return self.func(x.reshape(-1, x.shape[-1]))

    def transform_solution(self, x):
        n = np.ndim(x)
        x = _ensure_2d(x)
        x = np.clip(x, 0.0, 1.0)
        x = self.lb + x * (self.ub - self.lb)
        return x[0] if n == 1 else x

    def func(self, x):
        costs, constraints = self.calculate_result(x)
        return costs, constraints.T

    def get_optimal(self):
        return self.opt


class Constraint_Optimization_Problem(Basic_Problem):
    fun_id = None
    cons_num = None
    lb = -100.0
    ub = 100.0
    min_bound = lb
    max_bound = ub
    optimum = 0.0

    def __init_subclass__(cls, **kwargs):
        super().__init_subclass__(**kwargs)
        if "lb" in cls.__dict__:
            cls.min_bound = cls.lb
        elif "min_bound" in cls.__dict__:
            cls.lb = cls.min_bound

        if "ub" in cls.__dict__:
            cls.max_bound = cls.ub
        elif "max_bound" in cls.__dict__:
            cls.ub = cls.max_bound

    def __init__(
        self,
        dim,
        changeable=True,
        min_tensor=None,
        max_tensor=None,
        rng=None,
        rng_seed=None,
        dtype=np.float64,
    ):
        if self.fun_id not in CEC2017_NAMES:
            raise ValueError(f"{self.__class__.__name__}.fun_id must be in [1, 28].")

        self.dim = int(dim)
        self.changeable = changeable
        self.dtype = dtype
        self.rng = _make_rng(rng=rng, rng_seed=42 if rng_seed is None else rng_seed)
        self.aux_rng = _make_rng(rng_seed=47 if rng_seed is None else int(rng_seed) + 5)
        self.optimum = float(type(self).optimum)

        self.min_tensor = float(type(self).lb if min_tensor is None else min_tensor)
        self.max_tensor = float(type(self).ub if max_tensor is None else max_tensor)
        super().__init__(
            lb=np.full(self.dim, self.min_tensor, dtype=dtype),
            ub=np.full(self.dim, self.max_tensor, dtype=dtype),
            name=CEC2017_NAMES[self.fun_id],
        )

        self.rotated = rotate_gen(self.dim, self.rng, dtype=dtype).T
        self.rotated1 = rotate_gen(self.dim, self.aux_rng, dtype=dtype).T
        self.bias = self.rng.uniform(self.min_tensor, self.max_tensor, size=self.dim).astype(dtype)
        self.bias1 = self.aux_rng.uniform(self.min_tensor, self.max_tensor, size=self.dim).astype(dtype)
        self.opt = self.transform_solution(np.full(self.dim, 0.5, dtype=dtype))

    def __str__(self):
        return f"F{self.fun_id}: {self.name}"

    def z(self, x, rotate=False, bias=None):
        z = _ensure_2d(x, self.dtype) - (self.bias if bias is None else bias)
        return z @ self.rotated if rotate else z

    def zr1(self, x):
        return (_ensure_2d(x, self.dtype) - self.bias) @ self.rotated1

    def calculate_result(self, x, epsilon=FEASIBILITY_EPS):
        x = _ensure_2d(x, self.dtype)
        costs = _safe_nan(self.objective(x))
        constraints = _safe_nan(np.vstack(self.constraints(x, epsilon)))
        return costs, constraints

    def objective(self, x):
        raise NotImplementedError

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        raise NotImplementedError

    def cal_feasible_rate(self, penalties):
        penalties = _ensure_2d(penalties, self.dtype)
        return float(np.mean(np.all(penalties <= FEASIBILITY_EPS, axis=1)))


class C01(Constraint_Optimization_Problem):
    fun_id = 1
    cons_num = 1
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = _rotated_cumulative_sum(self.z(x))
        return np.sum(z**2, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.sum(z**2 - 5000.0 * np.cos(0.1 * math.pi * z) - 4000.0, axis=1)
        return [_positive_violation(c1, epsilon)]


class C02(C01):
    fun_id = 2
    lb = -100.0
    ub = 100.0

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        c1 = np.sum(z**2 - 5000.0 * np.cos(0.1 * math.pi * z) - 4000.0, axis=1)
        return [_positive_violation(c1, epsilon)]


class C03(Constraint_Optimization_Problem):
    fun_id = 3
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = _rotated_cumulative_sum(self.z(x))
        return np.sum(z**2, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.sum(z**2 - 5000.0 * np.cos(0.1 * math.pi * z) - 4000.0, axis=1)
        c2 = np.sum(-z * np.sin(0.1 * math.pi * z), axis=1)
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C04(Constraint_Optimization_Problem):
    fun_id = 4
    cons_num = 2
    lb = -10.0
    ub = 10.0

    def objective(self, x):
        z = self.z(x)
        return np.sum(z**2 - 10.0 * np.cos(2.0 * math.pi * z) + 10.0, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.sum(-z * np.sin(2.0 * z), axis=1)
        c2 = np.sum(z * np.sin(z), axis=1)
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon)]


class C05(Constraint_Optimization_Problem):
    fun_id = 5
    cons_num = 2
    lb = -10.0
    ub = 10.0

    def objective(self, x):
        z = self.z(x)
        return _rosenbrock_chain(z)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z1 = self.z(x, rotate=True)
        z2 = self.zr1(x)
        c1 = np.sum(z1**2 - 50.0 * np.cos(2.0 * math.pi * z1) - 40.0, axis=1)
        c2 = np.sum(z2**2 - 50.0 * np.cos(2.0 * math.pi * z2) - 40.0, axis=1)
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon)]


class C06(Constraint_Optimization_Problem):
    fun_id = 6
    cons_num = 6
    lb = -20.0
    ub = 20.0

    def objective(self, x):
        z = self.z(x)
        return np.sum(z**2 - 10.0 * np.cos(2.0 * math.pi * z) + 10.0, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        return [
            _equal_violation(np.sum(-z * np.sin(z), axis=1), epsilon),
            _equal_violation(np.sum(z * np.sin(math.pi * z), axis=1), epsilon),
            _equal_violation(np.sum(-z * np.cos(z), axis=1), epsilon),
            _equal_violation(np.sum(z * np.cos(math.pi * z), axis=1), epsilon),
            _equal_violation(np.sum(z * np.sin(2.0 * np.sqrt(np.abs(z))), axis=1), epsilon),
            _equal_violation(np.sum(-z * np.sin(2.0 * np.sqrt(np.abs(z))), axis=1), epsilon),
        ]


class C07(Constraint_Optimization_Problem):
    fun_id = 7
    cons_num = 2
    lb = -50.0
    ub = 50.0

    def objective(self, x):
        z = self.z(x)
        return np.sum(z * np.sin(z), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        value = np.sum(z - 100.0 * np.cos(0.5 * z) + 100.0, axis=1)
        return [_equal_violation(value, epsilon), _equal_violation(-value, epsilon)]


class C08(Constraint_Optimization_Problem):
    fun_id = 8
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        return np.max(self.z(x), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        even = _rotated_cumulative_sum(z[:, 0::2])
        odd = _rotated_cumulative_sum(z[:, 1::2])
        return [
            _equal_violation(np.sum(even**2, axis=1), epsilon),
            _equal_violation(np.sum(odd**2, axis=1), epsilon),
        ]


class C09(Constraint_Optimization_Problem):
    fun_id = 9
    cons_num = 2
    lb = -10.0
    ub = 10.0

    def objective(self, x):
        return np.max(self.z(x), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        odd = z[:, 1::2]
        even = z[:, 0::2]
        c1 = np.prod(odd, axis=1) if odd.shape[1] else np.zeros(z.shape[0], dtype=z.dtype)
        c2 = np.sum((even[:, :-1] ** 2 - even[:, 1:]) ** 2, axis=1)
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C10(Constraint_Optimization_Problem):
    fun_id = 10
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        return np.max(self.z(x), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.sum(_rotated_cumulative_sum(z) ** 2, axis=1)
        z_next = _append_zero_column(z)
        z0 = z.copy()
        z0[:, -1] = 0.0
        c2 = np.sum((z0 - z_next) ** 2, axis=1)
        return [_equal_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C11(Constraint_Optimization_Problem):
    fun_id = 11
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        return np.sum(self.z(x), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.prod(z, axis=1)
        z_next = _append_zero_column(z)
        z0 = z.copy()
        z0[:, -1] = 0.0
        c2 = np.sum((z0 - z_next) ** 2, axis=1)
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C12(Constraint_Optimization_Problem):
    fun_id = 12
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = self.z(x)
        return np.sum(z**2 - 10.0 * np.cos(2.0 * math.pi * z) + 10.0, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = 4.0 - np.sum(np.abs(z), axis=1)
        c2 = np.sum(z**2, axis=1) - 4.0
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon)]


class C13(Constraint_Optimization_Problem):
    fun_id = 13
    cons_num = 3
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = self.z(x)
        return _rosenbrock_chain(z)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.sum(z**2 - 10.0 * np.cos(2.0 * math.pi * z) + 10.0, axis=1) - 100.0
        c2 = np.sum(z, axis=1) - 2.0 * self.dim
        c3 = 5.0 - np.sum(z, axis=1)
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon), _positive_violation(c3, epsilon)]


class C14(Constraint_Optimization_Problem):
    fun_id = 14
    cons_num = 2
    lb = -100.0
    ub = 100.0
    optimum = 4.0

    def objective(self, x):
        z = self.z(x)
        return (
            -20.0 * np.exp(-0.2 * np.sqrt(np.sum(z**2, axis=1) / self.dim))
            + 20.0
            - np.exp(np.sum(np.cos(2.0 * math.pi * z), axis=1) / self.dim)
            + math.e
        )

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.sum(z[:, 1:] ** 2, axis=1) + 1.0 - np.abs(z[:, 0])
        c2 = np.sum(z**2, axis=1) - 4.0
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C15(Constraint_Optimization_Problem):
    fun_id = 15
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        return np.max(np.abs(self.z(x)), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        ob = self.objective(x)
        c1 = np.sum(z**2, axis=1) - 100.0 * self.dim
        c2 = np.cos(ob) + np.sin(ob)
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C16(Constraint_Optimization_Problem):
    fun_id = 16
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        return np.sum(np.abs(self.z(x)), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        ob = self.objective(x)
        c1 = np.sum(z**2, axis=1) - 100.0 * self.dim
        c2 = (np.cos(ob) + np.sin(ob)) ** 2 - np.exp(np.cos(ob) + np.sin(ob)) - 1.0 + math.e
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C17(Constraint_Optimization_Problem):
    fun_id = 17
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = self.z(x)
        a = np.sqrt(np.arange(1, self.dim + 1, dtype=self.dtype))
        return np.sum(z**2, axis=1) / 4000.0 + 1.0 - np.prod(np.cos(z / a), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = 1.0 - np.sum(np.sign(np.abs(z) - (np.sum(z**2, axis=1, keepdims=True) - z**2) - 1.0), axis=1)
        c2 = np.sum(z**2, axis=1) - 4.0 * self.dim
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C18(Constraint_Optimization_Problem):
    fun_id = 18
    cons_num = 3
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = _discontinuous_transform(self.z(x))
        return np.sum(z**2 - 10.0 * np.cos(2.0 * math.pi * z) + 10.0, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        y = self.z(x)
        y_next = np.concatenate((y[:, 1:], y[:, :1]), axis=1)
        c1 = 1.0 - np.sum(np.abs(y), axis=1)
        c2 = np.sum(y**2, axis=1) - 100.0 * self.dim
        c3 = np.sum(100.0 * (y**2 - y_next) ** 2, axis=1) + np.prod(np.sin((y - 1.0) * math.pi) ** 2, axis=1)
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon), _equal_violation(c3, epsilon)]


class C19(Constraint_Optimization_Problem):
    fun_id = 19
    cons_num = 2
    lb = -50.0
    ub = 50.0

    def objective(self, x):
        z = self.z(x)
        return np.sum(np.abs(z) ** 0.5 + 2.0 * np.sin(z**3), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        c1 = np.sum(-10.0 * np.exp(-0.2 * np.sqrt(z[:, :-1] ** 2 + z[:, 1:] ** 2)), axis=1) + (self.dim - 1.0) * 10.0 / math.exp(-5.0)
        c2 = np.sum(np.sin(2.0 * z) ** 2, axis=1) - 0.5 * self.dim
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon)]


class C20(Constraint_Optimization_Problem):
    fun_id = 20
    cons_num = 2
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = self.z(x)
        z_next = np.concatenate((z[:, 1:], z[:, :1]), axis=1)
        r2 = z**2 + z_next**2
        r = np.sqrt(r2)
        return np.sum(0.5 + (np.sin(r) ** 2 - 0.5) / (1.0 + 0.001 * r) ** 2, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x)
        s = np.sum(z, axis=1)
        c1 = np.cos(s) ** 2 - 0.25 * np.cos(s) - 0.125
        c2 = np.exp(np.cos(s)) - math.exp(0.25)
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon)]


class C21(C12):
    fun_id = 21
    lb = -100.0
    ub = 100.0

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        c1 = 4.0 - np.sum(np.abs(z), axis=1)
        c2 = np.sum(z**2, axis=1) - 4.0
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon)]


class C22(C13):
    fun_id = 22
    cons_num = 3
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = self.z(x, rotate=True)
        return _rosenbrock_chain(z)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        c1 = np.sum(z**2 - 10.0 * np.cos(2.0 * math.pi * z) + 10.0, axis=1) - 100.0
        c2 = np.sum(z, axis=1) - 2.0 * self.dim
        c3 = 5.0 - np.sum(z, axis=1)
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon), _positive_violation(c3, epsilon)]


class C23(C14):
    fun_id = 23
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = self.z(x, rotate=True)
        return (
            -20.0 * np.exp(-0.2 * np.sqrt(np.sum(z**2, axis=1) / self.dim))
            + 20.0
            - np.exp(np.sum(np.cos(2.0 * math.pi * z), axis=1) / self.dim)
            + math.e
        )

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        c1 = np.sum(z[:, 1:] ** 2, axis=1) + 1.0 - np.abs(z[:, 0])
        c2 = np.sum(z**2, axis=1) - 4.0
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C24(C15):
    fun_id = 24
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        return np.max(np.abs(self.z(x, rotate=True)), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        ob = self.objective(x)
        c1 = np.sum(z**2, axis=1) - 100.0 * self.dim
        c2 = np.cos(ob) + np.sin(ob)
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C25(C16):
    fun_id = 25
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        return np.sum(np.abs(self.z(x, rotate=True)), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        ob = self.objective(x)
        c1 = np.sum(z**2, axis=1) - 100.0 * self.dim
        c2 = (np.cos(ob) + np.sin(ob)) ** 2 - np.exp(np.cos(ob) + np.sin(ob)) - 1.0 + math.e
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C26(C17):
    fun_id = 26
    lb = -100.0
    ub = 100.0

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        c1 = 1.0 - np.sum(np.sign(np.abs(z) - (np.sum(z**2, axis=1, keepdims=True) - z**2) - 1.0), axis=1)
        c2 = np.sum(z**2, axis=1) - 4.0 * self.dim
        return [_positive_violation(c1, epsilon), _equal_violation(c2, epsilon)]


class C27(C18):
    fun_id = 27
    cons_num = 3
    lb = -100.0
    ub = 100.0

    def objective(self, x):
        z = _discontinuous_transform(self.z(x, rotate=True))
        return np.sum(z**2 - 10.0 * np.cos(2.0 * math.pi * z) + 10.0, axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        y = self.z(x, rotate=True)
        y_next = np.concatenate((y[:, 1:], y[:, :1]), axis=1)
        c1 = 1.0 - np.sum(np.abs(y), axis=1)
        c2 = np.sum(y**2, axis=1) - 100.0 * self.dim
        c3 = np.sum(100.0 * (y**2 - y_next) ** 2, axis=1) + np.prod(np.sin((y - 1.0) * math.pi) ** 2, axis=1)
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon), _equal_violation(c3, epsilon)]


class C28(C19):
    fun_id = 28
    lb = -50.0
    ub = 50.0

    def objective(self, x):
        z = self.z(x, rotate=True)
        return np.sum(np.abs(z) ** 0.5 + 2.0 * np.sin(z**3), axis=1)

    def constraints(self, x, epsilon=FEASIBILITY_EPS):
        z = self.z(x, rotate=True)
        c1 = np.sum(-10.0 * np.exp(-0.2 * np.sqrt(z[:, :-1] ** 2 + z[:, 1:] ** 2)), axis=1) + (self.dim - 1.0) * 10.0 / math.exp(-5.0)
        c2 = np.sum(np.sin(2.0 * z) ** 2, axis=1) - 0.5 * self.dim
        return [_positive_violation(c1, epsilon), _positive_violation(c2, epsilon)]


PROBLEM_CLASSES = {
    cls.fun_id: cls
    for cls in (
        C01,
        C02,
        C03,
        C04,
        C05,
        C06,
        C07,
        C08,
        C09,
        C10,
        C11,
        C12,
        C13,
        C14,
        C15,
        C16,
        C17,
        C18,
        C19,
        C20,
        C21,
        C22,
        C23,
        C24,
        C25,
        C26,
        C27,
        C28,
    )
}


def make_cec2017_problem(fun_id, *args, **kwargs):
    return PROBLEM_CLASSES[int(fun_id)](*args, **kwargs)


class CEC2017_Dataset:
    def __init__(self, data, batch_size=1, rng=None, rng_seed=None):
        self.data = list(data)
        self.batch_size = int(batch_size)
        self.N = len(self.data)
        self.ptr = list(range(0, self.N, self.batch_size))
        self.index = list(range(self.N))
        self.rng = _make_rng(rng=rng, rng_seed=rng_seed)

    @staticmethod
    def gen_exact_problem(ub, lb, dim, fun_id, rng=None, rng_seed=None, **kwargs):
        return make_cec2017_problem(fun_id, dim=dim, min_tensor=lb, max_tensor=ub, rng=rng, rng_seed=rng_seed, **kwargs)

    @staticmethod
    def gen_problem(ub, lb, dim, nums=100, rng=None, rng_seed=None, **kwargs):
        rng = _make_rng(rng=rng, rng_seed=rng_seed)
        return [
            make_cec2017_problem(int(rng.integers(1, 29)), dim=dim, min_tensor=lb, max_tensor=ub, rng=rng, **kwargs)
            for _ in range(nums)
        ]

    @staticmethod
    def get_datasets(
        suit,
        upperbound=None,
        train_batch_size=1,
        test_batch_size=1,
        difficulty=None,
        instance_seed=3849,
        user_train_list=None,
        user_test_list=None,
        **kwargs,
    ):
        match = re.fullmatch(r"cec2017[-_]?(\d+)", suit, flags=re.IGNORECASE)
        if match is None:
            raise ValueError(f"{suit} function suit is invalid or is not supported yet.")
        dim = int(match.group(1))
        rng = _make_rng(rng_seed=instance_seed)

        if difficulty not in ("easy", "difficult", "all", None):
            raise ValueError(f"{difficulty} difficulty is invalid.")
        if difficulty is None and user_test_list is None and user_train_list is None:
            difficulty = "all"

        easy_set = {1, 2, 3, 5, 12, 14, 17, 21}
        train_set, test_set, all_set = [], [], []
        for fun_id in range(1, 29):
            cls = PROBLEM_CLASSES[fun_id]
            ub = cls.ub if upperbound is None else upperbound
            lb = -ub if upperbound is not None else cls.lb
            problem = make_cec2017_problem(fun_id, dim=dim, min_tensor=lb, max_tensor=ub, rng=rng, **kwargs)
            if difficulty == "all":
                all_set.append(problem)
            elif user_train_list is None and user_test_list is None:
                (test_set if (difficulty == "easy") == (fun_id in easy_set) else train_set).append(problem)
            elif user_train_list is not None and user_test_list is not None:
                if fun_id in user_train_list:
                    train_set.append(problem)
                if fun_id in user_test_list:
                    test_set.append(problem)
            elif user_train_list is not None:
                (train_set if fun_id in user_train_list else test_set).append(problem)
            else:
                (test_set if fun_id in user_test_list else train_set).append(problem)

        if difficulty == "all":
            train_set = all_set.copy()
            test_set = all_set.copy()

        return CEC2017_Dataset(train_set, train_batch_size, rng=rng), CEC2017_Dataset(test_set, test_batch_size, rng=rng)

    def __getitem__(self, item):
        ptr = self.ptr[item]
        idx = self.index[ptr : min(ptr + self.batch_size, self.N)]
        return [self.data[i] for i in idx]

    def __len__(self):
        return len(self.ptr)

    def __add__(self, other):
        return CEC2017_Dataset(self.data + other.data, self.batch_size, rng=self.rng)

    def shuffle(self):
        self.index = self.rng.permutation(self.N).tolist()


def gen_exact_problem(ub, lb, dim, fun_id, rng=None, rng_seed=None, **kwargs):
    return CEC2017_Dataset.gen_exact_problem(ub, lb, dim, fun_id, rng=rng, rng_seed=rng_seed, **kwargs)