Handle downstream usage issues

.github/workflows/Downstream.yml

@@ -18,6 +18,8 @@ jobs:
         os: [ubuntu-latest]
         package:
           - {user: SciML, repo: ModelingToolkit.jl, group: All}
+          - {user: SciML, repo: OrdinaryDiffEq.jl, group: Interface}
+          - {user: SciML, repo: OrdinaryDiffEq.jl, group: Regression}
     steps:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -1,27 +1,31 @@
 name = "NonlinearSolve"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-authors = ["Kanav Gupta <[email protected]>"]
+authors = ["SciML"]
 version = "0.3.24"
 
 [deps]
 ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 ArrayInterfaceCore = "0.1.1"
 FiniteDiff = "2"
 ForwardDiff = "0.10.3"
+LinearSolve = "1"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1"
-SciMLBase = "1.32"
+SciMLBase = "1.73"
 Setfield = "0.7, 0.8, 1"
 StaticArrays = "0.12,1.0"
 UnPack = "1.0"

docs/src/basics/FAQ.md

@@ -15,15 +15,15 @@ myfun(x, lv) = x * sin(x) - lv
 
 function f(out, levels, u0)
     for i in 1:N
-        out[i] = solve(NonlinearProblem{false}(NonlinearFunction{false}(myfun),
+        out[i] = solve(IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun),
                 u0, levels[i]), Falsi()).u
     end
 end
 
 function f2(out, levels, u0)
     for i in 1:N
-        out[i] = solve(NonlinearProblem{false}(NonlinearFunction{false}(myfun),
-                u0, levels[i]), NewtonRaphson()).u
+        out[i] = solve(IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun),
+                u0, levels[i]), SimpleNewtonRaphson()).u
     end
 end
 

docs/src/basics/NonlinearFunctions.md

@@ -9,5 +9,6 @@ of pre-computed functions to speed up the calculations. This is offered via the
 ## Function Type Definitions
 
 ```@docs
+SciMLBase.IntervalNonlinearFunction
 SciMLBase.NonlinearFunction
 ```

docs/src/basics/NonlinearProblem.md

@@ -1,5 +1,32 @@
 # Nonlinear Problems
 
+## The Two Types of Nonlinear Problems
+
+NonlinearSolve.jl tackles two related types of nonlinear systems:
+
+1. Interval rootfinding problems. I.e., find the ``t in [t_0, t_f]`` such that `f(t) = 0`.
+2. Systems of nonlinear equations, i.e. find the `u` such that `f(u) = 0`.
+
+The former is for solving scalar rootfinding problems, i.e. finding a single number, and
+requires that a bracketing interval is known. For a bracketing interval, one must have that
+the sign of `f(t_0)` is opposite the sign of `f(t_f)`, thus guaranteeing a root in the
+interval.
+
+The latter type of nonlinear system can be multidimensional and thus no ordering nor
+boundaries are assumed to be known. For a system of nonlinear equations, `f` can return
+an array and the solver seeks to find the value of `u` for which all outputs of `f` are
+simultaniously zero.
+
+!!! note
+
+    Interval rootfinding problems allow for `f` to return an array, in which case the interval
+    rootfinding problem is interpreted as finding the first `t` such that any of the components
+    of the array hit zero.
+
+
+## Problem Construction Details
+
 ```@docs
+SciMLBase.IntervalNonlinearProblem
 SciMLBase.NonlinearProblem
 ```

docs/src/index.md

@@ -12,8 +12,7 @@ allow for automatically generating high-performance code.
 
 Performance is key: the current methods are made to be highly performant on
 scalar and statically sized small problems, with options for large-scale systems.
