Skip to content

Commit

Permalink
Run NLLS tests only in ≥ 1.10-
Browse files Browse the repository at this point in the history
  • Loading branch information
@avik-pal
avik-pal committed Nov 2, 2023
1 parent 35c0050 commit 5011657
Show file tree
Hide file tree
Showing 4 changed files with 117 additions and 109 deletions.
24 changes: 12 additions & 12 deletions src/solve/mirk.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -303,31 +303,31 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
jac_prototype = vcat(init_jacobian(cache_bc), init_jacobian(cache_collocation))

jac = if iip
(J, u, p) -> __mirk_mpoint_jacobian!(J, u, p, jac_alg.bc_diffmode,
(J, u, p) -> __mirk_mpoint_jacobian!(J, u, jac_alg.bc_diffmode,
jac_alg.nonbc_diffmode, cache_bc, cache_collocation, loss_bcₚ,
loss_collocationₚ, resid_bc, resid_collocation, cache.M, L)
loss_collocationₚ, resid_bc, resid_collocation, L)
else
(u, p) -> __mirk_mpoint_jacobian(u, p, jac_prototype, jac_alg.bc_diffmode,
(u, p) -> __mirk_mpoint_jacobian(jac_prototype, u, jac_alg.bc_diffmode,
jac_alg.nonbc_diffmode, cache_bc, cache_collocation, loss_bcₚ,
loss_collocationₚ, cache.M, L)
loss_collocationₚ, L)
end

nlf = NonlinearFunction{iip}(loss; resid_prototype = vcat(resid_bc, resid_collocation),
jac, jac_prototype)
return (L == cache.M ? NonlinearProblem : NonlinearLeastSquaresProblem)(nlf, y, cache.p)
end

function __mirk_mpoint_jacobian!(J, x, p, bc_diffmode, nonbc_diffmode, bc_diffcache,
function __mirk_mpoint_jacobian!(J, x, bc_diffmode, nonbc_diffmode, bc_diffcache,
nonbc_diffcache, loss_bc::BC, loss_collocation::C, resid_bc, resid_collocation,
M::Int, L::Int) where {BC, C}
L::Int) where {BC, C}
sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, resid_bc, x)
sparse_jacobian!(@view(J[(L + 1):end, :]), nonbc_diffmode, nonbc_diffcache,
loss_collocation, resid_collocation, x)
return nothing
end

function __mirk_mpoint_jacobian(x, p, J, bc_diffmode, nonbc_diffmode, bc_diffcache,
nonbc_diffcache, loss_bc::BC, loss_collocation::C, M::Int, L::Int) where {BC, C}
function __mirk_mpoint_jacobian(J, x, bc_diffmode, nonbc_diffmode, bc_diffcache,
nonbc_diffcache, loss_bc::BC, loss_collocation::C, L::Int) where {BC, C}
sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, x)
sparse_jacobian!(@view(J[(L + 1):end, :]), nonbc_diffmode, nonbc_diffcache,
loss_collocation, x)
Expand Down Expand Up @@ -358,10 +358,10 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
jac_prototype = init_jacobian(diffcache)

jac = if iip
(J, u, p) -> __mirk_2point_jacobian!(J, u, p, jac_alg.diffmode, diffcache, lossₚ,
(J, u, p) -> __mirk_2point_jacobian!(J, u, jac_alg.diffmode, diffcache, lossₚ,
resid)
else
(u, p) -> __mirk_2point_jacobian(u, p, jac_prototype, jac_alg.diffmode, diffcache,
(u, p) -> __mirk_2point_jacobian(u, jac_prototype, jac_alg.diffmode, diffcache,
lossₚ)
end

Expand All @@ -370,12 +370,12 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
return (L == cache.M ? NonlinearProblem : NonlinearLeastSquaresProblem)(nlf, y, cache.p)
end

function __mirk_2point_jacobian!(J, x, p, diffmode, diffcache, loss_fn::L, resid) where {L}
function __mirk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid) where {L}
sparse_jacobian!(J, diffmode, diffcache, loss_fn, resid, x)
return J
end

function __mirk_2point_jacobian(x, p, J, diffmode, diffcache, loss_fn::L) where {L}
function __mirk_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L) where {L}
sparse_jacobian!(J, diffmode, diffcache, loss_fn, x)
return J
end
11 changes: 9 additions & 2 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,11 @@ const GROUP = uppercase(get(ENV, "GROUP", "ALL"))
@time @safetestset "Orbital" begin
include("shooting/orbital.jl")
end
if VERSION ≥ v"1.10-"
@time @safetestset "Shooting NLLS Tests" begin
include("shooting/nonlinear_least_squares.jl")
end
end
end
end

