Setting precision for all calculations

Dockerfile

@@ -0,0 +1,11 @@
+FROM ubuntu:22.04
+
+ENV TZ="America/New_York"
+RUN ln -snf /usr/share/zoneinfo/${TZ} /etc/localtime && echo ${TZ} > /etc/timezone
+
+WORKDIR /opt/app
+ADD . /opt/app
+
+RUN apt update
+RUN apt install -y python3 python3-pip
+RUN python3 -m pip install -r requirements.txt

crfmnes/alg.py

@@ -3,19 +3,18 @@
 
 import math
 import numpy as np
-import copy
 
 # evaluation value of the infeasible solution
 INFEASIBLE = np.inf
 
 
 def get_h_inv(dim):
-    f = lambda a, b: ((1. + a * a) * math.exp(a * a / 2.) / 0.24) - 10. - dim
-    f_prime = lambda a: (1. / 0.24) * a * math.exp(a * a / 2.) * (3. + a * a)
-    h_inv = 6.0
+    f = lambda a, b: ((1. + a**2.) * math.exp(a**2. / 2.) / .24) - 10. - b
+    f_prime = lambda a: a * math.exp(a**2. / 2.) * (3. + a**2.) / .24
+    h_inv = 6.
     while abs(f(h_inv, dim)) > 1e-10:
         last = h_inv
-        h_inv = h_inv - 0.5 * (f(h_inv, dim) / f_prime(h_inv))
+        h_inv = h_inv - .5 * (f(h_inv, dim) / f_prime(h_inv))
         if abs(h_inv - last) < 1e-16:
             # Exit early since no further improvements are happening
             break
@@ -29,7 +28,7 @@ def sort_indices_by(evals, z):
     no_of_feasible_solutions = np.where(sorted_evals != INFEASIBLE)[0].size
     if no_of_feasible_solutions != lam:
         infeasible_z = z[:, np.where(evals == INFEASIBLE)[0]]
-        distances = np.sum(infeasible_z ** 2, axis=0)
+        distances = np.sum(infeasible_z**2, axis=0)
         infeasible_indices = sorted_indices[no_of_feasible_solutions:]
         indices_sorted_by_distance = np.argsort(distances)
         sorted_indices[no_of_feasible_solutions:] = infeasible_indices[indices_sorted_by_distance]
@@ -39,60 +38,66 @@ def sort_indices_by(evals, z):
 class CRFMNES:
     def __init__(self, dim, f, m, sigma, lamb, **kwargs):
 
+        self.t = kwargs.get('dtype', np.float64)
+        t = lambda v: self.t(v)
+        self._1 = t(1.)
+        self._2 = t(2.)
+        self._dim = t(dim)
+
         if 'seed' in kwargs.keys():
             np.random.seed(kwargs['seed'])
         self.dim = dim
         self.f = f
-        self.m = m
-        self.sigma = sigma
+        self.m = t(m)
+        self.sigma = t(sigma)
         self.lamb = lamb
         assert (lamb > 0 and lamb % 2 == 0), f"The value of 'lamb' must be an even, positive integer greater than 0"
 
-        self.v = kwargs.get('v', np.random.randn(dim, 1) / np.sqrt(dim))
-        self.D = np.ones([dim, 1])
-        self.constraint = kwargs.get('constraint', [[- np.inf, np.inf] for _ in range(dim)])
-        self.penalty_coef = kwargs.get('penalty_coef', 1e5)
+        self.v = kwargs.get('v', np.random.randn(dim, 1) / np.sqrt(dim)).astype(self.t)
+        self.D = np.ones([dim, 1], dtype=self.t)
+        self.constraint = np.asarray(kwargs.get('constraint', [[-t('inf'), +t('inf')] for _ in range(dim)]), dtype=self.t)
+        self.penalty_coef = t(kwargs.get('penalty_coef', 1e5))
         self.use_constraint_violation = kwargs.get('use_constraint_violation', True)
 
