SciML / NonlinearSolve.jl

High-performance and differentiation-enabled nonlinear solvers (Newton methods), bracketed rootfinding (bisection, Falsi), with sparsity and Newton-Krylov support.
https://docs.sciml.ai/NonlinearSolve/stable/
MIT License
226 stars 40 forks source link

Code restrictions for merit function not documented #298

Closed DavidSagan closed 9 months ago

DavidSagan commented 9 months ago

[Side note: You might want to consider a documentation problem category for issues.]

The following code fails:

using NonlinearSolve, LinearSolve, LinearAlgebra, Random, ForwardDiff, PyFormattedStrings

function fff(var, p)
  v_true = [1.0, 0.1, 2.0, 0.5]
  xx = [1.0, 2.0, 3.0, 4.0]
  println(f"{typeof(xx)}  {typeof(var)}  {typeof(v_true)}")
  xx[1] = var[1] - v_true[1]
  return var - v_true
end

v_true = [1.0, 0.1, 2.0, 0.5]
v_init = v_true .+ randn!(similar(v_true)) * 0.1

prob_oop = NonlinearLeastSquaresProblem{false}(fff, v_init)
sol = solve(prob_oop, LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)

The error generated is:

julia> include("z.jl")
Vector{Float64}  Vector{Float64}  Vector{Float64}
Vector{Float64}  Vector{Float64}  Vector{Float64}
Vector{Float64}  Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4}}  Vector{Float64}
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::VectorizationBase.Double{T}) where T<:Union{Float16, Float32, Float64, VectorizationBase.Vec{<:Any, <:Union{Float16, Float32, Float64}}, VectorizationBase.VecUnroll{var"#s45", var"#s44", var"#s43", V} where {var"#s45", var"#s44", var"#s43"<:Union{Float16, Float32, Float64}, V<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8, SIMDTypes.Bit, VectorizationBase.AbstractSIMD{var"#s44", var"#s43"}}}}
   @ VectorizationBase ~/.julia/packages/VectorizationBase/0dXyA/src/special/double.jl:111
  ...

Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4}, i1::Int64)
    @ Base ./array.jl:969
  [3] fff(var::Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4}}, p::SciMLBase.NullParameters)
    @ Main ~/Bmad/turn_by_turn_sextupole/z.jl:7
  [4] (::NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(fff), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 4}}, ::Vararg{Any})
    @ SciMLBase ~/.julia/packages/SciMLBase/ynHlA/src/scimlfunctions.jl:2408
... etc...

The error is generated in the line xx[1] = var[1] - v_true[1] in the function fff. I believe I understand this to be due to a call to fff inside solve with the argument var not being a vector but instead a construct using ForwardDiff. This is not covered in the documentation and I'm not sure how to get around this. [It is true that this line is completely irrelevant to the test case but in my real code I need to understand exactly what is happening.]

avik-pal commented 9 months ago

You will need https://github.com/SciML/PreallocationTools.jl

DavidSagan commented 9 months ago

@avik-pal Looking at the documentation for PreallocationTools it is not clear to me how to use this. Can you show me how to use this for implementing fff in the above example? That is, how to write xx[1] = var[1] - v_true[1] in a form that will not cause an error.

ChrisRackauckas commented 9 months ago

xx = eltype(var)[1.0, 2.0, 3.0, 4.0]

DavidSagan commented 9 months ago

@ChrisRackauckas Yes this works. Much Thanks!