SHARE
TWEET

Untitled

a guest Sep 22nd, 2019 90 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. export nelder_mead
  2. """
  3. nelder_mead(nlp)
  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. """
  7. function nelder_mead(nlp :: AbstractNLPModel;
  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. OK, I Understand
 
Top