-        self.w_rank_hat = (np.log(self.lamb / 2 + 1) - np.log(np.arange(1, self.lamb + 1))).reshape(self.lamb, 1)
-        self.w_rank_hat[np.where(self.w_rank_hat < 0)] = 0
-        self.w_rank = self.w_rank_hat / sum(self.w_rank_hat) - (1. / self.lamb)
-        self.mueff = 1 / ((self.w_rank + (1 / self.lamb)).T @ (self.w_rank + (1 / self.lamb)))[0][0]
-        self.cs = (self.mueff + 2.) / (self.dim + self.mueff + 5.)
-        self.cc = (4. + self.mueff / self.dim) / (self.dim + 4. + 2. * self.mueff / self.dim)
-        self.c1_cma = 2. / (math.pow(self.dim + 1.3, 2) + self.mueff)
+        self.w_rank_hat = (np.log(self.lamb / 2 + 1, dtype=self.t) - np.log(np.arange(1, self.lamb + 1), dtype=self.t)).reshape(self.lamb, 1)
+        self.w_rank_hat[np.where(self.w_rank_hat < 0.)] = 0.
+        self.w_rank = self.w_rank_hat / np.sum(self.w_rank_hat) - (1. / self.lamb)
+        self.mueff = np.reciprocal(((self.w_rank + (1. / self.lamb)).T @ (self.w_rank + (1. / self.lamb)))[0][0])
+        self.cs = (self.mueff + self._2) / (self.mueff + self._dim + t(5.))
+        self.cc = (self.mueff / self._dim + t(4.)) / (self.mueff * t(2.) / self._dim + self._dim + t(4.))
+        self.c1_cma = np.reciprocal(np.power(self._dim + 1.3, 2, dtype=self.t) + self.mueff) * self._2
         # initialization
-        self.chiN = np.sqrt(self.dim) * (1. - 1. / (4. * self.dim) + 1. / (21. * self.dim * self.dim))
-        self.pc = np.zeros([self.dim, 1])
-        self.ps = np.zeros([self.dim, 1])
+        self.chiN = np.sqrt(self._dim, dtype=self.t) * (np.reciprocal(21 * self.dim**2, dtype=self.t) - np.reciprocal(4 * self.dim, dtype=self.t) + self._1)
+        self.pc = np.zeros([self.dim, 1], dtype=self.t)
+        self.ps = np.zeros([self.dim, 1], dtype=self.t)
         # distance weight parameter
-        self.h_inv = get_h_inv(self.dim)
-        self.alpha_dist = lambda lambF: self.h_inv * min(1., math.sqrt(float(self.lamb) / self.dim)) * math.sqrt(
-            float(lambF) / self.lamb)
-        self.w_dist_hat = lambda z, lambF: math.exp(self.alpha_dist(lambF) * np.linalg.norm(z))
+        self.h_inv = t(get_h_inv(self.dim))
+        self.alpha_dist = lambda lambF: self.h_inv * min(1., np.sqrt(t(self.lamb) / self._dim)) * np.sqrt(
+            t(lambF) / self.lamb)
+        self.w_dist_hat = lambda z, lambF: np.exp(self.alpha_dist(lambF) * np.linalg.norm(z))
         # learning rate
-        self.eta_m = 1.0
-        self.eta_move_sigma = 1.
-        self.eta_stag_sigma = lambda lambF: math.tanh((0.024 * lambF + 0.7 * self.dim + 20.) / (self.dim + 12.))
-        self.eta_conv_sigma = lambda lambF: 2. * math.tanh((0.025 * lambF + 0.75 * self.dim + 10.) / (self.dim + 4.))
-        self.c1 = lambda lambF: self.c1_cma * (self.dim - 5) / 6 * (float(lambF) / self.lamb)
-        self.eta_B = lambda lambF: np.tanh((min(0.02 * lambF, 3 * np.log(self.dim)) + 5) / (0.23 * self.dim + 25))
+        self.eta_m = self._1
+        self.eta_move_sigma = self._1
+        self.eta_stag_sigma = lambda lambF: np.tanh((t(lambF) * .024 + t(.70) * self._dim + 20.) / t(self.dim + 12), dtype=self.t)
+        self.eta_conv_sigma = lambda lambF: np.tanh((t(lambF) * .025 + t(.75) * self._dim + 10.) / t(self.dim +  4), dtype=self.t) * self._2
+        self.c1 = lambda lambF: self.c1_cma * (self.dim - 5) / 6. * (t(lambF) / self.lamb)
+        self.eta_B = lambda lambF: np.tanh((min(t(lambF) * .02, np.log(self._dim, dtype=self.t) * 3.) + 5.) / t(.23 * self._dim + 25.), dtype=self.t)
 
         self.g = 0
         self.no_of_evals = 0
 
         self.idxp = np.arange(self.lamb / 2, dtype=int)
         self.idxm = np.arange(self.lamb / 2, self.lamb, dtype=int)
