Skip to content

The Nelder Mead optimisation method combined with winding numbers for root finding complex values functions

License

Notifications You must be signed in to change notification settings

jwscook/WindingNelderMead.jl

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

47 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

CI codecov DOI

WindingNelderMead.jl

Find the roots of a complex function, or even a function that returns two numbers converted into one that returns a complex number.

The following code

using WindingNelderMead, Plots
using WindingNelderMead: Vertex, Simplex, windingnumber, windingangle
using WindingNelderMead: centre, assessconvergence, position, value
using WindingNelderMead: bestvertex, issortedbyangle, hypervolume

function run()
  defaults = WindingNelderMead.convergenceconfig(2, Float64)
  positions = Vector{Float64}[]
  values = Vector{Float64}[]
  windingnumbers = Int[]
  initial_vertex_positions = [[0.0, 0.0], [0.1, 0.0], [0.0, 0.1]]
  objective(x) = (x[1] + im * x[2]) - (exp(1) + im * sqrt(2))
  simplex = Simplex(objective, initial_vertex_positions)
  simplices = [deepcopy(simplex)]
  function wrappedobjective(x)
    push!(positions, x)
    push!(windingnumbers, windingnumber(simplex))
    output = objective(x)
    push!(values, [reim(output)...])
    push!(simplices, simplex)
    return output
  end
  ics = rand(2)
  sizes = rand(2)
  solution = WindingNelderMead.optimise!(simplex, wrappedobjective)
  (s, n, returncode, its) = solution
  return positions, values, windingnumbers
end

positions, values, windingnumbers = run()

h = plot(layout=@layout [a b; c])
title!(h[1], "Input: complex x")
xlabel!(h[1], "Real x")
ylabel!(h[1], "Imag x")
title!(h[2], "Output: complex v")
xlabel!(h[2], "Real v")
ylabel!(h[2], "Imag v")
xlabel!(h[3], "Iteration")
ylabel!(h[3], "Log10 |v|")
anim = @animate for (i, w) in enumerate(windingnumbers) 
  x, y = positions[i]
  rv, iv = values[i]
  c = iszero(w) ? :blue : :red
  scatter!(h[1], [x], [y], mc=c, label=nothing,
           xlims=(-0.2, 4.2), ylims=(-0.2,2.2))
  scatter!(h[2], [rv], [iv], mc=c, label=nothing,
           xlims=(-3.2, 2.2), ylims=(-2.2,1.2))
  scatter!(h[3], [i], [log10(abs(rv + im * iv))], mc=c, label=nothing,
           xlims=(0, length(values) + 1), ylims=(-18,2))
  if i == 1
    scatter!(h[3], [-10], [0], mc=:blue, label="Nelder-Mead", legend=true)
    scatter!(h[3], [-10], [0], mc=:red, label="Winding", legend=true)
  else
    x_, y_ = positions[i-1]
    rv_, iv_ = values[i-1]
    plot!(h[1], [x_, x], [y_, y], lc=c, label=nothing)
    plot!(h[2], [rv_, rv], [iv_, iv], lc=c, label=nothing)
  end
end

gif(anim, "WindingNelderMead.gif", fps = 10)

results in

WindingNelderMead

Also, see the tests for more examples.

About

The Nelder Mead optimisation method combined with winding numbers for root finding complex values functions

Resources

License

Stars

Watchers

Forks

Packages

No packages published

Languages