From 75c0fb610f9b93d2070290d3859e7d8a25843f6a Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 5 Mar 2024 16:04:24 -0500 Subject: [PATCH 1/5] Export wrappers over LineSearches --- Project.toml | 2 +- src/NonlinearSolve.jl | 8 +++++--- src/globalization/line_search.jl | 16 ++++++++++++++++ 3 files changed, 22 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index bc1919c29..1294701b1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NonlinearSolve" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" authors = ["SciML"] -version = "3.7.3" +version = "3.8.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 991783672..e7fa9bd3e 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -9,8 +9,8 @@ import PrecompileTools: @recompile_invalidations, @compile_workload, @setup_work @recompile_invalidations begin using ADTypes, ConcreteStructs, DiffEqBase, FastBroadcast, FastClosures, LazyArrays, - LineSearches, LinearAlgebra, LinearSolve, MaybeInplace, Preferences, Printf, - SciMLBase, SimpleNonlinearSolve, SparseArrays, SparseDiffTools + LinearAlgebra, LinearSolve, MaybeInplace, Preferences, Printf, SciMLBase, + SimpleNonlinearSolve, SparseArrays, SparseDiffTools import ArrayInterface: undefmatrix, can_setindex, restructure, fast_scalar_indexing import DiffEqBase: AbstractNonlinearTerminationMode, @@ -20,6 +20,7 @@ import PrecompileTools: @recompile_invalidations, @compile_workload, @setup_work import FiniteDiff import ForwardDiff import ForwardDiff: Dual + import LineSearches import LinearSolve: ComposePreconditioner, InvPreconditioner, needs_concrete_A import RecursiveArrayTools: recursivecopy!, recursivefill! @@ -29,7 +30,7 @@ import PrecompileTools: @recompile_invalidations, @compile_workload, @setup_work import StaticArraysCore: StaticArray, SVector, SArray, MArray, Size, SMatrix, MMatrix end -@reexport using ADTypes, LineSearches, SciMLBase, SimpleNonlinearSolve +@reexport using ADTypes, SciMLBase, SimpleNonlinearSolve const AbstractSparseADType = Union{ADTypes.AbstractSparseFiniteDifferences, ADTypes.AbstractSparseForwardMode, ADTypes.AbstractSparseReverseMode} @@ -157,6 +158,7 @@ export NewtonDescent, SteepestDescent, Dogleg, DampedNewtonDescent, GeodesicAcce # Globalization ## Line Search Algorithms export LineSearchesJL, NoLineSearch, RobustNonMonotoneLineSearch, LiFukushimaLineSearch +export Static, HagerZhang, MoreThuente, StrongWolfe, BackTracking ## Trust Region Algorithms export RadiusUpdateSchemes diff --git a/src/globalization/line_search.jl b/src/globalization/line_search.jl index 9d4f52e71..56bb98cf6 100644 --- a/src/globalization/line_search.jl +++ b/src/globalization/line_search.jl @@ -53,6 +53,10 @@ end LineSearchesJL(method; kwargs...) = LineSearchesJL(; method, kwargs...) function LineSearchesJL(; method = LineSearches.Static(), autodiff = nothing, α = true) + if method isa LineSearchesJL # Prevent breaking old code + return LineSearchesJL(method.method, α, autodiff) + end + if method isa AbstractNonlinearSolveLineSearchAlgorithm Base.depwarn("Passing a native NonlinearSolve line search algorithm to \ `LineSearchesJL` or `LineSearch` is deprecated. Pass the method \ @@ -65,6 +69,18 @@ end Base.@deprecate_binding LineSearch LineSearchesJL true +Static(args...; kwargs...) = LineSearchesJL(LineSearches.Static(args...; kwargs...)) +HagerZhang(args...; kwargs...) = LineSearchesJL(LineSearches.HagerZhang(args...; kwargs...)) +function MoreThuente(args...; kwargs...) + return LineSearchesJL(LineSearches.MoreThuente(args...; kwargs...)) +end +function BackTracking(args...; kwargs...) + return LineSearchesJL(LineSearches.BackTracking(args...; kwargs...)) +end +function StrongWolfe(args...; kwargs...) + return LineSearchesJL(LineSearches.StrongWolfe(args...; kwargs...)) +end + # Wrapper over LineSearches.jl algorithms @concrete mutable struct LineSearchesJLCache <: AbstractNonlinearSolveLineSearchCache f From 59118352118cf08720533d534366293fc092dd08 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 5 Mar 2024 16:27:53 -0500 Subject: [PATCH 2/5] All polyalg in 23 test problems --- src/algorithms/lbroyden.jl | 4 ++-- test/core/23_test_problems_tests.jl | 12 +++++++++++- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/algorithms/lbroyden.jl b/src/algorithms/lbroyden.jl index 499cc7c34..f1959e7d0 100644 --- a/src/algorithms/lbroyden.jl +++ b/src/algorithms/lbroyden.jl @@ -159,8 +159,8 @@ function LinearAlgebra.mul!(y::AbstractVector, x::AbstractVector, J::BroydenLowR return y end -function LinearAlgebra.mul!( - J::BroydenLowRankJacobian, u, vᵀ::LinearAlgebra.AdjOrTransAbsVec, α::Bool, β::Bool) +function LinearAlgebra.mul!(J::BroydenLowRankJacobian, u::AbstractArray, + vᵀ::LinearAlgebra.AdjOrTransAbsVec, α::Bool, β::Bool) @assert α & β idx_update = mod1(J.idx + 1, size(J.U, 2)) copyto!(@view(J.U[:, idx_update]), _vec(u)) diff --git a/test/core/23_test_problems_tests.jl b/test/core/23_test_problems_tests.jl index 7924e372d..b0f17113a 100644 --- a/test/core/23_test_problems_tests.jl +++ b/test/core/23_test_problems_tests.jl @@ -40,6 +40,16 @@ end export test_on_library, problems, dicts end +@testitem "PolyAlgorithms" setup=[RobustnessTesting] begin + alg_ops = (RobustMultiNewton(), FastShortcutNonlinearPolyalg()) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [] + broken_tests[alg_ops[2]] = [] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + @testitem "NewtonRaphson" setup=[RobustnessTesting] begin alg_ops = (NewtonRaphson(),) @@ -91,7 +101,7 @@ end test_on_library(problems, dicts, alg_ops, broken_tests) end -@testitem "Broyden" retries=5 setup=[RobustnessTesting] begin +@testitem "Broyden" setup=[RobustnessTesting] begin alg_ops = (Broyden(), Broyden(; init_jacobian = Val(:true_jacobian)), Broyden(; update_rule = Val(:bad_broyden)), Broyden(; init_jacobian = Val(:true_jacobian), update_rule = Val(:bad_broyden))) From 6f68b2aac3bbc00669c6c15fc9b9b2f08794ba5d Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 5 Mar 2024 17:32:11 -0500 Subject: [PATCH 3/5] Use start index --- src/algorithms/klement.jl | 4 +- src/core/approximate_jacobian.jl | 4 +- src/default.jl | 68 +++++++++++++++++++++--------- test/core/forward_ad_tests.jl | 4 +- test/core/rootfind_tests.jl | 8 ++-- test/misc/matrix_resizing_tests.jl | 8 ++-- 6 files changed, 62 insertions(+), 34 deletions(-) diff --git a/src/algorithms/klement.jl b/src/algorithms/klement.jl index 428dc382c..66daec1b7 100644 --- a/src/algorithms/klement.jl +++ b/src/algorithms/klement.jl @@ -25,8 +25,8 @@ over this. differentiable problems. """ function Klement(; max_resets::Int = 100, linsolve = nothing, alpha = nothing, - linesearch = NoLineSearch(), precs = DEFAULT_PRECS, autodiff = nothing, - init_jacobian::Val{IJ} = Val(:identity)) where {IJ} + linesearch = NoLineSearch(), precs = DEFAULT_PRECS, + autodiff = nothing, init_jacobian::Val{IJ} = Val(:identity)) where {IJ} if !(linesearch isa AbstractNonlinearSolveLineSearchAlgorithm) Base.depwarn( "Passing in a `LineSearches.jl` algorithm directly is deprecated. \ diff --git a/src/core/approximate_jacobian.jl b/src/core/approximate_jacobian.jl index d0471e99b..4204b2db7 100644 --- a/src/core/approximate_jacobian.jl +++ b/src/core/approximate_jacobian.jl @@ -66,8 +66,8 @@ function ApproximateJacobianSolveAlgorithm{concrete_jac, name}(; linesearch = LineSearchesJL(; method = linesearch) end return ApproximateJacobianSolveAlgorithm{concrete_jac, name}( - linesearch, trustregion, descent, update_rule, reinit_rule, - max_resets, max_shrink_times, initialization) + linesearch, trustregion, descent, update_rule, + reinit_rule, max_resets, max_shrink_times, initialization) end @inline concrete_jac(::ApproximateJacobianSolveAlgorithm{CJ}) where {CJ} = CJ diff --git a/src/default.jl b/src/default.jl index 4d40fb995..d7f00edb9 100644 --- a/src/default.jl +++ b/src/default.jl @@ -1,6 +1,7 @@ # Poly Algorithms """ - NonlinearSolvePolyAlgorithm(algs, ::Val{pType} = Val(:NLS)) where {pType} + NonlinearSolvePolyAlgorithm(algs, ::Val{pType} = Val(:NLS); + start_index = 1) where {pType} A general way to define PolyAlgorithms for `NonlinearProblem` and `NonlinearLeastSquaresProblem`. This is a container for a tuple of algorithms that will be @@ -15,6 +16,10 @@ residual is returned. `NonlinearLeastSquaresProblem`. This is used to determine the correct problem type to dispatch on. +### Keyword Arguments + + - `start_index`: the index to start at. Defaults to `1`. + ### Example ```julia @@ -25,11 +30,14 @@ alg = NonlinearSolvePolyAlgorithm((NewtonRaphson(), Broyden())) """ struct NonlinearSolvePolyAlgorithm{pType, N, A} <: AbstractNonlinearSolveAlgorithm{:PolyAlg} algs::A + start_index::Int - function NonlinearSolvePolyAlgorithm(algs, ::Val{pType} = Val(:NLS)) where {pType} + function NonlinearSolvePolyAlgorithm( + algs, ::Val{pType} = Val(:NLS); start_index::Int = 1) where {pType} @assert pType ∈ (:NLS, :NLLS) + @assert 0 < start_index ≤ length(algs) algs = Tuple(algs) - return new{pType, length(algs), typeof(algs)}(algs) + return new{pType, length(algs), typeof(algs)}(algs, start_index) end end @@ -73,7 +81,7 @@ end function reinit_cache!(cache::NonlinearSolvePolyAlgorithmCache, args...; kwargs...) foreach(c -> reinit_cache!(c, args...; kwargs...), cache.caches) - cache.current = 1 + cache.current = cache.alg.start_index cache.nsteps = 0 cache.total_time = 0.0 end @@ -91,7 +99,7 @@ for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProb alg.algs), alg, -1, - 1, + alg.start_index, 0, 0.0, maxtime, @@ -140,6 +148,7 @@ end quote fus = tuple($(Tuple(resids)...)) minfu, idx = __findmin(cache.internalnorm, fus) + idx += cache.alg.start_index - 1 stats = __compile_stats(cache.caches[idx]) u = get_u(cache.caches[idx]) retcode = cache.caches[idx].retcode @@ -194,18 +203,22 @@ for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProb @eval begin @generated function SciMLBase.__solve( prob::$probType, alg::$algType{N}, args...; kwargs...) where {N} - calls = [] + calls = [:(current = alg.start_index)] sol_syms = [gensym("sol") for _ in 1:N] for i in 1:N cur_sol = sol_syms[i] push!(calls, quote - $(cur_sol) = SciMLBase.__solve(prob, alg.algs[$(i)], args...; kwargs...) - if SciMLBase.successful_retcode($(cur_sol)) - return SciMLBase.build_solution( - prob, alg, $(cur_sol).u, $(cur_sol).resid; - $(cur_sol).retcode, $(cur_sol).stats, - original = $(cur_sol), trace = $(cur_sol).trace) + if current == $i + $(cur_sol) = SciMLBase.__solve( + prob, alg.algs[$(i)], args...; kwargs...) + if SciMLBase.successful_retcode($(cur_sol)) + return SciMLBase.build_solution( + prob, alg, $(cur_sol).u, $(cur_sol).resid; + $(cur_sol).retcode, $(cur_sol).stats, + original = $(cur_sol), trace = $(cur_sol).trace) + end + current = $(i + 1) end end) end @@ -218,6 +231,7 @@ for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProb push!(calls, quote resids = tuple($(Tuple(resids)...)) minfu, idx = __findmin(DEFAULT_NORM, resids) + idx += alg.start_index - 1 end) for i in 1:N @@ -263,6 +277,7 @@ function RobustMultiNewton(::Type{T} = Float64; concrete_jac = nothing, linsolve algs = (TrustRegion(; concrete_jac, linsolve, precs, autodiff), TrustRegion(; concrete_jac, linsolve, precs, autodiff, radius_update_scheme = RadiusUpdateSchemes.Bastin), + NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), NewtonRaphson(; concrete_jac, linsolve, precs, linesearch = LineSearchesJL(; method = BackTracking()), autodiff), TrustRegion(; concrete_jac, linsolve, precs, @@ -276,7 +291,8 @@ end """ FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS, must_use_jacobian::Val = Val(false), - prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing) where {T} + prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing, + u0_len::Union{Int, Nothing} = nothing) where {T} A polyalgorithm focused on balancing speed and robustness. It first tries less robust methods for more performance and then tries more robust techniques if the faster ones fail. @@ -285,12 +301,19 @@ for more performance and then tries more robust techniques if the faster ones fa - `T`: The eltype of the initial guess. It is only used to check if some of the algorithms are compatible with the problem type. Defaults to `Float64`. + +### Keyword Arguments + + - `u0_len`: The length of the initial guess. If this is `nothing`, then the length of the + initial guess is not checked. If this is an integer and it is less than `25`, we use + jacobian based methods. """ function FastShortcutNonlinearPolyalg( ::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS, must_use_jacobian::Val{JAC} = Val(false), prefer_simplenonlinearsolve::Val{SA} = Val(false), - autodiff = nothing) where {T, JAC, SA} + u0_len::Union{Int, Nothing} = nothing, autodiff = nothing) where {T, JAC, SA} + start_index = 1 if JAC if __is_complex(T) algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff),) @@ -312,6 +335,7 @@ function FastShortcutNonlinearPolyalg( SimpleKlement(), NewtonRaphson(; concrete_jac, linsolve, precs, autodiff)) else + start_index = u0_len !== nothing ? (u0_len ≤ 25 ? 4 : 1) : 1 algs = (SimpleBroyden(), Broyden(; init_jacobian = Val(:true_jacobian), autodiff), SimpleKlement(), @@ -327,6 +351,8 @@ function FastShortcutNonlinearPolyalg( Klement(; linsolve, precs, autodiff), NewtonRaphson(; concrete_jac, linsolve, precs, autodiff)) else + # TODO: This number requires a bit rigorous testing + start_index = u0_len !== nothing ? (u0_len ≤ 25 ? 4 : 1) : 1 algs = (Broyden(; autodiff), Broyden(; init_jacobian = Val(:true_jacobian), autodiff), Klement(; linsolve, precs, autodiff), @@ -339,7 +365,7 @@ function FastShortcutNonlinearPolyalg( end end end - return NonlinearSolvePolyAlgorithm(algs, Val(:NLS)) + return NonlinearSolvePolyAlgorithm(algs, Val(:NLS); start_index) end """ @@ -392,17 +418,19 @@ end ## can use that! function SciMLBase.__init(prob::NonlinearProblem, ::Nothing, args...; kwargs...) must_use_jacobian = Val(prob.f.jac !== nothing) - return SciMLBase.__init( - prob, FastShortcutNonlinearPolyalg(eltype(prob.u0); must_use_jacobian), - args...; kwargs...) + return SciMLBase.__init(prob, + FastShortcutNonlinearPolyalg( + eltype(prob.u0); must_use_jacobian, u0_len = length(prob.u0)), + args...; + kwargs...) end function SciMLBase.__solve(prob::NonlinearProblem, ::Nothing, args...; kwargs...) must_use_jacobian = Val(prob.f.jac !== nothing) prefer_simplenonlinearsolve = Val(prob.u0 isa SArray) return SciMLBase.__solve(prob, - FastShortcutNonlinearPolyalg( - eltype(prob.u0); must_use_jacobian, prefer_simplenonlinearsolve), + FastShortcutNonlinearPolyalg(eltype(prob.u0); must_use_jacobian, + prefer_simplenonlinearsolve, u0_len = length(prob.u0)), args...; kwargs...) end diff --git a/test/core/forward_ad_tests.jl b/test/core/forward_ad_tests.jl index a13259939..24258e157 100644 --- a/test/core/forward_ad_tests.jl +++ b/test/core/forward_ad_tests.jl @@ -64,8 +64,8 @@ end @testitem "ForwardDiff.jl Integration" setup=[ForwardADTesting] begin for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), - PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(), DFSane(), nothing, - NLsolveJL(), CMINPACK(), KINSOL(; globalization_strategy = :LineSearch)) + PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(), DFSane(), + nothing, NLsolveJL(), CMINPACK(), KINSOL(; globalization_strategy = :LineSearch)) us = (2.0, @SVector[1.0, 1.0], [1.0, 1.0], ones(2, 2), @SArray ones(2, 2)) @testset "Scalar AD" begin diff --git a/test/core/rootfind_tests.jl b/test/core/rootfind_tests.jl index 86dc5ff2c..87ba7cd35 100644 --- a/test/core/rootfind_tests.jl +++ b/test/core/rootfind_tests.jl @@ -476,8 +476,8 @@ end @testitem "Broyden" setup=[CoreRootfindTesting] begin @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian) Update Rule: $(update_rule)" for lsmethod in ( - Static(), StrongWolfe(), BackTracking(), HagerZhang(), - MoreThuente(), LiFukushimaLineSearch()), + Static(), StrongWolfe(), BackTracking(), + HagerZhang(), MoreThuente(), LiFukushimaLineSearch()), ad in (AutoFiniteDiff(), AutoZygote()), init_jacobian in (Val(:identity), Val(:true_jacobian)), update_rule in (Val(:good_broyden), Val(:bad_broyden), Val(:diagonal)) @@ -575,8 +575,8 @@ end @testitem "LimitedMemoryBroyden" setup=[CoreRootfindTesting] begin @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in ( - Static(), StrongWolfe(), BackTracking(), HagerZhang(), - MoreThuente(), LiFukushimaLineSearch()), + Static(), StrongWolfe(), BackTracking(), + HagerZhang(), MoreThuente(), LiFukushimaLineSearch()), ad in (AutoFiniteDiff(), AutoZygote()) linesearch = LineSearchesJL(; method = lsmethod, autodiff = ad) diff --git a/test/misc/matrix_resizing_tests.jl b/test/misc/matrix_resizing_tests.jl index 1d17a696a..9210e3e6b 100644 --- a/test/misc/matrix_resizing_tests.jl +++ b/test/misc/matrix_resizing_tests.jl @@ -7,8 +7,8 @@ vecprob = NonlinearProblem(ff, vec(u0), p) prob = NonlinearProblem(ff, u0, p) - for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(), - RobustMultiNewton(), FastShortcutNonlinearPolyalg(), + for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), + PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(), Klement(), LimitedMemoryBroyden(; threshold = 2)) @test vec(solve(prob, alg).u) == solve(vecprob, alg).u end @@ -23,8 +23,8 @@ end vecprob = NonlinearProblem(fiip, vec(u0), p) prob = NonlinearProblem(fiip, u0, p) - for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(), - RobustMultiNewton(), FastShortcutNonlinearPolyalg(), + for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), + PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(), Klement(), LimitedMemoryBroyden(; threshold = 2)) @test vec(solve(prob, alg).u) == solve(vecprob, alg).u end From 390341f3ee3b6894aa0600710aef8901011590a8 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Wed, 6 Mar 2024 10:58:56 -0500 Subject: [PATCH 4/5] Fix tests --- src/default.jl | 6 ++---- src/utils.jl | 1 + test/misc/polyalg_tests.jl | 4 +--- 3 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/default.jl b/src/default.jl index d7f00edb9..cdcfb6825 100644 --- a/src/default.jl +++ b/src/default.jl @@ -142,13 +142,12 @@ end resids = map(x -> Symbol("$(x)_resid"), cache_syms) for (sym, resid) in zip(cache_syms, resids) - push!(calls, :($(resid) = get_fu($(sym)))) + push!(calls, :($(resid) = @isdefined($(sym)) ? get_fu($(sym)) : nothing)) end push!(calls, quote fus = tuple($(Tuple(resids)...)) minfu, idx = __findmin(cache.internalnorm, fus) - idx += cache.alg.start_index - 1 stats = __compile_stats(cache.caches[idx]) u = get_u(cache.caches[idx]) retcode = cache.caches[idx].retcode @@ -225,13 +224,12 @@ for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProb resids = map(x -> Symbol("$(x)_resid"), sol_syms) for (sym, resid) in zip(sol_syms, resids) - push!(calls, :($(resid) = $(sym).resid)) + push!(calls, :($(resid) = @isdefined($(sym)) ? $(sym).resid : nothing)) end push!(calls, quote resids = tuple($(Tuple(resids)...)) minfu, idx = __findmin(DEFAULT_NORM, resids) - idx += alg.start_index - 1 end) for i in 1:N diff --git a/src/utils.jl b/src/utils.jl index aec66bdd4..259ff945a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -99,6 +99,7 @@ function __findmin_caches(f, caches) end function __findmin(f, x) return findmin(x) do xᵢ + xᵢ === nothing && return Inf fx = f(xᵢ) return isnan(fx) ? Inf : fx end diff --git a/test/misc/polyalg_tests.jl b/test/misc/polyalg_tests.jl index 7b7fe6880..c5f1a8eb3 100644 --- a/test/misc/polyalg_tests.jl +++ b/test/misc/polyalg_tests.jl @@ -84,10 +84,8 @@ end Broyden(; autodiff = AutoFiniteDiff()), LimitedMemoryBroyden())) # Uses the `__solve` function - solver = solve(probN; abstol = 1e-9) - @test SciMLBase.successful_retcode(solver) + @test_throws MethodError solve(probN; abstol = 1e-9) @test_throws MethodError solve(probN, RobustMultiNewton(); abstol = 1e-9) - @test SciMLBase.successful_retcode(solver) solver = solve(probN, RobustMultiNewton(; autodiff = AutoFiniteDiff()); abstol = 1e-9) @test SciMLBase.successful_retcode(solver) solver = solve( From 81edf48dfbcdae96f86d381d94692e2ebff4a853 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Wed, 6 Mar 2024 11:39:07 -0500 Subject: [PATCH 5/5] Incorrect Jacobian Size for NLLS --- src/internal/jacobian.jl | 4 +++- test/core/nlls_tests.jl | 19 +++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/internal/jacobian.jl b/src/internal/jacobian.jl index fc4f12f33..3dfb904ee 100644 --- a/src/internal/jacobian.jl +++ b/src/internal/jacobian.jl @@ -83,7 +83,9 @@ function JacobianCache( JacobianOperator(prob, fu, u; jvp_autodiff, vjp_autodiff) else if has_analytic_jac - f.jac_prototype === nothing ? undefmatrix(u) : f.jac_prototype + f.jac_prototype === nothing ? + similar(fu, promote_type(eltype(fu), eltype(u)), length(fu), length(u)) : + copy(f.jac_prototype) elseif f.jac_prototype === nothing init_jacobian(jac_cache; preserve_immutable = Val(true)) else diff --git a/test/core/nlls_tests.jl b/test/core/nlls_tests.jl index 3722ca671..6b7b018cd 100644 --- a/test/core/nlls_tests.jl +++ b/test/core/nlls_tests.jl @@ -94,3 +94,22 @@ end @test maximum(abs, sol.resid) < 1e-6 end end + +@testitem "NLLS Analytic Jacobian" begin + dataIn = 1:10 + f(x, p) = x[1] * dataIn .^ 2 .+ x[2] * dataIn .+ x[3] + dataOut = f([1, 2, 3], nothing) + 0.1 * randn(10, 1) + + resid(x, p) = f(x, p) - dataOut + jac(x, p) = [dataIn .^ 2 dataIn ones(10, 1)] + x0 = [1, 1, 1] + + prob = NonlinearLeastSquaresProblem(resid, x0) + sol1 = solve(prob) + + nlfunc = NonlinearFunction(resid; jac) + prob = NonlinearLeastSquaresProblem(nlfunc, x0) + sol2 = solve(prob) + + @test sol1.u ≈ sol2.u +end