-        self.z = np.zeros([self.dim, self.lamb])
+        self.z = np.zeros([self.dim, self.lamb], dtype=self.t)
 
-        self.f_best = float('inf')
-        self.x_best = np.empty(self.dim)
+        self.f_best = t('inf')
+        self.x_best = np.empty(self.dim, dtype=self.t)
 
     def calc_violations(self, x):
-        violations = np.zeros(self.lamb)
+        violations = np.zeros(self.lamb, dtype=self.t)
         for i in range(self.lamb):
             for j in range(self.dim):
-                violations[i] += (- min(0, x[j][i] - self.constraint[j][0]) + max(0, x[j][i] - self.constraint[j][1])) * self.penalty_coef
+                violations[i] += (-min(0, x[j][i] - self.constraint[j][0]) + max(0, x[j][i] - self.constraint[j][1])) * self.penalty_coef
         return violations
 
     def optimize(self, iterations):
@@ -104,18 +109,20 @@ def optimize(self, iterations):
     def one_iteration(self):
         d = self.dim
         lamb = self.lamb
-        zhalf = np.random.randn(d, int(lamb / 2))  # dim x lamb/2
-        self.z[:, self.idxp] = zhalf
+        zhalf = np.random.randn(d, int(lamb / 2)).astype(self.t)  # dim x lamb/2
+        self.z[:, self.idxp] = +zhalf
         self.z[:, self.idxm] = -zhalf
-        normv = np.linalg.norm(self.v)
-        normv2 = normv ** 2
+        normv = self.t(np.linalg.norm(self.v))
+        normv2 = np.square(normv, dtype=self.t)
         vbar = self.v / normv
-        y = self.z + (np.sqrt(1 + normv2) - 1) * vbar @ (vbar.T @ self.z)
+        y = self.z + (np.sqrt(normv2 + self._1) - self._1) * vbar @ (vbar.T @ self.z)
         x = self.m + self.sigma * y * self.D
-        evals_no_sort = np.array([self.f(np.array(x[:, i].reshape(self.dim, 1))) for i in range(self.lamb)])
+        evals_no_sort = np.array(
+            [self.f(np.array(x[:, i].reshape(self.dim, 1), dtype=self.t)) for i in range(self.lamb)],
+            dtype=self.t)
         xs_no_sort = [x[:, i] for i in range(lamb)]
 
-        violations = np.zeros(lamb)
+        violations = np.zeros(lamb, dtype=self.t)
         if self.use_constraint_violation:
             violations = self.calc_violations(x)
             sorted_indices = sort_indices_by(evals_no_sort + violations, self.z)
@@ -135,58 +142,58 @@ def one_iteration(self):
             self.x_best = x_best
 
         # This operation assumes that if the solution is infeasible, infinity comes in as input.
-        lambF = np.sum(evals_no_sort < np.finfo(float).max)
+        lambF = np.sum(evals_no_sort < np.finfo(self.t).max)
 
         # evolution path p_sigma
-        self.ps = (1 - self.cs) * self.ps + np.sqrt(self.cs * (2. - self.cs) * self.mueff) * (self.z @ self.w_rank)
+        self.ps = (self._1 - self.cs) * self.ps + np.sqrt(self.cs * (self._2 - self.cs) * self.mueff) * (self.z @ self.w_rank)
         ps_norm = np.linalg.norm(self.ps)
         # distance weight
         w_tmp = np.array(
-            [self.w_rank_hat[i] * self.w_dist_hat(np.array(self.z[:, i]), lambF) for i in range(self.lamb)]).reshape(
-            self.lamb, 1)
-        weights_dist = w_tmp / sum(w_tmp) - 1. / self.lamb
+            [self.w_rank_hat[i] * self.w_dist_hat(np.array(self.z[:, i], dtype=self.t), lambF) for i in range(self.lamb)],
+            dtype=self.t).reshape(self.lamb, 1)
+        weights_dist = w_tmp / np.sum(w_tmp) - self._1 / self.lamb
         # switching weights and learning rate
         weights = weights_dist if ps_norm >= self.chiN else self.w_rank
         eta_sigma = self.eta_move_sigma if ps_norm >= self.chiN else self.eta_stag_sigma(
-            lambF) if ps_norm >= 0.1 * self.chiN else self.eta_conv_sigma(lambF)
+            lambF) if ps_norm >= .1 * self.chiN else self.eta_conv_sigma(lambF)
         # update pc, m
         wxm = (x - self.m) @ weights
