But it’s been abandoned for a long time.
Unless you are interested in algorithm development you should use a mainstream code- there are plenty of them.
2 Likes
I’m the other contributor, and I approve of this message 
1 Like
Hmmm…
I’m having an exact the same problem as OP: a transition state on a 2D surface.
I played around with Optim, but it failed to give a correct result, unless my initial guess is super good.
At the end, I wrote a simple piece of code implementing Newton’s method. The initial guess was from a linear interpolation of the two minima. It worked surprisingly, but I don’t know if it’s robust at all.
Maybe I should try NEB.
Thanks for sugeeting IntervalRootFinding. It should be able to solve my ccurent problem, but in the future each dimension might be extend to N dimension (2 N-D arrays), so I guess it won’t work any more.
What is MWE?
I meant a minimal example of a problem where you want to find a saddle point.
What do you mean by 2 N-D arrays?
If the problem has too many dimensions then indeed IntervalRootFinding.jl will not be able to solve it.
As I noted in another thread, given a local optimum computed by any method you can always find the Lagrange multipliers afterwards just by solving a small linear system.
2 Likes
I see.
In fact I just tested IntervalRootFinding. It works fine for a function, but not an object from Dierckx.Spline1D, right? Part of my target function is from an interpolation of a complicated surface. I could fit them and code the polynomials in julia, but it isn’t desirable for some reasons.
using DelimitedFiles
using IntervalArithmetic, IntervalRootFinding
using Dierckx
potential = readdlm("pes.txt")
dipole = readdlm("dm.txt")
pes = Spline1D(potential[:, 1], potential[:, 2])
dm = Spline1D(dipole[:, 1], dipole[:, 2])
pho(q) = 0.5 * q^2
inter(r, q) = dm(r) * q + dm(q)^2
utotal((r, q)) = pes(r) + inter(r, q) + pho(q)
∇u = ∇(utotal)
println(roots(∇u, IntervalBox(-20..20, 2), Newton, 1e-5)
With pes.txt being
-10 0.192779
-7.9798 0.0408572
-5.9596 0.000340054
-3.93939 0.0125455
-1.91919 0.0358005
0.10101 0.0454414
2.12121 0.033814
4.14141 0.0102734
6.16162 0.0011843
8.18182 0.0499205
10 0.192779
and dm.txt
-10 0.578831
-7.9798 2.15586
-5.9596 2.85126
-3.93939 2.59403
-1.91919 1.49945
0.10101 -0.0830599
2.12121 -1.63838
4.14141 -2.66116
6.16162 -2.82458
8.18182 -2.03555
10 -0.578831
Of course, you can feed any 2D array to Spline1D for a testing purpose.
The code above returns
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(uTotal),Interval{Float64}},Interval{Float64},2})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
Float64(::T) where T<:Number at boot.jl:715
Float64(::Int8) at float.jl:60
...
A normal 2-variable function like f((x, y)) = sin(x) + cos(y) works flawlessly, but even the univariate pes function doesn’t work in roots. I guess it’s the poblem with a spline’d function.
It’s basically just 2N dimensional, with each of the original dimension being described as N variable coupled.
1 Like
Yes, due to how it works internally (evaluating functions on intervals), IntervalRootFinding.jl will only work with functions that are written in Julia, whereas Dierckx is a wrapper around non-Julia code.
It may (should) work if you replace the interpolation with a Julia version.
Hmm… I thought Dierckx was a Julia re-implementation of the original Fortran code.
So Interpolations should work?
The Dierckx.jl README says:
This is a Julia wrapper for the dierckx Fortran library, the same library underlying the spline classes in scipy.interpolate. Some of the functionality here overlaps with Interpolations.jl, a pure-Julia interpolation package. Take a look at it if you have a use case not covered here."
I believe that Interpolations.jl should work, but I have never used the package so I can’t try it out.
I think I got the same problem.
Just did this simple test.
using Interpolations
using IntervalArithmetic, IntervalRootFinding
x = collect(-10:1.0:10.0)
y = x.^2
itp = LinearInterpolation(x, y)
roots(itp, -1.0..1.0)
Error message
ERROR: MethodError: no method matching roots(::Interpolations.Extrapolation{Float64,1,Interpolations.GriddedInterpolation{Float64,1,Float64,Gridded{Linear},Tuple{Array{Float64,1}}},Gridded{Linear},Throw{Nothing}}, ::Interval{Float64})
That’s not the same problem; this one should be fixable!
Basically roots is expecting a function I guess, not a LinearInterpolation object.
Can you turn it into a function?
No, unfortunately it’s more complicated than I thought:
julia> roots(x -> itp(x), -1.0..1.0)
ERROR: TypeError: non-boolean (Interval{Float64}) used in boolean context
Stacktrace:
[1] signless(::Float64, ::Interval{Float64}) at ./operators.jl:126
I believe that roots is trying to evaluate the Interpolations object over an interval, which is failing. Fixing that won’t be easy unfortunately.
I think it wouldn’t be too hard to get a valid interval version of a linear interpolation since you just need to evaluate at the endpoints and any knot that lies inside the interval and then take a convex hull. So just a searchsorted on the knots field: Interpolations.jl/gridded.jl at bcd05a3f0843661104589c31da8d257fecdbe265 · JuliaMath/Interpolations.jl · GitHub.
Yes agreed, but I expect the OP really wants to use something like cubic splines, and in higher dimensions.
This should also be doable using the piecewise approach I outlined here: Latest recommendations for global optimization - #45 by dpsanders
1 Like