class BPP(Problem):
    def __init__(
        self,
        items: Iterable[int],
        capacity: int,
        seed: int = 42,
        *args,
        **kwargs,
    ):
        self._items = tuple(items)
        self._capacity = capacity
        dim = len(self._items)
        assert len(self._items) > 0
        assert self._capacity > 0

        bounds = list((0, dim - 1) for _ in range(dim))
        super().__init__(dimension=dim, bounds=bounds, name="BPP", seed=seed)

    def evaluate(self, individual: Sequence | Solution) -> tuple[float]:
        """Evaluates the candidate individual with the information of the Bin Packing.
        The fitness of the solution is the amount of unused space, as well as the
        number of bins for a specific solution. Falkenauer (1998) performance metric
        defined as:
            (x) = \\frac{\\sum_{k=1}^{N} \\left(\\frac{fill_k}{C}\\right)^2}{N}

        Args:
            individual (Sequence | Solution): Individual to evaluate

        Returns:
            Tuple[float]: Falkenauer Fitness
        """
        if len(individual) != self._dimension:
            msg = f"Mismatch between individual variables ({len(individual)}) and instance variables ({self._dimension}) in {self.__class__.__name__}"
            raise ValueError(msg)

        used_bins = np.max(individual).astype(int) + 1
        fill_i = np.zeros(used_bins)

        for item_idx, bin in enumerate(individual):
            fill_i[bin] += self._items[item_idx]

        fitness = (
            sum(((f_i / self._capacity) * (f_i / self._capacity)) for f_i in fill_i)
            / used_bins
        )
        if isinstance(individual, Solution):
            individual.fitness = fitness
            individual.objectives = (fitness,)

        return (fitness,)

    def __call__(self, individual: Sequence | Solution) -> tuple[float]:
        return self.evaluate(individual)

    def __repr__(self):
        return f"BPP<n={self._dimension},C={self._capacity},I={self._items}>"

    def __len__(self):
        return self._dimension

    def __array__(self, dtype=np.int32, copy: Optional[bool] = False) -> npt.ArrayLike:
        return np.asarray([self._capacity, *self._items], dtype=dtype, copy=copy)

    def create_solution(self) -> Solution:
        items = list(range(self._dimension))
        return Solution(variables=items)

    def to_file(self, filename: str = "instance.bpp"):
        with open(filename, "w") as file:
            file.write(f"{len(self)}\t{self._capacity}\n\n")
            content = "\n".join(str(i) for i in self._items)
            file.write(content)

    @classmethod
    def from_file(cls, filename: str):
        with open(filename) as f:
            lines = f.readlines()
            lines = [line.rstrip() for line in lines]

        (_, capacity) = lines[0].split()
        items = list(int(i) for i in lines[2:])

        return cls(items=items, capacity=int(capacity))

    def to_instance(self) -> Instance:
        _vars = [self._capacity, *self._items]
        return Instance(variables=_vars)

