Use Generators with a global prng within pymoo.

pymoo/__init__.py

@@ -1,3 +1,19 @@
 from pymoo.version import __version__
+import numpy as np
+import threading
 
+# Create the global random state Singleton
+class PymooPRNG(object):
+    _lock = threading.Lock()
+    _instance = None
 
+    def __new__(cls, seed=None):
+        if cls._instance is None: 
+            with cls._lock:
+                # Another thread could have created the instance
+                # before we acquired the lock. So check that the
+                # instance is still nonexistent.
+                if not cls._instance:
+                    cls._instance = np.random.default_rng(seed)
+
+        return cls._instance

pymoo/algorithms/hyperparameters.py

@@ -5,6 +5,7 @@
 from pymoo.core.parameters import get_params, flatten, set_params, hierarchical
 from pymoo.core.problem import ElementwiseProblem
 from pymoo.optimize import minimize
+import pymoo
 
 
 def create(algorithm, params):
@@ -77,7 +78,7 @@ def __init__(self, problem, n_runs=None, seeds=None, func_stats=stats_single_obj
             if n_runs is None:
                 raise Exception("Either provide number of runs or seeds directly.")
 
-            seeds = np.random.randint(1, 1000000, size=n_runs)
+            seeds = pymoo.PymooPRNG().integers(1, 1000000, size=n_runs)
 
         self.seeds = seeds
         self.func_stats = func_stats

pymoo/algorithms/moo/ctaea.py

@@ -17,6 +17,7 @@
 from pymoo.util.misc import has_feasible, random_permuations
 from pymoo.util.nds.non_dominated_sorting import NonDominatedSorting
 
+import pymoo
 
 # =========================================================================================================
 # Implementation
@@ -37,13 +38,13 @@ def comp_by_cv_dom_then_random(pop, P, **kwargs):
             elif rel == -1:
                 S[i] = b
             else:
-                S[i] = np.random.choice([a, b])
+                S[i] = pymoo.PymooPRNG().choice([a, b])
         elif pop[a].CV <= 0.0:
             S[i] = a
         elif pop[b].CV <= 0.0:
             S[i] = b
         else:
-            S[i] = np.random.choice([a, b])
+            S[i] = pymoo.PymooPRNG().choice([a, b])
 
     return S[:, None].astype(int)
 
@@ -68,7 +69,7 @@ def _do(self, problem, Hm, n_select, n_parents, **kwargs):
         if Pc <= Pd:
             # Choose from DA
             P[::n_parents, :] += n_pop
-        pf = np.random.random(n_select)
+        pf = pymoo.PymooPRNG().random(n_select)
         P[1::n_parents, :][pf >= Pc] += n_pop
 
         # compare using tournament function
@@ -175,7 +176,7 @@ def _updateCA(self, pop, n_survive):
                         if (delta_d[min_d_i] < 0) or (
                                 delta_d[min_d_i] == 0 and (FV[crowdest[list(min_d_i)]] > niche_worst).any()):
                             min_d_i = list(min_d_i)
-                            np.random.shuffle(min_d_i)
+                            pymoo.PymooPRNG().shuffle(min_d_i)
                             closest = crowdest[min_d_i]
                             niche_worst = closest[np.argmax(FV[closest])]
                         if FV[niche_worst] > worst_fit:

pymoo/algorithms/moo/dnsga2.py

@@ -2,7 +2,7 @@
 
 from pymoo.algorithms.moo.nsga2 import NSGA2
 from pymoo.core.population import Population
-
+import pymoo
 
 class DNSGA2(NSGA2):
 
@@ -32,7 +32,7 @@ def _advance(self, **kwargs):
         n_samples = int(np.ceil(len(pop) * self.perc_detect_change))
 
         # choose randomly some individuals of the current population to test if there was a change
-        I = np.random.choice(np.arange(len(pop)), size=n_samples)
+        I = pymoo.PymooPRNG().choice(np.arange(len(pop)), size=n_samples)
         samples = self.evaluator.eval(self.problem, Population.new(X=X[I]))
 
         # calculate the differences between the old and newly evaluated pop
@@ -47,7 +47,7 @@ def _advance(self, **kwargs):
             pop = Population.new(X=X)
 
             # find indices to be replaced (introduce diversity)
-            I = np.where(np.random.random(len(pop)) < self.perc_diversity)[0]
+            I = np.where(pymoo.PymooPRNG().random(len(pop)) < self.perc_diversity)[0]
 
             # replace with randomly sampled individuals
             if self.version == "A":

pymoo/algorithms/moo/moead.py

@@ -14,6 +14,7 @@
 from pymoo.util.display.multi import MultiObjectiveOutput
 from pymoo.util.reference_direction import default_ref_dirs
 
+import pymoo
 
 class NeighborhoodSelection(Selection):
 
@@ -28,10 +29,10 @@ def _do(self, problem, pop, n_select, n_parents, neighbors=None, **kwargs):
         prob = get(self.prob, size=n_select)
 
         for k in range(n_select):
-            if np.random.random() < prob[k]:
-                P[k] = np.random.choice(neighbors[k], n_parents, replace=False)
+            if pymoo.PymooPRNG().random() < prob[k]:
+                P[k] = pymoo.PymooPRNG().choice(neighbors[k], n_parents, replace=False)
             else:
-                P[k] = np.random.permutation(len(pop))[:n_parents]
+                P[k] = pymoo.PymooPRNG().permutation(len(pop))[:n_parents]
 
         return P
 
@@ -98,12 +99,13 @@ def _next(self):
         pop = self.pop
 
         # iterate for each member of the population in random order
-        for k in np.random.permutation(len(pop)):
+        for k in pymoo.PymooPRNG().permutation(len(pop)):
             # get the parents using the neighborhood selection
             P = self.selection.do(self.problem, pop, 1, self.mating.crossover.n_parents, neighbors=[self.neighbors[k]])
 
             # perform a mating using the default operators - if more than one offspring just pick the first
-            off = np.random.choice(self.mating.do(self.problem, pop, 1, parents=P, n_max_iterations=1))
+            # TODO: this is not just taking the first it's drawing one randomly - is this intended?
+            off = pymoo.PymooPRNG().choice(self.mating.do(self.problem, pop, 1, parents=P, n_max_iterations=1))
 
             # evaluate the offspring
             off = yield off
@@ -143,7 +145,7 @@ def _infill(self):
         pop_size, cross_parents, cross_off = self.pop_size, self.mating.crossover.n_parents, self.mating.crossover.n_offsprings
 
         # do the mating in a random order
-        indices = np.random.permutation(len(self.pop))[:self.n_offsprings]
+        indices = pymoo.PymooPRNG().permutation(len(self.pop))[:self.n_offsprings]
 
         # get the parents using the neighborhood selection
         P = self.selection.do(self.problem, self.pop, self.n_offsprings, cross_parents,
@@ -153,7 +155,7 @@ def _infill(self):
         off = self.mating.do(self.problem, self.pop, 1e12, n_max_iterations=1, parents=P)
 
         # select a random offspring from each mating
-        off = Population.create(*[np.random.choice(pool) for pool in np.reshape(off, (self.n_offsprings, -1))])
+        off = Population.create(*[pymoo.PymooPRNG().choice(pool) for pool in np.reshape(off, (self.n_offsprings, -1))])
 
         # store the indices because of the neighborhood matching in advance
         self.indices = indices

pymoo/algorithms/moo/nsga3.py

@@ -15,6 +15,7 @@
 from pymoo.util.misc import intersect, has_feasible
 from pymoo.util.nds.non_dominated_sorting import NonDominatedSorting
 
+import pymoo
 
 # =========================================================================================================
 # Implementation
@@ -32,7 +33,7 @@ def comp_by_cv_then_random(pop, P, **kwargs):
 
         # both solutions are feasible just set random
         else:
-            S[i] = np.random.choice([a, b])
+            S[i] = pymoo.PymooPRNG().choice([a, b])
 
     return S[:, None].astype(int)
 
@@ -214,15 +215,15 @@ def niching(pop, n_remaining, niche_count, niche_of_individuals, dist_to_niche):
 
         # all niches with the minimum niche count (truncate if randomly if more niches than remaining individuals)
         next_niches = next_niches_list[np.where(next_niche_count == min_niche_count)[0]]
-        next_niches = next_niches[np.random.permutation(len(next_niches))[:n_select]]
+        next_niches = next_niches[pymoo.PymooPRNG().permutation(len(next_niches))[:n_select]]
 
         for next_niche in next_niches:
 
             # indices of individuals that are considered and assign to next_niche
             next_ind = np.where(np.logical_and(niche_of_individuals == next_niche, mask))[0]
 
             # shuffle to break random tie (equal perp. dist) or select randomly
-            np.random.shuffle(next_ind)
+            pymoo.PymooPRNG().shuffle(next_ind)
 
             if niche_count[next_niche] == 0:
                 next_ind = next_ind[np.argmin(dist_to_niche[next_ind])]

pymoo/algorithms/moo/sms.py

@@ -17,6 +17,7 @@
 from pymoo.util.nds.non_dominated_sorting import NonDominatedSorting
 from pymoo.util.normalization import normalize
 
+import pymoo
 
 # ---------------------------------------------------------------------------------------------------------
 # Environmental Survival - Remove the solution with the least HV contribution
@@ -122,7 +123,7 @@ def cv_and_dom_tournament(pop, P, *args, **kwargs):
 
             # if rank or domination relation didn't make a decision compare by crowding
             if np.isnan(S[i]):
-                S[i] = np.random.choice([a, b])
+                S[i] = pymoo.PymooPRNG().choice([a, b])
 
     return S[:, None].astype(int, copy=False)
 

pymoo/algorithms/moo/unsga3.py

@@ -3,6 +3,7 @@
 from pymoo.algorithms.moo.nsga3 import NSGA3
 from pymoo.operators.selection.tournament import compare, TournamentSelection
 
+import pymoo
 
 # =========================================================================================================
 # Implementation
@@ -33,7 +34,7 @@ def comp_by_rank_and_ref_line_dist(pop, P, **kwargs):
                                    method='smaller_is_better')
 
         if np.isnan(S[i]):
-            S[i] = np.random.choice([a, b])
+            S[i] = pymoo.PymooPRNG().choice([a, b])
 
     return S[:, None].astype(int)
 

pymoo/algorithms/soo/nonconvex/brkga.py

@@ -15,6 +15,7 @@
 from pymoo.util.display.single import SingleObjectiveOutput
 from pymoo.util.nds.non_dominated_sorting import NonDominatedSorting
 
+import pymoo
 
 # =========================================================================================================
 # Implementation
@@ -68,8 +69,8 @@ def _do(self, problem, pop, n_select, n_parents, **kwargs):
             non_elites = elites
 
         # do the mating selection - always one elite and one non-elites
-        s_elite = np.random.choice(elites, size=n_select)
-        s_non_elite = np.random.choice(non_elites, size=n_select)
+        s_elite = pymoo.PymooPRNG().choice(elites, size=n_select)
+        s_non_elite = pymoo.PymooPRNG().choice(non_elites, size=n_select)
 
         return np.column_stack([s_elite, s_non_elite])
 

pymoo/algorithms/soo/nonconvex/cmaes.py

@@ -13,6 +13,7 @@
 from pymoo.util.optimum import filter_optimum
 from pymoo.vendor.vendor_cmaes import my_fmin
 
+import pymoo
 
 # =========================================================================================================
 # Implementation
@@ -282,6 +283,7 @@ def __init__(self,
         popsize : 4+int(3*np.log(N))
             Population size, AKA lambda, number of new solution per iteration
 
+        # TODO is randn actually something that we can pass? doesn't appear to be used in code base.
         randn : np.random.randn
             Randn(lam, N) must return an np.array of shape (lam, N), see also cma.utilities.math.randhss
 

pymoo/algorithms/soo/nonconvex/de.py

@@ -38,6 +38,7 @@
 from pymoo.util.display.single import SingleObjectiveOutput
 from pymoo.util.misc import where_is_what
 
+import pymoo
 
 # =========================================================================================================
 # Crossover
@@ -54,7 +55,7 @@ def de_differential(X, F, jitter, alpha=0.001):
     for i in range(1, n_parents, 2):
         # create the weight vectors with jitter to give some variation
         _F = F[:, None].repeat(n_var, axis=1)
-        _F[jitter] *= (1 + alpha * (np.random.random((jitter.sum(), n_var)) - 0.5))
+        _F[jitter] *= (1 + alpha * (pymoo.PymooPRNG().random((jitter.sum(), n_var)) - 0.5))
 
         # add the difference to the vector
         delta += _F * (X[i] - X[i + 1])
@@ -129,7 +130,7 @@ def do(self, problem, pop, n_offsprings, algorithm=None, **kwargs):
 
             itself = np.array(targets)[:, None]
 
-            best = lambda: np.random.choice(np.where(pop.get("rank") == 0)[0], replace=True, size=n_matings)
+            best = lambda: pymoo.PymooPRNG().choice(np.where(pop.get("rank") == 0)[0], replace=True, size=n_matings)
 
             if sel_type == "rand":
                 fast_fill_random(P, len(pop), columns=range(n_parents), Xp=itself)
@@ -175,10 +176,10 @@ def do(self, problem, pop, n_offsprings, algorithm=None, **kwargs):
                 _trial = np.copy(_target)
                 _trial[M] = _donor[M]
             elif name == "line":
-                w = np.random.random((len(K), 1)) * _CR[:, None]
+                w = pymoo.PymooPRNG().random((len(K), 1)) * _CR[:, None]
                 _trial = _target + w * (_donor - _target)
             elif name == "hypercube":
-                w = np.random.random((len(K), _target.shape[1])) * _CR[:, None]
+                w = pymoo.PymooPRNG().random((len(K), _target.shape[1])) * _CR[:, None]
                 _trial = _target + w * (_donor - _target)
             else:
                 raise Exception(f"Unknown crossover variant: {name}")
@@ -252,7 +253,7 @@ def _infill(self):
 
         # if number of offsprings is set lower than pop_size - randomly select
         if self.n_offsprings < self.pop_size:
-            index = np.random.permutation(len(infills))[:self.n_offsprings]
+            index = pymoo.PymooPRNG().permutation(len(infills))[:self.n_offsprings]
             infills = infills[index]
 
         infills.set("index", index)

pymoo/algorithms/soo/nonconvex/es.py

@@ -9,6 +9,7 @@
 from pymoo.util.display.single import SingleObjectiveOutput
 from pymoo.util.optimum import filter_optimum
 
+import pymoo
 
 class ES(GeneticAlgorithm):
 
@@ -96,7 +97,7 @@ def _infill(self):
         sigmap = np.minimum(self.sigma_max, es_sigma(sigmap, self.tau, self.taup))
 
         # execute the evolutionary strategy to calculate the offspring solutions
-        Xp = X + sigmap * np.random.normal(size=sigmap.shape)
+        Xp = X + sigmap * pymoo.PymooPRNG().normal(size=sigmap.shape)
 
         # if gamma is not none do the differential variation overwrite Xp and sigmap for the first mu-1 individuals
         if self.gamma is not None:
@@ -127,7 +128,7 @@ def _set_optimum(self):
 
 def es_sigma(sigma, tau, taup):
     _lambda, _n = sigma.shape
-    return sigma * np.exp(taup * np.random.normal(size=(_lambda, 1)) + tau * np.random.normal(size=(_lambda, _n)))
+    return sigma * np.exp(taup * pymoo.PymooPRNG().normal(size=(_lambda, 1)) + tau * pymoo.PymooPRNG().normal(size=(_lambda, _n)))
 
 
 def es_intermediate_recomb(sigma):
@@ -136,7 +137,7 @@ def es_intermediate_recomb(sigma):
 
     for i in range(_lambda):
         for j in range(_n):
-            k = np.random.randint(_lambda)
+            k = pymoo.PymooPRNG().integers(_lambda)
             sigma_hat[i, j] = (sigma[i, j] + sigma[k, j]) / 2.0
 
     return sigma_hat
@@ -160,7 +161,7 @@ def es_mut_repair(Xp, X, sigma, xl, xu, n_trials):
             break
         else:
             # do the mutation again vectored for all values not in bound
-            Xp[i, j] = X[i, j] + sigma[i, j] * np.random.normal(size=len(i))
+            Xp[i, j] = X[i, j] + sigma[i, j] * pymoo.PymooPRNG().normal(size=len(i))
 
     # if there are still solutions which boundaries are violated, set them to the original X
     if not all_in_bounds:
@@ -189,7 +190,7 @@ def es_mut_loop(X, sigmap, xl, xu, n_trials=10):
             for _ in range(n_trials):
 
                 # calculate the mutated value
-                x = X[i, j] + sigmap[i, j] * np.random.normal()
+                x = X[i, j] + sigmap[i, j] * pymoo.PymooPRNG().normal()
 
                 # if it is inside the bounds accept it - otherwise try again
                 if xl[j] <= x <= xu[j]:

pymoo/algorithms/soo/nonconvex/g3pcx.py

@@ -14,6 +14,7 @@
 from pymoo.operators.selection.rnd import fast_fill_random
 from pymoo.util.display.single import SingleObjectiveOutput
 
+import pymoo
 
 # =========================================================================================================
 # Implementation
@@ -77,7 +78,7 @@ def _next(self, **kwargs):
 
             pop, family_size = self.pop, get(self.family_size)
 
-            rnd = np.random.choice(np.arange(len(pop)), size=family_size, replace=False)
+            rnd = pymoo.PymooPRNG().choice(np.arange(len(pop)), size=family_size, replace=False)
             family = Population.merge(pop[rnd], off)
             pop[rnd] = FitnessSurvival().do(self.problem, family, n_survive=family_size)