Expand All @@ -31,8 +36,10 @@ const GROUP = uppercase(get(ENV, "GROUP", "ALL"))
@time @safetestset "Interpolation Tests" begin
include("mirk/interpolation_test.jl")
end
@time @safetestset "MIRK Nonlinear Least Squares Tests" begin
include("mirk/nonlinear_least_squares.jl")
if VERSION ≥ v"1.10-"
@time @safetestset "MIRK NLLS Tests" begin
include("mirk/nonlinear_least_squares.jl")
end
end
end
end
Expand Down
96 changes: 96 additions & 0 deletions test/shooting/nonlinear_least_squares.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test

@testset "Overconstrained BVP" begin
SOLVERS = [
Shooting(Tsit5(); nlsolve = LevenbergMarquardt()),
Shooting(Tsit5(); nlsolve = GaussNewton()),
MultipleShooting(10, Tsit5(); nlsolve = LevenbergMarquardt()),
MultipleShooting(10, Tsit5(); nlsolve = GaussNewton())]

# OOP MP-BVP
f1(u, p, t) = [u[2], -u[1]]

function bc1(sol, p, t)
t₁, t₂ = extrema(t)
solₜ₁ = sol(t₁)
solₜ₂ = sol(t₂)
solₜ₃ = sol((t₁ + t₂) / 2)
# We know that this overconstrained system has a solution
return [solₜ₁[1], solₜ₂[1] - 1, solₜ₃[1] - 0.51735, solₜ₃[2] + 1.92533]
end

tspan = (0.0, 100.0)
u0 = [0.0, 1.0]

bvp1 = BVProblem(BVPFunction{false}(f1, bc1; bcresid_prototype = zeros(4)), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp1, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
@test norm(bc1(sol, nothing, sol.t)) < 1e-4
end

# IIP MP-BVP
function f1!(du, u, p, t)
du[1] = u[2]
du[2] = -u[1]
return nothing
end

function bc1!(resid, sol, p, t)
(t₁, t₂) = extrema(t)
solₜ₁ = sol(t₁)
solₜ₂ = sol(t₂)
solₜ₃ = sol((t₁ + t₂) / 2)
# We know that this overconstrained system has a solution
resid[1] = solₜ₁[1]
resid[2] = solₜ₂[1] - 1
resid[3] = solₜ₃[1] - 0.51735
resid[4] = solₜ₃[2] + 1.92533
return nothing
end

bvp2 = BVProblem(BVPFunction{true}(f1!, bc1!; bcresid_prototype = zeros(4)), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp2, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
resid_f = Array{Float64}(undef, 4)
bc1!(resid_f, sol, nothing, sol.t)
@test norm(resid_f) < 1e-4
end

# OOP TP-BVP
bc1a(ua, p) = [ua[1]]
bc1b(ub, p) = [ub[1] - 1, ub[2] + 1.729109]

bvp3 = TwoPointBVProblem(BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true),
bcresid_prototype = (zeros(1), zeros(2))), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp3, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
@test norm(vcat(bc1a(sol(0.0), nothing), bc1b(sol(100.0), nothing))) < 1e-4
end

# IIP TP-BVP
bc1a!(resid, ua, p) = (resid[1] = ua[1])
bc1b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109)