class BPPDomain(Domain):
    __capacity_approaches = ("evolved", "percentage", "fixed")
    __feat_names = names = "mean,std,median,max,min,tiny,small,medium,large,huge".split(
        ","
    )

    def __init__(
        self,
        dimension: int = 50,
        min_i: int = 1,
        max_i: int = 1000,
        capacity_approach: str = "fixed",
        max_capacity: int = 100,
        capacity_ratio: float = 0.8,
        seed: int = 42,
    ):
        if dimension < 0:
            raise ValueError(f"Expected dimension > 0 got {dimension}")
        if min_i < 0:
            raise ValueError(f"Expected min_i > 0 got {min_i}")
        if max_i < 0:
            raise ValueError(f"Expected max_i > 0 got {max_i}")
        if min_i > max_i:
            raise ValueError(
                f"Expected min_i to be less than max_i got ({min_i}, {max_i})"
            )

        self._dimension = dimension
        self._min_i = min_i
        self._max_i = max_i
        self._max_capacity = max_capacity

        if capacity_ratio < 0.0 or capacity_ratio > 1.0 or not float(capacity_ratio):
            self.capacity_ratio = 0.8  # Default
            msg = "The capacity ratio must be a float number in the range [0.0-1.0]. Set as 0.8 as default."
            print(msg)
        else:
            self.capacity_ratio = capacity_ratio

        if capacity_approach not in self.__capacity_approaches:
            msg = f"The capacity approach {capacity_approach} is not available. Please choose between {self.__capacity_approaches}. Evolved approach set as default."
            print(msg)
            self._capacity_approach = "fixed"
        else:
            self._capacity_approach = capacity_approach

        bounds = [(1.0, self._max_capacity)] + [
            (self._min_i, self._max_i) for _ in range(self._dimension)
        ]
        super().__init__(dimension=dimension, bounds=bounds, name="BPP", seed=seed)

    @property
    def capacity_approach(self):
        return self._capacity_approach

    @capacity_approach.setter
    def capacity_approach(self, app):
        """Setter for the Maximum capacity generator approach.
        It forces to update the variable to one of the specify values

        Args:
            app (str): Approach for setting the capacity. It should be fixed, evolved or percentage.
        """
        if app not in self.__capacity_approaches:
            msg = f"The capacity approach {app} is not available. Please choose between {self.__capacity_approaches}. Evolved approach set as default."
            print(msg)
            self._capacity_approach = "fixed"
        else:
            self._capacity_approach = app

    def generate_instances(self, n: int = 1) -> List[Instance]:
        """Generates N instances for the domain.

        Args:
            n (int, optional): Number of instances to generate. Defaults to 1.

        Returns:
            List[Instance]: A list of Instance objects created from the raw numpy generation
        """
        instances = np.empty(shape=(n, self.dimension + 1), dtype=np.int32)
        instances = self._rng.integers(
            low=self._min_i, high=self._max_i, size=(n, self._dimension + 1), dtype=int
        )
        # Sets the capacity according to the method
        match self.capacity_approach:
            case "evolved":
                instances[:, 0] = self._rng.integers(1, self._max_capacity, size=n)
            case "percentage":
                instances[:, 0] = (
                    np.sum(instances[:, 1:], axis=1, dtype=int) * self.capacity_ratio
                )
            case "fixed":
                instances[:, 0] = self._max_capacity
        return list(Instance(i) for i in instances)

    def extract_features(self, instances: Sequence[Instance]) -> np.ndarray:
        """Extract the features of the instance based on the BPP domain.
           For the BPP the features are:
           N, Capacity, MeanWeights, MedianWeights, VarianceWeights, MaxWeight,
           MinWeight, Huge, Large, Medium, Small, Tiny

        Args:
            instances (Instance): Instances to extract the features from

        Returns:
            np.ndarray: Values of each feature
        """
        if not isinstance(instances, np.ndarray):
            instances = np.asarray(instances)

        norm_variables = np.asarray(instances, copy=True)
        norm_variables[:, 1:] = norm_variables[:, 1:] / norm_variables[:, [0]]

        return np.column_stack(
            [
                np.mean(norm_variables, axis=1),
                np.std(norm_variables, axis=1),
                np.median(norm_variables, axis=1),
                np.max(norm_variables, axis=1),
                np.min(norm_variables, axis=1),
                np.mean(norm_variables > 0.5, axis=1),  # Huge
                np.mean(
                    (0.5 >= norm_variables) & (norm_variables > 0.33333333333), axis=1
                ),
                np.mean(
                    (0.33333333333 >= norm_variables) & (norm_variables > 0.25), axis=1
                ),
                np.mean(0.25 >= norm_variables, axis=1),  # Small
                np.mean(0.1 >= norm_variables, axis=1),  # Tiny
            ],
        ).astype(np.float32)

    def extract_features_as_dict(
        self, instances: Sequence[Instance]
    ) -> List[Dict[str, np.float32]]:
        """Creates a dictionary with the features of the instances.
        The key are the names of each feature and the values are
        the values extracted from instance.

        Args:
            instances (Sequence[Instance]): Instances to extract the features from.
        Returns:
            Dict[str, np.float32]: Dictionary with the names/values of each feature
        """
        features = self.extract_features(instances)
        named_features: list[dict[str, np.float32]] = [{}] * len(features)
        for i, feats in enumerate(features):
            named_features[i] = {k: v for k, v in zip(BPPDomain.__feat_names, feats)}

        return named_features

    def generate_problems_from_instances(
        self, instances: Sequence[Instance]
    ) -> List[Problem]:
        if not isinstance(instances, np.ndarray):
            instances = np.asarray(instances)

        # Assume evolved capacities
        capacities = instances[:, 0].astype(np.int32)
        match self.capacity_approach:
            case "percentage":
                capacities[:] = (
                    np.sum(instances[:, 1:], axis=1) * self.capacity_ratio
                ).astype(np.int32)
            case "fixed":
                capacities[:] = self._max_capacity
        return list(
            BPP(items=instances[i, 1:], capacity=capacities[i])
            for i in range(len(instances))
        )