-        self.pc = (1. - self.cc) * self.pc + np.sqrt(self.cc * (2. - self.cc) * self.mueff) * wxm / self.sigma
+        self.pc = (self._1 - self.cc) * self.pc + np.sqrt(self.cc * (self._2 - self.cc) * self.mueff) * wxm / self.sigma
         self.m += self.eta_m * wxm
         # calculate s, t
         # step1
-        normv4 = normv2 ** 2
+        normv4 = np.square(normv2, dtype=self.t)
         exY = np.append(y, self.pc / self.D, axis=1)  # dim x lamb+1
         yy = exY * exY  # dim x lamb+1
         ip_yvbar = vbar.T @ exY
         yvbar = exY * vbar  # dim x lamb+1. exYのそれぞれの列にvbarがかかる
-        gammav = 1. + normv2
+        gammav = normv2 + self._1
         vbarbar = vbar * vbar
         alphavd = np.min(
-            [1, np.sqrt(normv4 + (2 * gammav - np.sqrt(gammav)) / np.max(vbarbar)) / (2 + normv2)])  # scalar
-        t = exY * ip_yvbar - vbar * (ip_yvbar ** 2 + gammav) / 2  # dim x lamb+1
-        b = -(1 - alphavd ** 2) * normv4 / gammav + 2 * alphavd ** 2
-        H = np.ones([self.dim, 1]) * 2 - (b + 2 * alphavd ** 2) * vbarbar  # dim x 1
-        invH = H ** (-1)
-        s_step1 = yy - normv2 / gammav * (yvbar * ip_yvbar) - np.ones([self.dim, self.lamb + 1])  # dim x lamb+1
+            [self._1, np.sqrt(normv4 + (gammav * self._2 - np.sqrt(gammav)) / np.max(vbarbar)) / (normv2 + self._2)])  # scalar
+        t = exY * ip_yvbar - vbar * (ip_yvbar ** self._2 + gammav) / self._2  # dim x lamb+1
+        b = (alphavd ** self._2 - self._1) * normv4 / gammav + self._2 * alphavd ** self._2
+        H = np.full([self.dim, 1], self._2, dtype=self.t) - (b + self._2 * alphavd ** self._2) * vbarbar  # dim x 1
+        invH = H ** -self._1
+        s_step1 = yy - normv2 / gammav * (yvbar * ip_yvbar) - np.ones([self.dim, self.lamb + 1], dtype=self.t)  # dim x lamb+1
         ip_vbart = vbar.T @ t  # 1 x lamb+1
-        s_step2 = s_step1 - alphavd / gammav * ((2 + normv2) * (t * vbar) - normv2 * vbarbar @ ip_vbart)  # dim x lamb+1
+        s_step2 = s_step1 - alphavd / gammav * ((self._2 + normv2) * (t * vbar) - normv2 * vbarbar @ ip_vbart)  # dim x lamb+1
         invHvbarbar = invH * vbarbar
         ip_s_step2invHvbarbar = invHvbarbar.T @ s_step2  # 1 x lamb+1
         s = (s_step2 * invH) - b / (
-                    1 + b * vbarbar.T @ invHvbarbar) * invHvbarbar @ ip_s_step2invHvbarbar  # dim x lamb+1
+                    self._1 + b * vbarbar.T @ invHvbarbar) * invHvbarbar @ ip_s_step2invHvbarbar  # dim x lamb+1
         ip_svbarbar = vbarbar.T @ s  # 1 x lamb+1
-        t = t - alphavd * ((2 + normv2) * (s * vbar) - vbar @ ip_svbarbar)  # dim x lamb+1
+        t = t - alphavd * ((self._2 + normv2) * (s * vbar) - vbar @ ip_svbarbar)  # dim x lamb+1
         # update v, D
