1import time 2 3import numpy as np 4from numpy import sin, cos, pi, exp, sqrt, abs 5from scipy.optimize import rosen 6 7 8class SimpleQuadratic: 9 10 def fun(self, x): 11 return np.dot(x, x) 12 13 def der(self, x): 14 return 2. * x 15 16 def hess(self, x): 17 return 2. * np.eye(x.size) 18 19 20class AsymmetricQuadratic: 21 22 def fun(self, x): 23 return np.dot(x, x) + x[0] 24 25 def der(self, x): 26 d = 2. * x 27 d[0] += 1 28 return d 29 30 def hess(self, x): 31 return 2. * np.eye(x.size) 32 33 34class SlowRosen: 35 36 def fun(self, x): 37 time.sleep(40e-6) 38 return rosen(x) 39 40 41class LJ: 42 """ 43 The Lennard Jones potential 44 45 a mathematically simple model that approximates the interaction between a 46 pair of neutral atoms or molecules. 47 https://en.wikipedia.org/wiki/Lennard-Jones_potential 48 49 E = sum_ij V(r_ij) 50 51 where r_ij is the cartesian distance between atom i and atom j, and the 52 pair potential has the form 53 54 V(r) = 4 * eps * ( (sigma / r)**12 - (sigma / r)**6 55 56 Notes 57 ----- 58 the double loop over many atoms makes this *very* slow in Python. If it 59 were in a compiled language it would be much faster. 60 """ 61 62 def __init__(self, eps=1.0, sig=1.0): 63 self.sig = sig 64 self.eps = eps 65 66 def vij(self, r): 67 return 4. * self.eps * ((self.sig / r)**12 - (self.sig / r)**6) 68 69 def dvij(self, r): 70 p7 = 6. / self.sig * (self.sig / r)**7 71 p13 = -12. / self.sig * (self.sig / r)**13 72 return 4. * self.eps * (p7 + p13) 73 74 def fun(self, coords): 75 natoms = coords.size // 3 76 coords = np.reshape(coords, [natoms, 3]) 77 energy = 0. 78 for i in range(natoms): 79 for j in range(i + 1, natoms): 80 dr = coords[j, :] - coords[i, :] 81 r = np.linalg.norm(dr) 82 energy += self.vij(r) 83 return energy 84 85 def der(self, coords): 86 natoms = coords.size // 3 87 coords = np.reshape(coords, [natoms, 3]) 88 energy = 0. 89 grad = np.zeros([natoms, 3]) 90 for i in range(natoms): 91 for j in range(i + 1, natoms): 92 dr = coords[j, :] - coords[i, :] 93 r = np.linalg.norm(dr) 94 energy += self.vij(r) 95 g = self.dvij(r) 96 grad[i, :] += -g * dr/r 97 grad[j, :] += g * dr/r 98 grad = grad.reshape([natoms * 3]) 99 return grad 100 101 def get_random_configuration(self): 102 rnd = np.random.uniform(-1, 1, [3 * self.natoms]) 103 return rnd * float(self.natoms)**(1. / 3) 104 105 106class LJ38(LJ): 107 natoms = 38 108 target_E = -173.928427 109 110 111class LJ30(LJ): 112 natoms = 30 113 target_E = -128.286571 114 115 116class LJ20(LJ): 117 natoms = 20 118 target_E = -77.177043 119 120 121class LJ13(LJ): 122 natoms = 13 123 target_E = -44.326801 124 125 126class Booth: 127 target_E = 0. 128 solution = np.array([1., 3.]) 129 xmin = np.array([-10., -10.]) 130 xmax = np.array([10., 10.]) 131 132 def fun(self, coords): 133 x, y = coords 134 return (x + 2. * y - 7.)**2 + (2. * x + y - 5.)**2 135 136 def der(self, coords): 137 x, y = coords 138 dfdx = 2. * (x + 2. * y - 7.) + 4. * (2. * x + y - 5.) 139 dfdy = 4. * (x + 2. * y - 7.) + 2. * (2. * x + y - 5.) 140 return np.array([dfdx, dfdy]) 141 142 143class Beale: 144 target_E = 0. 145 solution = np.array([3., 0.5]) 146 xmin = np.array([-4.5, -4.5]) 147 xmax = np.array([4.5, 4.5]) 148 149 def fun(self, coords): 150 x, y = coords 151 p1 = (1.5 - x + x * y)**2 152 p2 = (2.25 - x + x * y**2)**2 153 p3 = (2.625 - x + x * y**3)**2 154 return p1 + p2 + p3 155 156 def der(self, coords): 157 x, y = coords 158 dfdx = (2. * (1.5 - x + x * y) * (-1. + y) + 159 2. * (2.25 - x + x * y**2) * (-1. + y**2) + 160 2. * (2.625 - x + x * y**3) * (-1. + y**3)) 161 dfdy = (2. * (1.5 - x + x * y) * (x) + 162 2. * (2.25 - x + x * y**2) * (2. * y * x) + 163 2. * (2.625 - x + x * y**3) * (3. * x * y**2)) 164 return np.array([dfdx, dfdy]) 165 166 167""" 168Global Test functions for minimizers. 169 170HolderTable, Ackey and Levi have many competing local minima and are suited 171for global minimizers such as basinhopping or differential_evolution. 172(https://en.wikipedia.org/wiki/Test_functions_for_optimization) 173 174See also https://mpra.ub.uni-muenchen.de/2718/1/MPRA_paper_2718.pdf 175 176""" 177 178 179class HolderTable: 180 target_E = -19.2085 181 solution = [8.05502, 9.66459] 182 xmin = np.array([-10, -10]) 183 xmax = np.array([10, 10]) 184 stepsize = 2. 185 temperature = 2. 186 187 def fun(self, x): 188 return - abs(sin(x[0]) * cos(x[1]) * exp(abs(1. - sqrt(x[0]**2 + 189 x[1]**2) / pi))) 190 191 def dabs(self, x): 192 """derivative of absolute value""" 193 if x < 0: 194 return -1. 195 elif x > 0: 196 return 1. 197 else: 198 return 0. 199 200#commented out at the because it causes FloatingPointError in 201#basinhopping 202# def der(self, x): 203# R = sqrt(x[0]**2 + x[1]**2) 204# g = 1. - R / pi 205# f = sin(x[0]) * cos(x[1]) * exp(abs(g)) 206# E = -abs(f) 207# 208# dRdx = x[0] / R 209# dgdx = - dRdx / pi 210# dfdx = cos(x[0]) * cos(x[1]) * exp(abs(g)) + f * self.dabs(g) * dgdx 211# dEdx = - self.dabs(f) * dfdx 212# 213# dRdy = x[1] / R 214# dgdy = - dRdy / pi 215# dfdy = -sin(x[0]) * sin(x[1]) * exp(abs(g)) + f * self.dabs(g) * dgdy 216# dEdy = - self.dabs(f) * dfdy 217# return np.array([dEdx, dEdy]) 218 219 220class Ackley: 221 # note: this function is not smooth at the origin. the gradient will never 222 # converge in the minimizer 223 target_E = 0. 224 solution = [0., 0.] 225 xmin = np.array([-5, -5]) 226 xmax = np.array([5, 5]) 227 228 def fun(self, x): 229 E = (-20. * exp(-0.2 * sqrt(0.5 * (x[0]**2 + x[1]**2))) + 20. + np.e - 230 exp(0.5 * (cos(2. * pi * x[0]) + cos(2. * pi * x[1])))) 231 return E 232 233 def der(self, x): 234 R = sqrt(x[0]**2 + x[1]**2) 235 term1 = -20. * exp(-0.2 * R) 236 term2 = -exp(0.5 * (cos(2. * pi * x[0]) + cos(2. * pi * x[1]))) 237 238 deriv1 = term1 * (-0.2 * 0.5 / R) 239 240 dfdx = 2. * deriv1 * x[0] - term2 * pi * sin(2. * pi * x[0]) 241 dfdy = 2. * deriv1 * x[1] - term2 * pi * sin(2. * pi * x[1]) 242 243 return np.array([dfdx, dfdy]) 244 245 246class Levi: 247 target_E = 0. 248 solution = [1., 1.] 249 xmin = np.array([-10, -10]) 250 xmax = np.array([10, 10]) 251 252 def fun(self, x): 253 E = (sin(3. * pi * x[0])**2 + (x[0] - 1.)**2 * 254 (1. + sin(3 * pi * x[1])**2) + 255 (x[1] - 1.)**2 * (1. + sin(2 * pi * x[1])**2)) 256 return E 257 258 def der(self, x): 259 260 dfdx = (2. * 3. * pi * 261 cos(3. * pi * x[0]) * sin(3. * pi * x[0]) + 262 2. * (x[0] - 1.) * (1. + sin(3 * pi * x[1])**2)) 263 264 dfdy = ((x[0] - 1.)**2 * 2. * 3. * pi * cos(3. * pi * x[1]) * sin(3. * 265 pi * x[1]) + 2. * (x[1] - 1.) * 266 (1. + sin(2 * pi * x[1])**2) + (x[1] - 1.)**2 * 267 2. * 2. * pi * cos(2. * pi * x[1]) * sin(2. * pi * x[1])) 268 269 return np.array([dfdx, dfdy]) 270 271 272class EggHolder: 273 target_E = -959.6407 274 solution = [512, 404.2319] 275 xmin = np.array([-512., -512]) 276 xmax = np.array([512., 512]) 277 278 def fun(self, x): 279 a = -(x[1] + 47) * np.sin(np.sqrt(abs(x[1] + x[0]/2. + 47))) 280 b = -x[0] * np.sin(np.sqrt(abs(x[0] - (x[1] + 47)))) 281 return a + b 282 283 284class CrossInTray: 285 target_E = -2.06261 286 solution = [1.34941, -1.34941] 287 xmin = np.array([-10., -10]) 288 xmax = np.array([10., 10]) 289 290 def fun(self, x): 291 arg = abs(100 - sqrt(x[0]**2 + x[1]**2)/pi) 292 val = np.power(abs(sin(x[0]) * sin(x[1]) * exp(arg)) + 1., 0.1) 293 return -0.0001 * val 294 295 296class Schaffer2: 297 target_E = 0 298 solution = [0., 0.] 299 xmin = np.array([-100., -100]) 300 xmax = np.array([100., 100]) 301 302 def fun(self, x): 303 num = np.power(np.sin(x[0]**2 - x[1]**2), 2) - 0.5 304 den = np.power(1 + 0.001 * (x[0]**2 + x[1]**2), 2) 305 return 0.5 + num / den 306 307 308class Schaffer4: 309 target_E = 0.292579 310 solution = [0, 1.253131828927371] 311 xmin = np.array([-100., -100]) 312 xmax = np.array([100., 100]) 313 314 def fun(self, x): 315 num = cos(sin(abs(x[0]**2 - x[1]**2)))**2 - 0.5 316 den = (1+0.001*(x[0]**2 + x[1]**2))**2 317 return 0.5 + num / den 318