pymodules/optimize.py at master · yaoyuki/pymodules · GitHub

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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
import numpy as np
import numba as nb
__epsilon__ = np.finfo(float).eps
@nb.jit(nopython = True)
# def brent(func, ax, cx, args = (), bx = None, tol=1.0e-8, itmax=500):
def brent(func, ax, bx, cx, args = (), tol=1.0e-8, itmax=500):

    #parameters
    CGOLD=0.3819660

    ZEPS=1.0e-3*__epsilon__
    #*np.finfo(float).eps

    brent = 1.0e20
    xmin = 1.0e20

    a=min(ax,cx)
    b=max(ax,cx)

#     bx = (a+b)/2.

    v=bx
    w=v
    x=v
    e=0.0
    # fx=func(x)
    fx=func(*((x, ) + args))
    fv=fx
    fw=fx

    d = 0.0

    it = 0
    for it in range(itmax):

        xm=0.5*(a+b)
        tol1=tol*abs(x)+ZEPS
        tol2=2.0*tol1
        if (abs(x-xm) <= (tol2-0.5*(b-a))):
            xmin=x
            brent=fx
            return xmin, brent

        if (abs(e) > tol1):
            r=(x-w)*(fx-fv)
            q=(x-v)*(fx-fw)
            p=(x-v)*q-(x-w)*r
            q=2.0*(q-r)

            if (q > 0.0):
                p=-p

            q=abs(q)
            etemp=e
            e=d

            if abs(p) >= abs(0.5*q*etemp) or  p <= q*(a-x) or p >= q*(b-x):

                #e=merge(a-x,b-x, x >= xm )
                if x >= xm:
                    e = a-x
                else:
                    e = b-x
                d=CGOLD*e

            else:
                d=p/q
                u=x+d

                if (u-a < tol2 or b-u < tol2):
                    d= abs(tol1)*np.sign(xm - x)  #sign(tol1,xm-x)

        else:

            if x >= xm:
                e = a-x
            else:
                e = b-x

            d=CGOLD*e

        u = 0.  #merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
        if abs(d) >= tol1:
            u = x+d
        else:
            u = x+abs(tol1)*np.sign(d)

        ###put your objective function also here###
        # fu = func(u)
        fu = func(*((u, ) + args))

        if (fu <= fx):
            if (u >= x):
                a=x
            else:
                b=x

            #shft(v,w,x,u)
            v = w
            w = x
            x = u
            #shft(fv,fw,fx,fu)
            fv = fw
            fw = fx
            fx = fu


        else:
            if (u < x):
                a=u
            else:
                b=u

            if fu <= fw or w == x:
                v=w
                fv=fw
                w=u
                fw=fu

            elif fu <= fv or v == x or v == w:
                v=u
                fv=fu

    if it == itmax-1:
        print('brent: exceed maximum iterations')

    return x, fx