Use NaNMath

[anriseth]

In case something goes wrong during the interpolation step.

[codecov]

Codecov Report

Merging #71 into master will not change coverage.
The diff coverage is 100%.

@@          Coverage Diff           @@
##           master     #71   +/-   ##
======================================
  Coverage    61.9%   61.9%           
======================================
  Files           7       7           
  Lines         546     546           
======================================
  Hits          338     338           
  Misses        208     208

Impacted Files	Coverage Δ
src/backtracking.jl	94.28% <100%> (ø)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 675f5b5...ca92160. Read the comment docs.

[andreasnoack]

This causes the following problem

julia> NaNMath.min(ForwardDiff.Dual{Float64}(1), ForwardDiff.Dual{Float64}(1))
ERROR: StackOverflowError:
Stacktrace:
 [1] min(::ForwardDiff.Dual{Float64,Int64,0}, ::ForwardDiff.Dual{Float64,Int64,0}) at /Users/andreasnoack/.julia/packages/NaNMath/pEdac/src/NaNMath.jl:291 (repeats 80000 times)

julia> min(ForwardDiff.Dual{Float64}(1), ForwardDiff.Dual{Float64}(1))
Dual{Float64}(1)

which can happen with autodiff=:forward.

[pkofod]

@anriseth do you remember what part "could go wrong"?

[anriseth]

Hmm, no unfortunately not. I remember running into an issue where NaNs propagated through and the culprit may have been alphatmp in some weird edgecase.

We can either revert the Namath operations which are unlikely to affect as many people as the Forwarddiff problems, or add an explicit check on NaN instead
(I believe the problem was alphatmp, and not alpha, rhohi, or rholo)

[andreasnoack]

I tried to change NaNMath.min and max in Backtracking locally and I immediately got gradients full of NaNs so I probably hit the issue that made you do this PR.

[anriseth]

Ah, that's annoying. Do any of the other line search methods work for you?
I assume Static would work. We could make a vanilla backtracking as well, to get something, potentially, a little better than Static but which support Duals.

[pkofod]

The other searches don’t really work for Andreas, because they don’t have the initial finiteness backtracking step (it should just be included in the others as well)

[anriseth]

because they don’t have the initial finiteness backtracking step

I may have misunderstood, but I thought they did? :)

https://github.com/JuliaNLSolvers/LineSearches.jl/blob/master/src/hagerzhang.jl#L142
https://github.com/JuliaNLSolvers/LineSearches.jl/blob/master/src/morethuente.jl#L216
https://github.com/JuliaNLSolvers/LineSearches.jl/blob/master/src/static.jl#L20

[pkofod]

Oh that was added... did you try morethuente Andreas?

[andreasnoack]

I've locally defined DiffRules for NaNMath.min/max but then I hit the NaN issue right after. The culprit seems to be

LineSearches.jl/src/backtracking.jl

Line 32 in fae0768

α_0 = min(α_0, min(alphamax, ls.maxstep / norm(s, Inf)))

.
Since maxstep is initialized to Inf we get things like

julia> Inf/ForwardDiff.Dual{Float64}(1.0, 0.0)
Dual{Float64}(Inf,NaN)

Initially, I thought I could work around this by flipping the NANSAFE_MODE_ENABLED flag in ForwardDiff but it didn't work, see JuliaDiff/ForwardDiff.jl#179 (comment).

I guess it might be difficult to get correct behavior for Duals when we hit Infs. Maybe another solution could be to pick set ls.maxstep to a large but finite value.

[andreasnoack]

I can of course just set maxstep but it might be difficult to figure out.

[anriseth]

If I remember correctly maxstep is something @cortner found useful (correct me if I'm wrong). By setting the default to Inf we thought it wouldn't make any difference to people anyway, but you've proved that wrong :p

Two possible ways around this:

1. If you've defined DiffRules for NanMath.min, then would it work to make the line
   α_0 = NaNMath.min(α_0, NaNMath.min(alphamax, ls.maxstep / norm(s, Inf))) ?
2. We can explicitly check for NaN and skip the min operation when there are NaNs (or does that cause issues elsewhere?)

If you remove the line checking for maxstep locally, does the optimization run to convergence?

[cortner]

Correct - I introduced this As a stabilising mechanism.

Could one just check whether maxstep is Inf and in that case skip this operation?

[andreasnoack]

Regarding the initial comment then I think it should be fixed by JuliaDiff/DiffRules.jl#31.

Regarding the other issue then JuliaDiff/ForwardDiff.jl#386 might help here but it's still an option that is off by default.

I'm not sure what the better solution is but at least I've been able to get things working now for my use case by specifying a finite maxstep.