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
229 stars 41 forks source link

Remake defined but doesn't seem to do anything #168

Closed aabills closed 8 months ago

aabills commented 1 year ago

MWE:

using NonlinearSolve

function mwe(out, u, p)
    out[1] = u[1] - p[1]
end

p = [5.0]

prob = NonlinearProblem(mwe, [5.2], p=p)
u = solve(prob, NewtonRaphson())

prob2 = remake(prob, p=[6.0])
u2 = solve(prob2, NewtonRaphson())

Results in:

julia> u
u: 1-element Vector{Float64}:
 5.0

julia> u2
u: 1-element Vector{Float64}:
 5.0

But when looking at ?remake:

 remake(prob::NonlinearProblem; f = missing, u0 = missing, p = missing,
      problem_type = missing, kwargs = missing, _kwargs...)

  Remake the given NonlinearProblem. If u0 or p are given as symbolic maps ModelingToolkit.jl has to be loaded.

I don't know if this is expected behavior (and maybe should result in an error or a warning) or should just be fixed, but I felt like it should be reported as it can result in a silent wrong answer.

DanielVandH commented 1 year ago

The problem is that p is not a keyword argument for NonlinearProblem, so passing p as you are doing it gets pushed into prob.kwargs instead:

julia> prob = NonlinearProblem(mwe, [5.2], p=p)
NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 1-element Vector{Float64}:
 5.2

julia> prob.p
SciMLBase.NullParameters()

julia> prob.kwargs
pairs(::NamedTuple) with 1 entry:
  :p => [5.0]

It works if you don't do this -

using NonlinearSolve

function mwe(out, u, p)
    out[1] = u[1] - p[1]
end

p = [5.0]

prob = NonlinearProblem(mwe, [5.2], p)
prob.p 
prob.kwargs 

u = solve(prob, NewtonRaphson())

prob2 = remake(prob, p=[6.0])
u2 = solve(prob2, NewtonRaphson())
julia> u
u: 1-element Vector{Float64}:
 5.0

julia> u2
u: 1-element Vector{Float64}:
 6.0
ChrisRackauckas commented 1 year ago

We should probably catch and throw an error if someone does p = ....