class Knapsack(Problem):
    def __init__(
        self,
        profits: Sequence[int],
        weights: Sequence[int],
        capacity: int = 0,
        seed: int = 42,
        *args,
        **kwargs,
    ):
        if len(profits) != len(weights):
            raise ValueError(
                f"The number of weights and profits is different in Knapsack. Got {len(weights)} weights and {len(profits)} profits"
            )
        if capacity <= 0:
            raise ValueError(f"Capacity must be a positive integer. Got {capacity}")

        super().__init__(dimension=len(profits), bounds=[], name="KP", seed=seed)

        self.weights = weights
        self.profits = profits
        self.capacity = capacity
        self.penalty_factor = 100.0

    def get_bounds_at(self, i: int) -> tuple:
        if i < 0 or i > self._dimension:
            raise ValueError(
                f"Index {i} out-of-range. The bounds are 0-{self._dimension} "
            )
        return (0, 1)

    @property
    def bounds(self):
        return list((0, 1) for _ in range(self._dimension))

    def evaluate(self, individual: Sequence | Solution | np.ndarray) -> Tuple[float]:
        """Evaluates the candidate individual with the information of the Knapsack

        Args:
            individual (Sequence | Solution): Individual to evaluate

        Raises:
            ValueError: Raises an error if the len(individual) != len(profits or weights)

        Returns:
            Tuple[float]: Profit
        """

        if len(individual) != self._dimension:
            msg = f"Mismatch between individual variables and instance variables in {self.__class__.__name__}"
            raise ValueError(msg)

        profit = np.dot(individual, self.profits)
        packed = np.dot(individual, self.weights)
        difference = max(0, packed - self.capacity)
        penalty = self.penalty_factor * difference
        profit -= penalty

        return (profit,)

    def __call__(self, individual: Sequence | Solution | np.ndarray) -> Tuple[float]:
        return self.evaluate(individual)

    def __array__(self, dtype=np.int32, copy: Optional[bool] = None) -> npt.ArrayLike:
        """Creates a numpy array from the Knapsack instance description.

        Returns:
            npt.ArrayLike: 1d numpy array of size 1 + (2 * dimension)
        """
        return np.asarray(
            [
                self.capacity,
                *list(
                    itertools.chain.from_iterable([*zip(self.weights, self.profits)])
                ),
            ],
            dtype=dtype,
            copy=copy,
        )

    def __repr__(self):
        return f"KP<n={len(self.profits)},C={self.capacity}>"

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

    def create_solution(self) -> Solution:
        chromosome = self._rng.integers(low=0, high=1, size=self._dimension)
        return Solution(variables=chromosome)

    def to_file(self, filename: str = "instance.kp"):
        with open(filename, "w") as file:
            file.write(f"{len(self)}\t{self.capacity}\n\n")
            content = "\n".join(
                f"{w_i}\t{p_i}" for w_i, p_i in zip(self.weights, self.profits)
            )
            file.write(content)

    @classmethod
    def from_file(cls, filename: str):
        content = np.loadtxt(filename, dtype=int)
        capacity = content[0][1]
        weights, profits = content[1:, 0], content[1:, 1]
        return cls(profits=profits, weights=weights, capacity=capacity)

    def to_instance(self) -> Instance:
        _vars = [self.capacity] + list(
            itertools.chain.from_iterable([*zip(self.weights, self.profits)])
        )
        return Instance(variables=_vars)

