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