-        exw = np.append(self.eta_B(lambF) * weights, np.array([self.c1(lambF)]).reshape(1, 1),
+        exw = np.append(self.eta_B(lambF) * weights, np.array([self.c1(lambF)], dtype=self.t).reshape(1, 1),
                         axis=0)  # lamb+1 x 1
         self.v = self.v + (t @ exw) / normv
         self.D = self.D + (s @ exw) * self.D
         # calculate detA
-        nthrootdetA = np.exp(np.sum(np.log(self.D)) / self.dim + np.log(1 + self.v.T @ self.v) / (2 * self.dim))[0][0]
+        nthrootdetA = np.exp(np.sum(np.log(self.D)) / self._dim + np.log(1 + self.v.T @ self.v) / (2 * self.dim))[0][0]
         self.D = self.D / nthrootdetA
         # update sigma
-        G_s = np.sum((self.z * self.z - np.ones([self.dim, self.lamb])) @ weights) / self.dim
-        self.sigma = self.sigma * np.exp(eta_sigma / 2 * G_s)
+        G_s = np.sum((self.z * self.z - np.ones([self.dim, self.lamb], dtype=self.t)) @ weights) / self._dim
+        self.sigma = self.sigma * np.exp(eta_sigma / self._2 * G_s, dtype=self.t)
 
         return xs_no_sort, evals_no_sort, violations

examples/demo_sphere.py

@@ -11,16 +11,17 @@ def sphere(x):
 
 
 def main():
+    dtype = np.float128 # np.float64 / np.float32
     dim = 40
-    mean = np.ones([dim, 1]) * 2
+    mean = np.full([dim, 1], 2.)
     sigma = 0.5
     lamb = 16 # note that lamb (sample size) should be even number
     iteration_number = 500
 
-    cr = CRFMNES(dim, sphere, mean, sigma, lamb)
+    cr = CRFMNES(dim, sphere, mean, sigma, lamb, dtype=dtype)
     x_best, f_best = cr.optimize(iteration_number)
     print("x_best:{}, f_best:{}".format(x_best, f_best))
 
 
 if __name__ == '__main__':
-    main()
+    main()

requirements.txt

@@ -1,2 +1,2 @@
-numpy==1.15.2
-pytest==3.9.1
+numpy>==1.24.1
+pytest==7.2.1

tests/test_constraint_sphere.py

@@ -21,7 +21,7 @@ def test_run_d40_const_sphere():
     allowable_evals = (19.4 + 1.1*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, const_sphere, mean, sigma, lamb)
+    cr = CRFMNES(dim, const_sphere, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12
@@ -36,7 +36,7 @@ def test_run_d80_const_sphere():
     allowable_evals = (48.8 + 1.4*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, const_sphere, mean, sigma, lamb)
+    cr = CRFMNES(dim, const_sphere, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12

tests/test_ellipsoid.py

@@ -20,7 +20,7 @@ def test_run_d40_ellipsoid():
     allowable_evals = (8.8 + 0.5*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, ellipsoid, mean, sigma, lamb)
+    cr = CRFMNES(dim, ellipsoid, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12
@@ -35,7 +35,7 @@ def test_run_d80_ellipsoid():
     allowable_evals = (17.5 + 0.6*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, ellipsoid, mean, sigma, lamb)
+    cr = CRFMNES(dim, ellipsoid, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12

tests/test_ktablet.py

@@ -21,7 +21,7 @@ def test_run_d40_ktablet():
     allowable_evals = (9.1 + 0.6*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, ktablet, mean, sigma, lamb)
+    cr = CRFMNES(dim, ktablet, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12
@@ -36,7 +36,7 @@ def test_run_d80_ktablet():
     allowable_evals = (18.6 + 0.8*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, ktablet, mean, sigma, lamb)
+    cr = CRFMNES(dim, ktablet, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12

tests/test_rastrigin.py

@@ -21,7 +21,7 @@ def test_run_d40_rastrigin():
     allowable_evals = (148 + 8.8*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, rastrigin, mean, sigma, lamb)
+    cr = CRFMNES(dim, rastrigin, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12
@@ -36,7 +36,7 @@ def test_run_d80_rastrigin():
     allowable_evals = (296 + 8.0*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, rastrigin, mean, sigma, lamb)
+    cr = CRFMNES(dim, rastrigin, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12

tests/test_rosenbrock_chain.py

@@ -18,7 +18,7 @@ def test_run_d40_rosen():
     allowable_evals = (52.2 + 1.5*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, rosenbrock_chain, mean, sigma, lamb)
+    cr = CRFMNES(dim, rosenbrock_chain, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12
@@ -33,7 +33,7 @@ def test_run_d80_rosen():
     allowable_evals = (192 + 2.0*3) * 1e3 # 2 sigma
     iteration_number = int(allowable_evals / lamb) + 1
 
-    cr = CRFMNES(dim, rosenbrock_chain, mean, sigma, lamb)
+    cr = CRFMNES(dim, rosenbrock_chain, mean, sigma, lamb, dtype=np.float128)
     x_best, f_best = cr.optimize(iteration_number)
     print("f_best:{}".format(f_best))
     assert f_best < 1e-12