class KnapsackDomain(Domain):
    __capacity_approaches = ("evolved", "percentage", "fixed")

    def __init__(
        self,
        dimension: int = 50,
        min_p: int = 1,
        min_w: int = 1,
        max_p: int = 1_000,
        max_w: int = 1_000,
        capacity_approach: str = "evolved",
        max_capacity: int = int(1e5),
        capacity_ratio: float = 0.8,
        seed: Optional[int] = None,
    ):
        self.min_p = min_p
        self.min_w = min_w
        self.max_p = max_p
        self.max_w = max_w
        self.max_capacity = max_capacity

        if capacity_ratio < 0.0 or capacity_ratio > 1.0 or not float(capacity_ratio):
            self.capacity_ratio = 0.8  # Default
            msg = "The capacity ratio must be a float number in the range [0.0-1.0]. Set as 0.8 as default."
            print(msg)
        else:
            self.capacity_ratio = capacity_ratio

        if capacity_approach not in self.__capacity_approaches:
            msg = f"The capacity approach {capacity_approach} is not available. Please choose between {self.__capacity_approaches}. Evolved approach set as default."
            print(msg)
            self._capacity_approach = "evolved"
        else:
            self._capacity_approach = capacity_approach

        bounds = [(1.0, self.max_capacity)] + [
            (min_w, max_w) if i % 2 == 0 else (min_p, max_p)
            for i in range(2 * dimension)
        ]
        super().__init__(
            dimension=dimension,
            bounds=bounds,
            name="KP",
            feat_names="capacity,max_p,max_w,min_p,min_w,avg_eff,mean,std".split(","),
            seed=seed,
        )

    @property
    def capacity_approach(self):
        return self._capacity_approach

    @capacity_approach.setter
    def capacity_approach(self, app):
        """Setter for the Maximum capacity generator approach.
        It forces to update the variable to one of the specify values

        Args:
            app (str): Approach for setting the capacity. It should be fixed, evolved or percentage.
        """
        if app not in self.__capacity_approaches:
            msg = f"The capacity approach {app} is not available. Please choose between {self.__capacity_approaches}. Evolved approach set as default."
            print(msg)
            self._capacity_approach = "evolved"
        else:
            self._capacity_approach = app

    def generate_instances(self, n: int = 1) -> List[Instance]:
        """Generates N instances for the domain.

        Args:
            n (int, optional): Number of instances to generate. Defaults to 1.

        Returns:
            List[Instance]: A list of Instance objects created from the raw numpy generation
        """
        weights_and_profits = np.empty(shape=(n, self.dimension * 2), dtype=np.int32)
        weights_and_profits[:, 0::2] = self._rng.integers(
            low=self.min_w, high=self.max_w, size=(n, self.dimension)
        )
        weights_and_profits[:, 1::2] = self._rng.integers(
            low=self.min_p, high=self.max_p, size=(n, self.dimension)
        )
        # Assume fixed
        capacities = np.full(n, fill_value=self.max_capacity, dtype=np.int32)
        match self.capacity_approach:
            case "evolved":
                capacities[:] = self._rng.integers(1, self.max_capacity, size=n)
            case "percentage":
                capacities[:] = (
                    np.sum(weights_and_profits[:, 1::2], axis=1) * self.capacity_ratio
                ).astype(np.int32)
        return list(
            Instance(i) for i in np.column_stack((capacities, weights_and_profits))
        )

    def extract_features(
        self, instances: Sequence[Instance] | np.ndarray
    ) -> np.ndarray:
        """Extract the features of the instance based on the domain

        Args:
            instances (Sequence[Instance]): Instances to extract the features from.

        Returns:
            ArrayLike: 2d array with the features of each instance
        """

        if not isinstance(instances, np.ndarray):
            instances = np.asarray(instances, copy=True)

        features = np.empty(shape=(len(instances), 8), dtype=np.float32)
        weights = instances[:, 1::2]
        profits = instances[:, 2::2]
        features[:, 0] = instances[:, 0]  # Qs
        features[:, 1] = np.max(profits, axis=1)
        features[:, 2] = np.max(weights, axis=1)
        features[:, 3] = np.min(profits, axis=1)
        features[:, 4] = np.min(weights, axis=1)
        features[:, 5] = np.mean(profits / weights)
        features[:, 6] = np.mean(instances[:, 1:], axis=1)
        features[:, 7] = np.std(instances[:, 1:], axis=1)

        return features

    def extract_features_as_dict(
        self, instances: Sequence[Instance] | np.ndarray
    ) -> List[Dict[str, np.float32]]:
        """Creates a dictionary with the features of the instance.
        The key are the names of each feature and the values are
        the values extracted from instance.

        Args:
            instances (Sequence[Instance]): Instances to extract the features from. They should in the an array form.

        Returns:
            Dict[str, float]: Dictionary with the names/values of each feature
        """
        features = self.extract_features(instances)
        named_features: list[dict[str, np.float32]] = [{}] * len(features)
        for i, feats in enumerate(features):
            named_features[i] = {k: v for k, v in zip(self.feat_names, feats)}

        return named_features

    def generate_problems_from_instances(
        self, instances: Sequence[Instance] | np.ndarray
    ) -> List:
        """Generates a List of Knapsack objects from the instances

        Args:
            instances (Sequence[Instance]): Instances to create the problems from

        Returns:
            List: List containing len(instances) objects of type Knapsack
        """
        if not isinstance(instances, np.ndarray):
            instances = np.asarray(instances)

        capacities = instances[:, 0].astype(int)
        weights = instances[:, 1::2].astype(int)
        profits = instances[:, 2::2].astype(int)
        # Sets the capacity according to the method
        if self.capacity_approach == "percentage":
            capacities[:] = (np.sum(weights, axis=1) * self.capacity_ratio).astype(
                np.int32
            )
        elif self.capacity_approach == "fixed":
            capacities[:] = self.max_capacity
        return list(
            Knapsack(profits=profits[i], weights=weights[i], capacity=capacities[i])
            for i in range(len(instances))
        )

