```import sys

import numpy as np
from scipy import special

from pymoo.util.misc import find_duplicates, cdist

# =========================================================================================================
# Model
# =========================================================================================================

class ReferenceDirectionFactory:

def __init__(self, n_dim, scaling=None, lexsort=True, verbose=False, seed=None, **kwargs) -> None:
super().__init__()
self.n_dim = n_dim
self.scaling = scaling
self.lexsort = lexsort
self.verbose = verbose
self.seed = seed

def do(self):

# set the random seed if it is provided
if self.seed is not None:
np.random.seed(self.seed)

if self.n_dim == 1:
return np.array([[1.0]])
else:

val = self._do()
if isinstance(val, tuple):
ref_dirs, other = val[0], val[1:]
else:
ref_dirs = val

if self.scaling is not None:
ref_dirs = scale_reference_directions(ref_dirs, self.scaling)

# do ref_dirs is desired
if self.lexsort:
I = np.lexsort([ref_dirs[:, j] for j in range(ref_dirs.shape[1])][::-1])
ref_dirs = ref_dirs[I]

return ref_dirs

def _do(self):
return None

# =========================================================================================================
# Das Dennis Reference Directions (Uniform)
# =========================================================================================================

def get_number_of_uniform_points(n_partitions, n_dim):
"""
Returns the number of uniform points that can be created uniformly.
"""
return int(special.binom(n_dim + n_partitions - 1, n_partitions))

def get_partition_closest_to_points(n_points, n_dim):
"""
Returns the corresponding partition number which create the desired number of points
or less!
"""

if n_dim == 1:
return 0

n_partitions = 1
_n_points = get_number_of_uniform_points(n_partitions, n_dim)
while _n_points <= n_points:
n_partitions += 1
_n_points = get_number_of_uniform_points(n_partitions, n_dim)
return n_partitions - 1

def das_dennis(n_partitions, n_dim):
if n_partitions == 0:
return np.full((1, n_dim), 1 / n_dim)
else:
ref_dirs = []
ref_dir = np.full(n_dim, np.nan)
das_dennis_recursion(ref_dirs, ref_dir, n_partitions, n_partitions, 0)
return np.concatenate(ref_dirs, axis=0)

def das_dennis_recursion(ref_dirs, ref_dir, n_partitions, beta, depth):
if depth == len(ref_dir) - 1:
ref_dir[depth] = beta / (1.0 * n_partitions)
ref_dirs.append(ref_dir[None, :])
else:
for i in range(beta + 1):
ref_dir[depth] = 1.0 * i / (1.0 * n_partitions)
das_dennis_recursion(ref_dirs, np.copy(ref_dir), n_partitions, beta - i, depth + 1)

class UniformReferenceDirectionFactory(ReferenceDirectionFactory):

def __init__(self, n_dim, scaling=None, n_points=None, n_partitions=None, **kwargs) -> None:
super().__init__(n_dim, scaling=scaling, **kwargs)

if n_points is not None:
n_partitions = get_partition_closest_to_points(n_points, n_dim)
results_in = get_number_of_uniform_points(n_partitions, n_dim)

# the number of points are not matching to any partition number
if results_in != n_points:
results_in_next = get_number_of_uniform_points(n_partitions + 1, n_dim)
raise Exception("The number of points (n_points = %s) can not be created uniformly.\n"
"Either choose n_points = %s (n_partitions = %s) or "
"n_points = %s (n_partitions = %s)." %
(n_points, results_in, n_partitions, results_in_next, n_partitions + 1))

self.n_partitions = n_partitions

elif n_partitions is not None:
self.n_partitions = n_partitions

else:
raise Exception("Either provide number of partitions or number of points.")

def _do(self):
return das_dennis(self.n_partitions, self.n_dim)

# =========================================================================================================
# Multi Layer
# =========================================================================================================

class MultiLayerReferenceDirectionFactory:

def __init__(self, *args) -> None:
self.layers = []
self.layers.extend(args)

self.layers.extend(args)

def do(self):
ref_dirs = []
for factory in self.layers:
ref_dirs.append(factory)
ref_dirs = np.concatenate(ref_dirs, axis=0)
is_duplicate = find_duplicates(ref_dirs)
return ref_dirs[np.logical_not(is_duplicate)]

# =========================================================================================================
# Util
# =========================================================================================================

def sample_on_unit_simplex(n_points, n_dim, unit_simplex_mapping="kraemer"):

if unit_simplex_mapping == "sum":
rnd = map_onto_unit_simplex(np.random.random((n_points, n_dim)), "sum")

elif unit_simplex_mapping == "kraemer":
rnd = map_onto_unit_simplex(np.random.random((n_points, n_dim)), "kraemer")

elif unit_simplex_mapping == "das-dennis":
n_partitions = get_partition_closest_to_points(n_points, n_dim)
rnd = UniformReferenceDirectionFactory(n_dim, n_partitions=n_partitions).do()

else:
raise Exception("Please define a valid sampling on unit simplex strategy!")

return rnd

def map_onto_unit_simplex(rnd, method):
n_points, n_dim = rnd.shape

if method == "sum":
ret = rnd / rnd.sum(axis=1)[:, None]

elif method == "kraemer":
M = sys.maxsize

rnd *= M
rnd = rnd[:, :n_dim - 1]
rnd = np.column_stack([np.zeros(n_points), rnd, np.full(n_points, M)])

rnd = np.sort(rnd, axis=1)

ret = np.full((n_points, n_dim), np.nan)
for i in range(1, n_dim + 1):
ret[:, i - 1] = rnd[:, i] - rnd[:, i - 1]
ret /= M

else:
raise Exception("Invalid unit simplex mapping!")

return ret

def scale_reference_directions(ref_dirs, scaling):
return ref_dirs * scaling + ((1 - scaling) / ref_dirs.shape[1])

def select_points_with_maximum_distance(X, n_select, selected=[]):
n_points, n_dim = X.shape

# calculate the distance matrix
D = cdist(X, X)

# if no selection provided pick randomly in the beginning
if len(selected) == 0:
selected = [np.random.randint(len(X))]

# create variables to store what selected and what not
not_selected = [i for i in range(n_points) if i not in selected]

# remove unnecessary points
dist_to_closest_selected = D[:, selected].min(axis=1)

# now select the points until sufficient ones are found
while len(selected) < n_select:
# find point that has the maximum distance to all others
index_in_not_selected = dist_to_closest_selected[not_selected].argmax()
I = not_selected[index_in_not_selected]

# add the closest distance to selected point
is_closer = D[I] < dist_to_closest_selected
dist_to_closest_selected[is_closer] = D[I][is_closer]

# add it to the selected and remove from not selected
selected.append(I)
not_selected = np.delete(not_selected, index_in_not_selected)

return selected```