-
Notifications
You must be signed in to change notification settings - Fork 0
/
torczon.py
100 lines (73 loc) · 2.35 KB
/
torczon.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
import numpy as np
# This is a simple implementation of a modified version of the torczon
# derivative free optimization algorithm.
def improve_S(S=[], device_f=[]):
# S is a matrix of dimension (n+1)x(n)
# each line is a argument values
# compute the best and make it the first
y = np.asarray([device_f(S[i, :]) for i in range(S.shape[0])])
i = np.argmin(y)
if i > 0:
S[[i, 0]] = S[[0, i]]
D = np.asarray([S[0, :]-S[i, :]
for i in range(S.shape[0])])
Sc = S+0.5*D
Sr = S+2*D
Se = S+4*D
yr = np.asarray([device_f(Sr[i, :]) for i in range(S.shape[0])]).min()
if yr < y[i]:
ye = np.asarray([device_f(Se[i, :]) for i in range(S.shape[0])]).min()
if ye < yr:
Snew = Se
ynew = ye
else:
Snew = Sr
ynew = yr
else:
Snew = Sc
ynew = y[i]
return Snew, y[i], ynew
def solve(f_user=[], par=[], x0=[], xmin=[],
xmax=[], Nguess=[], Niter=[], initial_box_width=0.1):
best_y = 1e19
n = len(xmin)
def f(x):
cond1 = np.all([x[i] <= xmax[i] for i in range(n)])
cond2 = np.all([x[i] >= xmin[i] for i in range(n)])
if cond1 & cond2:
return f_user(x, par)
else:
return 2e19
for ir in range(Nguess):
if (len(x0) > 0) and (ir == 0):
S = np.zeros((n+1, n))
S[0, :] = x0
for ix in range(n):
S[ix+1, :] = x0
S[ix+1, ix] = x0[ix]+initial_box_width
else:
Xmin = np.asarray([xmin for i in range(n+1)])
Xmax = np.asarray([xmax for i in range(n+1)])
S = Xmin + np.random.rand(n+1, n)*(Xmax-Xmin)
for i in range(Niter):
S, y, ynew = improve_S(S=S, device_f=f)
if ynew < best_y:
best_y = ynew
best_x = S[0, :]
class R:
pass
R.f, R.x = best_y, best_x
return R
if __name__ == '__main__':
class Data():
def __init__(self):
self.a = 10
self.b = -1
self.c = 0.0
p = Data()
def fdex(x, p):
return (x[0]-p.a)**2 +p.c*(x[1]*x[0])**4+(x[1]-p.b)**2
Nguess = 5
Niter = 30
R = solve(fdex, p, [0.0,0.0], [-10,-10], [10, 10], Nguess, Niter)
print(R.x)