bvp4 = TwoPointBVProblem(BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true),
bcresid_prototype = (zeros(1), zeros(2))), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp4, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
resida = Array{Float64}(undef, 1)
residb = Array{Float64}(undef, 2)
bc1a!(resida, sol(0.0), nothing)
bc1b!(residb, sol(100.0), nothing)
@test norm(vcat(resida, residb)) < 1e-4
end
end
95 changes: 0 additions & 95 deletions test/shooting/shooting_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,101 +80,6 @@ using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test
end
end

@testset "Overconstrained BVP" begin
SOLVERS = [
Shooting(Tsit5(); nlsolve = LevenbergMarquardt()),
Shooting(Tsit5(); nlsolve = GaussNewton()),
MultipleShooting(10, Tsit5(); nlsolve = LevenbergMarquardt()),
MultipleShooting(10, Tsit5(); nlsolve = GaussNewton())]

# OOP MP-BVP
f1(u, p, t) = [u[2], -u[1]]

function bc1(sol, p, t)
t₁, t₂ = extrema(t)
solₜ₁ = sol(t₁)
solₜ₂ = sol(t₂)
solₜ₃ = sol((t₁ + t₂) / 2)
# We know that this overconstrained system has a solution
return [solₜ₁[1], solₜ₂[1] - 1, solₜ₃[1] - 0.51735, solₜ₃[2] + 1.92533]
end

tspan = (0.0, 100.0)
u0 = [0.0, 1.0]

bvp1 = BVProblem(BVPFunction{false}(f1, bc1; bcresid_prototype = zeros(4)), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp1, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
@test norm(bc1(sol, nothing, sol.t)) < 1e-4
end

# IIP MP-BVP
function f1!(du, u, p, t)
du[1] = u[2]
du[2] = -u[1]
return nothing
end

function bc1!(resid, sol, p, t)
(t₁, t₂) = extrema(t)
solₜ₁ = sol(t₁)
solₜ₂ = sol(t₂)
solₜ₃ = sol((t₁ + t₂) / 2)
# We know that this overconstrained system has a solution
resid[1] = solₜ₁[1]
resid[2] = solₜ₂[1] - 1
resid[3] = solₜ₃[1] - 0.51735
resid[4] = solₜ₃[2] + 1.92533
return nothing
end

bvp2 = BVProblem(BVPFunction{true}(f1!, bc1!; bcresid_prototype = zeros(4)), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp2, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
resid_f = Array{Float64}(undef, 4)
bc1!(resid_f, sol, nothing, sol.t)
@test norm(resid_f) < 1e-4
end

# OOP TP-BVP
bc1a(ua, p) = [ua[1]]
bc1b(ub, p) = [ub[1] - 1, ub[2] + 1.729109]

bvp3 = TwoPointBVProblem(BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true),
bcresid_prototype = (zeros(1), zeros(2))), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp3, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
@test norm(vcat(bc1a(sol(0.0), nothing), bc1b(sol(100.0), nothing))) < 1e-4
end

# IIP TP-BVP
bc1a!(resid, ua, p) = (resid[1] = ua[1])
bc1b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109)

bvp4 = TwoPointBVProblem(BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true),
bcresid_prototype = (zeros(1), zeros(2))), u0, tspan)

for solver in SOLVERS
@time sol = solve(bvp4, solver;
nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000),
verbose = false)
resida = Array{Float64}(undef, 1)
residb = Array{Float64}(undef, 2)
bc1a!(resida, sol(0.0), nothing)
bc1b!(residb, sol(100.0), nothing)
@test norm(vcat(resida, residb)) < 1e-4
end
end

@testset "Shooting with Complex Values" begin
# Test for complex values
function f1!(du, u, p, t)
Expand Down

0 comments on commit 5011657

Please sign in to comment.