class TSP(Problem):
    """Symmetric Travelling Salesman Problem"""

    def __init__(
        self,
        nodes: int,
        coords: np.ndarray,
        seed: int = 42,
        *args,
        **kwargs,
    ):
        """Creates a new Symmetric Travelling Salesman Problem

        Args:
            nodes (int): Number of nodes/cities in the instance to solve
            coords (np.ndarray(N, 2)): Coordinates of each node/city.
        """
        self._nodes = nodes
        if coords.shape[1] != 2:
            raise ValueError(
                f"Expected coordinates shape to be (N, 2). Instead coords has the following shape: {coords.shape}"
            )
        if not isinstance(coords, np.ndarray):
            coords = np.asarray(coords)

        self._coords = coords
        x_min, y_min = np.min(self._coords, axis=0)
        x_max, y_max = np.max(self._coords, axis=0)
        bounds = list(((x_min, y_min), (x_max, y_max)) for _ in range(self._nodes))
        super().__init__(dimension=nodes, bounds=bounds, name="TSP", seed=seed)

        self._distances = np.zeros((self._nodes, self._nodes))
        differences = self._coords[:, np.newaxis, :] - self._coords[np.newaxis, :, :]
        self._distances = np.sqrt(np.sum(differences**2, axis=-1))

    def __evaluate_constraints(self, individual: Sequence | Solution) -> bool:
        counter = Counter(individual)
        if any(counter[c] != 1 for c in counter if c != 0) or (
            individual[0] != 0 or individual[-1] != 0
        ):
            return False
        return True

    def evaluate(self, individual: Sequence | Solution) -> tuple[float]:
        """Evaluates the candidate individual with the information of the Travelling Salesmas Problem.

        The fitness of the solution is the fraction of the sum of the distances of the tour
        Args:
            individual (Sequence | Solution): Individual to evaluate

        Returns:
            Tuple[float]: Fitness
        """
        if len(individual) != self._nodes + 1:
            msg = f"Mismatch between individual variables ({len(individual)}) and instance variables ({self._nodes}) in {self.__class__.__name__}. A solution for the TSP must be a sequence of len {self._nodes + 1}"
            raise ValueError(msg)

        penalty: np.float64 = np.float64(0)

        if self.__evaluate_constraints(individual):
            distance: float = 0.0
            for i in range(len(individual) - 2):
                distance += self._distances[individual[i]][individual[i + 1]]

            fitness = 1.0 / distance
        else:
            fitness = 2.938736e-39  # --> 1.0 / np.float.max
            penalty = np.finfo(np.float64).max

        if isinstance(individual, Solution):
            individual.fitness = fitness
            individual.objectives = (fitness,)
            individual.constraints = (penalty,)

        return (fitness,)

    def __call__(self, individual: Sequence | Solution) -> tuple[float]:
        return self.evaluate(individual)

    def __repr__(self):
        return f"TSP<n={self._nodes}>"

    def __len__(self):
        return self._nodes

    def __array__(self, dtype=np.float32, copy: Optional[bool] = True) -> npt.ArrayLike:
        return np.asarray(self._coords, dtype=dtype, copy=copy)

    def create_solution(self) -> Solution:
        items = [0] + list(range(1, self._nodes)) + [0]
        return Solution(variables=items)

    def to_file(self, filename: str = "instance.tsp"):
        with open(filename, "w") as file:
            file.write(f"{len(self)}\n\n")
            content = "\n".join(f"{x}\t{y}" for (x, y) in self._coords)
            file.write(content)

    @classmethod
    def from_file(cls, filename: str) -> Self:
        # TODO: Improve using np.loadtxt
        with open(filename) as f:
            lines = f.readlines()
            lines = [line.rstrip() for line in lines]

        nodes = int(lines[0])
        coords = np.zeros(shape=(nodes, 2), dtype=np.float32)
        for i, line in enumerate(lines[2:]):
            x, y = line.split()
            coords[i] = [np.float32(x), np.float32(y)]

        return cls(nodes=nodes, coords=coords)

    def to_instance(self) -> Instance:
        return Instance(variables=self._coords.flatten())

