• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Sep 22nd, 2019 90 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
2. """
4. This implementation follows the algorithm described in chapter 9 of [1].
5. [1] Numerical Optimization (Jorge Nocedal and Stephen J. Wright), Springer, 2006.
6. """
8.                     x :: AbstractVector = copy(nlp.meta.x0),
9.                     vertices :: AbstractArray = Array{typeof(x)}(undef, nlp.meta.nvar+1),
10.             tol :: Real = √eps(eltype(x)),
11.                     max_time :: Float64 = 30.0,
12.                     max_eval :: Int = -1)
13.
14.   # Initial simplex
15.   n = nlp.meta.nvar
16.   # If the user doesn't provide any vertices points, should it be initialized with undef?
17.   # Another alternative is to make sure the user provides a matrix instead of an Array
18.   if !(isassigned(vertices))
19.   vertices[1] = x
20.     for j in 1:n
21.       xt = copy(x)
22.       xt[j] += xt[j] == 0 ? 0.00025 : 0.05
23.       vertices[j+1] = xt
24.     end
25.   end
26.
27.   # Really unsure about both these lines, looks concise, but feels wrong
28.   # Since we're storing an array of different types, the data resulting datatype is Any
29.   # Another alternative is to instead of using a single array, use two, one for the points, another for the obj value.
30.   vertices = hcat(vertices, map(x->obj(nlp, x), vertices))
31.   vertices = vertices[sortperm(vertices[:, 2]), :]
32.   println("Vertices Datatype : " typeof(vertices))
33.
34.   k = 0
35.   el_time = 0.0
36.   start_time = time()
37.   tired = neval_obj(nlp) > max_eval >= 0 || el_time > max_time
38.   norm_opt = norm(vertices[1, 2] - vertices[n+1, 2])
39.   optimal =  norm_opt < tol
40.   T = eltype(x)
41.   status =:unknown
42.   @info log_header([:iter, :f, :nrm], [Int, T, T], hdr_override=Dict(:f => "f(x)", :nrm => "‖x₁ - xₙ+1‖"))
43.
44.   while !(optimal || tired)
45.
46.     shrink = true
47.     x_cen = sum(vertices[i, 1] for i in 1:n)/n
48.     μ(t) = x_cen + t * (vertices[n+1, 1] - x_cen)
49.     x_ref = μ(-1)
50.     f_ref = obj(nlp, x_ref)
51.     f_bver = vertices[1, 2]
52.
53.     @info log_row(Any[k, f_bver, norm_opt])
54.
55.     # Reflection
56.     if vertices[1,2] <= f_ref < vertices[n,2]
57.       vertices[n+1,1] = x_ref
58.       shrink = false
59.     # Expansion
60.     elseif f_ref < vertices[1, 2]
61.       x_exp = μ(-2)
62.       f_exp = obj(nlp,x_exp)
63.       f_exp < f_ref ? vertices[n+1, 1] = x_exp : vertices[n+1, 1] = x_ref
64.       shrink = false
65.     # Contraction
66.     elseif f_ref >= vertices[n, 2]
67.       # Outside
68.       if vertices[n, 2] <= f_ref < vertices[n+1, 2]
69.         x_ocn = μ(-0.5)
70.         f_ocn = obj(nlp, x_ocn)
71.         if f_ocn <= f_ref
72.           vertices[n+1, 1] = x_ocn
73.           shrink = false
74.         end
75.       # Inside
76.     else
77.         x_icn = μ(0.5)
78.         f_icn = obj(nlp, x_icn)
79.         if f_icn < vertices[n+1, 2]
80.           vertices[n+1, 1]= x_icn
81.           shrink = false
82.         end
83.       end
84.     end
85.
86.     if shrink
87.       for i in 2:n+1
88.         vertices[i, 1] = 0.5 * (vertices[1, 1] + vertices[i, 1])
89.       end
90.     end
91.
92.     vertices[1:n+1, 2] = map(x->obj(nlp, x), vertices[1:n+1, 1])
93.     vertices = vertices[sortperm(vertices[:, 2]), :]
94.     k += 1
95.     x = vertices[1, 1]
96.     tired = neval_obj(nlp) > max_eval >= 0 || el_time > max_time
97.     norm_opt = norm(vertices[1, 2] - vertices[n+1, 2])
98.     optimal =  norm_opt < tol
99.     el_time = time() - start_time
100.
101.     # As the algorithm goes, Float32 becomes Float64 and the same happens with Float16.
102.     # No clue why
103.     println("Vertices eltype : " eltype(vertices[1]))
104.     for i in 1:nlp.meta.nvar+1
105.       println("Vertice " ,i , " : " vertices[i,:])
106.     end
107.
108.   end
109.
110.   if optimal
111.     status = :acceptable
112.   elseif tired
113.     if neval_obj(nlp) > max_eval ≥  0
114.       status = :max_eval
115.     elseif el_time >= max_time
116.       status = :max_time
117.     end
118.   end
119.
120.   return GenericExecutionStats(status, nlp, solution=vertices[1, 1], objective=vertices[1, 2],
121.                               iter = k, elapsed_time = el_time)
122. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top