-If you run into any performance issues, please file an issue. Note that this
-package is distinct from [SciMLNLSolve.jl](https://github.com/SciML/SciMLNLSolve.jl).
+If you run into any performance issues, please file an issue.
 Consult the [NonlinearSystemSolvers](@ref nonlinearsystemsolvers) page for
 information on how to import solvers from different packages.
 
@@ -79,7 +78,7 @@ Pkg.status(;mode = PKGMODE_MANIFEST) # hide
 </details>
 ```
 ```@raw html
-You can also download the 
+You can also download the
 <a href="
 ```
 ```@eval
@@ -100,4 +99,4 @@ link = "https://github.com/SciML/"*name*".jl/tree/gh-pages/v"*version*"/assets/P
 ```
 ```@raw html
 ">project</a> file.
-```
+```

docs/src/solvers/BracketingSolvers.md

@@ -1,8 +1,8 @@
-# Bracketing Solvers
+# Interval Rootfinding Methods (Bracketing Solvers)
 
-`solve(prob::NonlinearProblem,alg;kwargs)`
+`solve(prob::IntervalNonlinearProblem,alg;kwargs)`
 
-Solves for ``f(u)=0`` in the problem defined by `prob` using the algorithm
+Solves for ``f(t)=0`` in the problem defined by `prob` using the algorithm
 `alg`. If no algorithm is given, a default algorithm will be chosen.
 
 This page is solely focused on the bracketing methods for scalar nonlinear equations.
@@ -14,7 +14,10 @@ less stable than `Bisection`.
 
 ## Full List of Methods
 
-### NonlinearSolve.jl
+### SimpleNonlinearSolve.jl
+
+These methods are automatically included as part of NonlinearSolve.jl. Though one can use
+SimpleNonlinearSolve.jl directly to decrease the dependencies and improve load time.
 
 - `Falsi`: A non-allocating regula falsi method
 - `Bisection`: A common bisection method

docs/src/solvers/NonlinearSystemSolvers.md

@@ -9,13 +9,15 @@ This page is solely focused on the methods for nonlinear systems.
 
 ## Recommended Methods
 
-`NewtonRaphson` is a good choice for most problems. It is non-allocating on
-static arrays and thus really well-optimized for small systems, while for large
+`NewtonRaphson` is a good choice for most problems.  For large
 systems it can make use of sparsity patterns for sparse automatic differentiation
 and sparse linear solving of very large systems. That said, as a classic Newton
-method, its stability region can be smaller than other methods. `NLSolveJL`'s
+method, its stability region can be smaller than other methods. Meanwhile, `SimpleNewtonRaphson`
+is an implementation which is specialized for small equations. It is non-allocating on
+static arrays and thus really well-optimized for small systems, thus usually outperforming
+the other methods when such types are used for `u0`. `NLSolveJL`'s
 `:trust_region` method can be a good choice for high stability, along with
-`CMINPACK`.
+`CMINPACK`.s
 
 For a system which is very non-stiff (i.e., the condition number of the Jacobian
 is small, or the eigenvalues of the Jacobian are within a few orders of magnitude),
@@ -29,8 +31,17 @@ These are the core solvers.
 
 - `NewtonRaphson(;autodiff=true,chunk_size=12,diff_type=Val{:forward},linsolve=DEFAULT_LINSOLVE)`:
   A Newton-Raphson method with swappable nonlinear solvers and autodiff methods
-  for high performance on large and sparse systems. When used on objects like
-  static arrays, this method is non-allocating.
+  for high performance on large and sparse systems.
+
+### SimpleNonlinearSolve.jl
+
+These methods are included with NonlinearSolve.jl by default, though SimpleNonlinearSolve.jl
+can be used directly to reduce dependencies and improve load times.
+
+- `SimpleNewtonRaphson()`: A simplified implementation of the Newton-Raphson method. Has the
+  property that when used with states `u` as a `Number` or `StaticArray`, the solver is
+  very efficient and non-allocating. Thus this implmentation is well-suited for small
+  systems of equations.
 
 ### SciMLNLSolve.jl
 

docs/src/tutorials/iterator_interface.md

@@ -4,11 +4,8 @@ There is an iterator form of the nonlinear solver which mirrors the DiffEq integ
 
 ```julia
 f(u, p) = u .* u .- 2.0
-u0 = (1.0, 2.0) # brackets
+u0 = 1.5 # brackets
 probB = NonlinearProblem(f, u0)
-solver = init(probB, Falsi()) # Can iterate the solver object
-solver = solve!(solver)
+cache = init(probB, NewtonRaphson()) # Can iterate the solver object
+solver = solve!(cache)
 ```
-
-Note that the `solver` object is actually immutable since we want to make it
-live on the stack for the sake of performance.

docs/src/tutorials/nonlinear.md

@@ -27,6 +27,6 @@ a bracket instead of an initial condition, for example:
 ```julia
 f(u, p) = u .* u .- 2.0
 u0 = (1.0, 2.0) # brackets
-probB = NonlinearProblem(f, u0)
+probB = IntervalNonlinearProblem(f, u0)
 sol = solve(probB, Falsi())
 ```

src/NonlinearSolve.jl

@@ -9,25 +9,28 @@ using StaticArrays
 using RecursiveArrayTools
 using LinearAlgebra
 import ArrayInterfaceCore
+import LinearSolve
+using DiffEqBase
 
 @reexport using SciMLBase
+@reexport using SimpleNonlinearSolve
 
 abstract type AbstractNonlinearSolveAlgorithm <: SciMLBase.AbstractNonlinearAlgorithm end
-abstract type AbstractBracketingAlgorithm <: AbstractNonlinearSolveAlgorithm end
 abstract type AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ} <:
               AbstractNonlinearSolveAlgorithm end
-abstract type AbstractImmutableNonlinearSolver <: AbstractNonlinearSolveAlgorithm end
+
+function SciMLBase.__solve(prob::NonlinearProblem,
+                           alg::AbstractNonlinearSolveAlgorithm, args...;
+                           kwargs...)
+    cache = init(prob, alg, args...; kwargs...)
+    sol = solve!(cache)
+end
 
 include("utils.jl")
 include("jacobian.jl")
-include("types.jl")
-include("solve.jl")
-include("bisection.jl")
-include("falsi.jl")
 include("raphson.jl")
-include("scalar.jl")
+include("ad.jl")
 
-# DiffEq styled algorithms
-export Bisection, Falsi, NewtonRaphson
+export NewtonRaphson
 
 end # module

src/ad.jl

@@ -0,0 +1,39 @@
+function scalar_nlsolve_ad(prob, alg, args...; kwargs...)
+    f = prob.f
+    p = value(prob.p)
+
+    u0 = value(prob.u0)
+    newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
+
+    sol = solve(newprob, alg, args...; kwargs...)
+
+    uu = sol.u
+    if p isa Number
+        f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p)
+    else
+        f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p)
+    end
+
+    f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu)
+    pp = prob.p
+    sumfun = let f_x′ = -f_x
+        ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p)
+    end
+    partials = sum(sumfun, zip(f_p, pp))
+    return sol, partials
+end
+
+function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector}, iip,
+                                                <:Dual{T, V, P}}, alg::NewtonRaphson,
+                         args...; kwargs...) where {iip, T, V, P}
+    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
+    return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;
+                                    retcode = sol.retcode)
+end
+function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector}, iip,
+                                                <:AbstractArray{<:Dual{T, V, P}}},
+                         alg::NewtonRaphson, args...; kwargs...) where {iip, T, V, P}
+    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
+    return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;
+                                    retcode = sol.retcode)
+end