class TSPDomain(Domain):
    """Domain to generate instances for the Symmetric Travelling Salesman Problem."""

    __FEAT_NAMES = "size,std_distances,centroid_x,centroid_y,radius,fraction_distances,area,variance_nnNds,variation_nnNds,cluster_ratio,mean_cluster_radius".split(
        ","
    )

    def __init__(
        self,
        dimension: int = 100,
        x_range: Tuple[int, int] = (0, 1000),
        y_range: Tuple[int, int] = (0, 1000),
        seed: int = 42,
    ):
        """Creates a new TSPDomain to generate instances for the Symmetric Travelling Salesman Problem

        Args:
            dimension (int, optional): Dimension of the instances to generate. Defaults to 100.
            x_range (Tuple[int, int], optional): Ranges for the Xs coordinates of each node/city. Defaults to (0, 1000).
            y_range (Tuple[int, int], optional): Ranges for the ys coordinates of each node/city. Defaults to (0, 1000).

        Raises:
            ValueError: If dimension is < 0
            ValueError: If x_range OR y_range does not have 2 dimensions each
            ValueError: If minimum ranges are greater than maximum ranges
        """
        if dimension < 0:
            raise ValueError(f"Expected dimension > 0 got {dimension}")
        if len(x_range) != 2 or len(y_range) != 2:
            raise ValueError(
                f"Expected x_range and y_range to be a tuple with only to integers. Got: x_range = {x_range} and y_range = {y_range}"
            )
        x_min, x_max = x_range
        y_min, y_max = y_range
        if x_min < 0 or x_max <= x_min:
            raise ValueError(
                f"Expected x_range to be (x_min, x_max) where x_min >= 0 and x_max > x_min. Got: x_range {x_range}"
            )
        if y_min < 0 or y_max <= y_min:
            raise ValueError(
                f"Expected y_range to be (y_min, y_max) where y_min >= 0 and y_max > y_min. Got: y_range {y_range}"
            )

        self._x_range = x_range
        self._y_range = y_range
        __bounds = [
            (x_min, x_max) if i % 2 == 0 else (y_min, y_max)
            for i in range(dimension * 2)
        ]

        super().__init__(dimension=dimension, bounds=__bounds, name="TSP", seed=seed)

    def generate_instances(self, n: int = 1) -> List[Instance]:
        """Generates N instances using numpy. It can return the instances in two formats:
        1. A numpy ndarray with the definition of the instances
        2. A list of Instance objects created from the raw numpy generation

        Args:
            n (int, optional): Number of instances to generate. Defaults to 1.
            cast (bool, optional): Whether to cast the raw data to Instance objects. Defaults to False.

        Returns:
            List[Instance]: Sequence of instances
        """
        instances = np.empty(shape=(n, self.dimension * 2), dtype=np.float32)
        instances[:, 0::2] = self._rng.uniform(
            low=self._x_range[0],
            high=self._x_range[1],
            size=(n, (self.dimension)),
        )
        instances[:, 1::2] = self._rng.uniform(
            low=self._y_range[0],
            high=self._y_range[1],
            size=(n, (self.dimension)),
        )
        return list(Instance(coords) for coords in instances)

    def extract_features(self, instances: Sequence[Instance]) -> np.ndarray:
        """Extract the features of the instance based on the TSP domain.
           For the TSP the features are:
            - Size
            - Standard deviation of the distances
            - Centroid coordinates
            - Radius of the instance
            - Fraction of distinct distances
            - Rectangular area
            - Variance of the normalised nearest neighbours distances
            - Coefficient of variation of the nearest neighbours distances
            - Cluster ratio
            - Mean cluster radius
        Args:
            instance (Instance): Instance to extract the features from

        Returns:
            Tuple[float]: Values of each feature
        """
        _instances = np.asarray(instances, copy=True)
        N_INSTANCES = len(_instances)
        N_CITIES = len(_instances[0]) // 2  # self.dimension // 2
        assert _instances is not instances
        coords = np.asarray(_instances, copy=True).reshape((N_INSTANCES, N_CITIES, 2))
        xs = coords[:, :, 0]
        ys = coords[:, :, 1]
        areas = (
            (np.max(xs, axis=1) - np.min(xs, axis=1))
            * (np.max(ys, axis=1) - np.min(ys, axis=1))
        ).astype(np.float64)

        # Compute distances for all instances
        distances = np.zeros((N_INSTANCES, N_CITIES, N_CITIES))
        differences = coords[:, :, np.newaxis, :] - coords[:, np.newaxis, :, :]
        distances = np.sqrt(np.sum(differences**2, axis=-1))
        mask = ~np.eye(N_CITIES, dtype=bool)
        std_distances = np.std(distances[:, mask], axis=1)

        centroids = np.mean(coords, axis=1)
        expanded_centroids = centroids[:, np.newaxis, :]
        centroids_distances = np.linalg.norm(coords - expanded_centroids, axis=-1)
        radius = np.mean(centroids_distances, axis=1)

        fractions = np.array(
            [
                np.unique(d[np.triu_indices_from(d, k=1)]).size
                / (N_CITIES * (N_CITIES - 1) / 2)
                for d in distances
            ]
        )
        # Top five only
        norm_distances = np.sort(distances, axis=2)[:, :, ::-1][:, :, :5] / np.max(
            distances, axis=(1, 2), keepdims=True
        )

        variance_nnds = np.var(norm_distances, axis=(1, 2))
        variation_nnds = variance_nnds / np.mean(norm_distances, axis=(1, 2))

        cluster_ratio = np.empty(shape=N_INSTANCES, dtype=np.float64)
        mean_cluster_radius = np.empty(shape=N_INSTANCES, dtype=np.float64)

        for i in range(N_INSTANCES):
            scale = np.mean(np.std(coords[i], axis=0))
            dbscan = DBSCAN(eps=0.2 * scale, min_samples=1)
            labels = dbscan.fit_predict(coords[i])
            unique_labels = [label for label in set(labels) if label != -1]
            cluster_ratio[i] = len(unique_labels) / N_CITIES
            # Cluster radius
            cluster_radius = np.empty(shape=len(unique_labels), dtype=np.float64)
            for j, label_id in enumerate(unique_labels):
                points_in_cluster = coords[i][labels == label_id]
                cluster_centroid = np.mean(points_in_cluster, axis=0)
                cluster_radius[j] = np.mean(
                    np.linalg.norm(points_in_cluster - cluster_centroid, axis=1)
                )

            mean_cluster_radius[i] = (
                np.mean(cluster_radius) if cluster_radius.size > 0 else 0.0
            )
        return np.column_stack(
            [
                np.full(shape=len(_instances), fill_value=N_CITIES),
                std_distances,
                centroids[:, 0],
                centroids[:, 1],
                radius,
                fractions,
                areas,
                variance_nnds,
                variation_nnds,
                cluster_ratio,
                mean_cluster_radius,
            ]
        ).astype(np.float64)

    def extract_features_as_dict(
        self, instances: Sequence[Instance]
    ) -> List[Dict[str, np.float32]]:
        """Creates a dictionary with the features of the instance.
        The key are the names of each feature and the values are
        the values extracted from instance.

        Args:
            instance (Instance): Instance to extract the features from

        Returns:
            Mapping[str, float]: Dictionary with the names/values of each feature
        """
        features = self.extract_features(instances)
        named_features: list[dict[str, np.float32]] = [{}] * len(features)
        for i, feats in enumerate(features):
            named_features[i] = {k: v for k, v in zip(TSPDomain.__FEAT_NAMES, feats)}
        return named_features

    def generate_problem_from_instance(self, instance: Instance) -> TSP:
        n_nodes = len(instance) // 2
        coords = np.array([*zip(instance[::2], instance[1::2])])
        return TSP(nodes=n_nodes, coords=coords)

    def generate_problems_from_instances(
        self, instances: Sequence[Instance]
    ) -> List[Problem]:
        if not isinstance(instances, np.ndarray):
            instances = np.asarray(instances)

        dimension = instances.shape[1] // 2
        return list(
            TSP(
                nodes=dimension, coords=np.array([*zip(instance[0::2], instance[1::2])])
            )
            for instance in instances
        )