|
1 | 1 | using ModelingToolkit, NonlinearSolve, OrdinaryDiffEq, Sundials, SciMLBase, Test
|
2 | 2 | using SymbolicIndexingInterface
|
3 | 3 | using ModelingToolkit: t_nounits as t, D_nounits as D
|
4 |
| - |
| 4 | +using StochasticDiffEq, OrdinaryDiffEq, NonlinearSolve, SymbolicIndexingInterface, |
| 5 | + LinearAlgebra, Test |
5 | 6 | @testset "CheckInit" begin
|
6 | 7 | abstol = 1e-10
|
7 | 8 | @testset "Sundials + DAEProblem" begin
|
|
76 | 77 | prob = ODEProblem(pend, [x => 1, y => 0, λ => 2], (0.0, 1.5), [g => 1])
|
77 | 78 | @test occursin("Initialization status: OVERDETERMINED", sprint(summary, prob))
|
78 | 79 | end
|
| 80 | + |
| 81 | + |
| 82 | +@testset "CheckInit" begin |
| 83 | + @testset "ODEProblem" begin |
| 84 | + function rhs(u, p, t) |
| 85 | + return [u[1] * t, u[1]^2 - u[2]^2] |
| 86 | + end |
| 87 | + function rhs!(du, u, p, t) |
| 88 | + du[1] = u[1] * t |
| 89 | + du[2] = u[1]^2 - u[2]^2 |
| 90 | + end |
| 91 | + |
| 92 | + oopfn = ODEFunction{false}(rhs, mass_matrix = [1 0; 0 0]) |
| 93 | + iipfn = ODEFunction{true}(rhs!, mass_matrix = [1 0; 0 0]) |
| 94 | + |
| 95 | + @testset "Inplace = $(SciMLBase.isinplace(f))" for f in [oopfn, iipfn] |
| 96 | + prob = ODEProblem(f, [1.0, 1.0], (0.0, 1.0)) |
| 97 | + integ = init(prob) |
| 98 | + u0, _, success = SciMLBase.get_initial_values( |
| 99 | + prob, integ, f, SciMLBase.CheckInit(), |
| 100 | + Val(SciMLBase.isinplace(f)); abstol = 1e-10) |
| 101 | + @test success |
| 102 | + @test u0 == prob.u0 |
| 103 | + |
| 104 | + integ.u[2] = 2.0 |
| 105 | + @test_throws SciMLBase.CheckInitFailureError SciMLBase.get_initial_values( |
| 106 | + prob, integ, f, SciMLBase.CheckInit(), |
| 107 | + Val(SciMLBase.isinplace(f)); abstol = 1e-10) |
| 108 | + end |
| 109 | + |
| 110 | + @testset "With I mass matrix" begin |
| 111 | + function rhs(u, p, t) |
| 112 | + return u |
| 113 | + end |
| 114 | + prob = ODEProblem(ODEFunction(rhs; mass_matrix = I), ones(2), (0.0, 1.0)) |
| 115 | + integ = init(prob) |
| 116 | + u0, _, success = SciMLBase.get_initial_values( |
| 117 | + prob, integ, prob.f, SciMLBase.CheckInit(), |
| 118 | + Val(false); abstol = 1e-10 |
| 119 | + ) |
| 120 | + @test success |
| 121 | + @test u0 == prob.u0 |
| 122 | + end |
| 123 | + end |
| 124 | + |
| 125 | + @testset "DAEProblem" begin |
| 126 | + function daerhs(du, u, p, t) |
| 127 | + return [du[1] - u[1] * t - p, u[1]^2 - u[2]^2] |
| 128 | + end |
| 129 | + function daerhs!(resid, du, u, p, t) |
| 130 | + resid[1] = du[1] - u[1] * t - p |
| 131 | + resid[2] = u[1]^2 - u[2]^2 |
| 132 | + end |
| 133 | + |
| 134 | + oopfn = DAEFunction{false}(daerhs) |
| 135 | + iipfn = DAEFunction{true}(daerhs!) |
| 136 | + |
| 137 | + @testset "Inplace = $(SciMLBase.isinplace(f))" for f in [oopfn, iipfn] |
| 138 | + prob = DAEProblem(f, [1.0, 0.0], [1.0, 1.0], (0.0, 1.0), 1.0) |
| 139 | + integ = init(prob, DImplicitEuler()) |
| 140 | + u0, _, success = SciMLBase.get_initial_values( |
| 141 | + prob, integ, f, SciMLBase.CheckInit(), |
| 142 | + Val(SciMLBase.isinplace(f)); abstol = 1e-10) |
| 143 | + @test success |
| 144 | + @test u0 == prob.u0 |
| 145 | + |
| 146 | + integ.u[2] = 2.0 |
| 147 | + @test_throws SciMLBase.CheckInitFailureError SciMLBase.get_initial_values( |
| 148 | + prob, integ, f, SciMLBase.CheckInit(), |
| 149 | + Val(SciMLBase.isinplace(f)); abstol = 1e-10) |
| 150 | + |
| 151 | + integ.u[2] = 1.0 |
| 152 | + integ.du[1] = 2.0 |
| 153 | + @test_throws SciMLBase.CheckInitFailureError SciMLBase.get_initial_values( |
| 154 | + prob, integ, f, SciMLBase.CheckInit(), |
| 155 | + Val(SciMLBase.isinplace(f)); abstol = 1e-10) |
| 156 | + end |
| 157 | + end |
| 158 | + |
| 159 | + @testset "SDEProblem" begin |
| 160 | + mm_A = [1 0 0; 0 1 0; 0 0 0] |
| 161 | + function sdef!(du, u, p, t) |
| 162 | + du[1] = u[1] |
| 163 | + du[2] = u[2] |
| 164 | + du[3] = u[1] + u[2] + u[3] - 1 |
| 165 | + end |
| 166 | + function sdef(u, p, t) |
| 167 | + du = similar(u) |
| 168 | + sdef!(du, u, p, t) |
| 169 | + du |
| 170 | + end |
| 171 | + |
| 172 | + function g!(du, u, p, t) |
| 173 | + @. du = 0.1 |
| 174 | + end |
| 175 | + function g(u, p, t) |
| 176 | + du = similar(u) |
| 177 | + g!(du, u, p, t) |
| 178 | + du |
| 179 | + end |
| 180 | + iipfn = SDEFunction{true}(sdef!, g!; mass_matrix = mm_A) |
| 181 | + oopfn = SDEFunction{false}(sdef, g; mass_matrix = mm_A) |
| 182 | + |
| 183 | + @testset "Inplace = $(SciMLBase.isinplace(f))" for f in [oopfn, iipfn] |
| 184 | + prob = SDEProblem(f, [1.0, 1.0, -1.0], (0.0, 1.0)) |
| 185 | + integ = init(prob, ImplicitEM()) |
| 186 | + u0, _, success = SciMLBase.get_initial_values( |
| 187 | + prob, integ, f, SciMLBase.CheckInit(), |
| 188 | + Val(SciMLBase.isinplace(f)); abstol = 1e-10) |
| 189 | + @test success |
| 190 | + @test u0 == prob.u0 |
| 191 | + |
| 192 | + integ.u[2] = 2.0 |
| 193 | + @test_throws SciMLBase.CheckInitFailureError SciMLBase.get_initial_values( |
| 194 | + prob, integ, f, SciMLBase.CheckInit(), |
| 195 | + Val(SciMLBase.isinplace(f)); abstol = 1e-10) |
| 196 | + end |
| 197 | + end |
| 198 | +end |
| 199 | + |
| 200 | +@testset "OverrideInit" begin |
| 201 | + function rhs2(u, p, t) |
| 202 | + return [u[1] * t + p, u[1]^2 - u[2]^2] |
| 203 | + end |
| 204 | + |
| 205 | + @testset "No-op without `initialization_data`" begin |
| 206 | + prob = ODEProblem(rhs2, [1.0, 2.0], (0.0, 1.0), 1.0) |
| 207 | + @test SciMLBase.initialization_status(prob) === nothing |
| 208 | + integ = init(prob) |
| 209 | + integ.u[2] = 3.0 |
| 210 | + u0, p, success = SciMLBase.get_initial_values( |
| 211 | + prob, integ, prob.f, SciMLBase.OverrideInit(), Val(false)) |
| 212 | + @test u0 ≈ [1.0, 3.0] |
| 213 | + @test success |
| 214 | + end |
| 215 | + |
| 216 | + # unknowns are u[2], p. Parameter is u[1] |
| 217 | + initprob = NonlinearProblem([1.0, 1.0], [1.0]) do x, _u1 |
| 218 | + u2, p = x |
| 219 | + u1 = _u1[1] |
| 220 | + return [u1^2 - u2^2, p^2 - 2p + 1] |
| 221 | + end |
| 222 | + update_initializeprob! = function (iprob, integ) |
| 223 | + iprob.p[1] = integ.u[1] |
| 224 | + end |
| 225 | + initprobmap = function (nlsol) |
| 226 | + return [parameter_values(nlsol)[1], nlsol.u[1]] |
| 227 | + end |
| 228 | + initprobpmap = function (_, nlsol) |
| 229 | + return nlsol.u[2] |
| 230 | + end |
| 231 | + initialization_data = SciMLBase.OverrideInitData( |
| 232 | + initprob, update_initializeprob!, initprobmap, initprobpmap) |
| 233 | + fn = ODEFunction(rhs2; initialization_data) |
| 234 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 235 | + @test SciMLBase.initialization_status(prob) == SciMLBase.FULLY_DETERMINED |
| 236 | + integ = init(prob; initializealg = NoInit()) |
| 237 | + |
| 238 | + @testset "Errors without `nlsolve_alg`" begin |
| 239 | + @test_throws SciMLBase.OverrideInitMissingAlgorithm SciMLBase.get_initial_values( |
| 240 | + prob, integ, fn, SciMLBase.OverrideInit(), Val(false)) |
| 241 | + end |
| 242 | + |
| 243 | + abstol = 1e-10 |
| 244 | + reltol = 1e-10 |
| 245 | + @testset "Solves" begin |
| 246 | + @testset "with explicit alg" begin |
| 247 | + u0, p, success = SciMLBase.get_initial_values( |
| 248 | + prob, integ, fn, SciMLBase.OverrideInit(), |
| 249 | + Val(false); nlsolve_alg = NewtonRaphson(), abstol, reltol) |
| 250 | + |
| 251 | + @test u0 ≈ [2.0, 2.0] |
| 252 | + @test p ≈ 1.0 |
| 253 | + @test success |
| 254 | + |
| 255 | + initprob.p[1] = 1.0 |
| 256 | + end |
| 257 | + @testset "with alg in `OverrideInit`" begin |
| 258 | + u0, p, success = SciMLBase.get_initial_values( |
| 259 | + prob, integ, fn, |
| 260 | + SciMLBase.OverrideInit(; nlsolve = NewtonRaphson(), abstol, reltol), |
| 261 | + Val(false)) |
| 262 | + |
| 263 | + @test u0 ≈ [2.0, 2.0] |
| 264 | + @test p ≈ 1.0 |
| 265 | + @test success |
| 266 | + |
| 267 | + initprob.p[1] = 1.0 |
| 268 | + end |
| 269 | + @testset "with trivial problem and no alg" begin |
| 270 | + iprob = NonlinearProblem((u, p) -> 0.0, nothing, 1.0) |
| 271 | + iprobmap = (_) -> [1.0, 1.0] |
| 272 | + initdata = SciMLBase.OverrideInitData(iprob, nothing, iprobmap, nothing) |
| 273 | + _fn = ODEFunction(rhs2; initialization_data = initdata) |
| 274 | + _prob = ODEProblem(_fn, [2.0, 0.0], (0.0, 1.0), 1.0) |
| 275 | + _integ = init(_prob; initializealg = NoInit()) |
| 276 | + |
| 277 | + u0, p, success = SciMLBase.get_initial_values( |
| 278 | + _prob, _integ, _fn, SciMLBase.OverrideInit(), Val(false); abstol, reltol) |
| 279 | + |
| 280 | + @test u0 ≈ [1.0, 1.0] |
| 281 | + @test p ≈ 1.0 |
| 282 | + @test success |
| 283 | + end |
| 284 | + @testset "with kwargs provided to `get_initial_values`" begin |
| 285 | + u0, p, success = SciMLBase.get_initial_values( |
| 286 | + prob, integ, fn, SciMLBase.OverrideInit(), |
| 287 | + Val(false); nlsolve_alg = NewtonRaphson(), abstol, reltol, u0 = [-1.0, 1.0]) |
| 288 | + @test u0 ≈ [2.0, -2.0] |
| 289 | + @test p ≈ 1.0 |
| 290 | + @test success |
| 291 | + end |
| 292 | + end |
| 293 | + |
| 294 | + @testset "Solves with non-integrator value provider" begin |
| 295 | + _integ = ProblemState(; u = integ.u, p = parameter_values(integ), t = integ.t) |
| 296 | + u0, p, success = SciMLBase.get_initial_values( |
| 297 | + prob, _integ, fn, SciMLBase.OverrideInit(), |
| 298 | + Val(false); nlsolve_alg = NewtonRaphson(), abstol, reltol) |
| 299 | + |
| 300 | + @test u0 ≈ [2.0, 2.0] |
| 301 | + @test p ≈ 1.0 |
| 302 | + @test success |
| 303 | + |
| 304 | + initprob.p[1] = 1.0 |
| 305 | + end |
| 306 | + |
| 307 | + @testset "Solves without `update_initializeprob!`" begin |
| 308 | + initdata = SciMLBase.@set initialization_data.update_initializeprob! = nothing |
| 309 | + fn = ODEFunction(rhs2; initialization_data = initdata) |
| 310 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 311 | + integ = init(prob; initializealg = NoInit()) |
| 312 | + |
| 313 | + u0, p, success = SciMLBase.get_initial_values( |
| 314 | + prob, integ, fn, SciMLBase.OverrideInit(), |
| 315 | + Val(false); nlsolve_alg = NewtonRaphson(), abstol, reltol) |
| 316 | + @test u0 ≈ [1.0, 1.0] |
| 317 | + @test p ≈ 1.0 |
| 318 | + @test success |
| 319 | + end |
| 320 | + |
| 321 | + @testset "Solves without `initializeprobmap`" begin |
| 322 | + initdata = SciMLBase.@set initialization_data.initializeprobmap = nothing |
| 323 | + fn = ODEFunction(rhs2; initialization_data = initdata) |
| 324 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 325 | + integ = init(prob; initializealg = NoInit()) |
| 326 | + |
| 327 | + u0, p, success = SciMLBase.get_initial_values( |
| 328 | + prob, integ, fn, SciMLBase.OverrideInit(), |
| 329 | + Val(false); nlsolve_alg = NewtonRaphson(), abstol, reltol) |
| 330 | + |
| 331 | + @test u0 ≈ [2.0, 0.0] |
| 332 | + @test p ≈ 1.0 |
| 333 | + @test success |
| 334 | + end |
| 335 | + |
| 336 | + @testset "Solves without `initializeprobpmap`" begin |
| 337 | + initdata = SciMLBase.@set initialization_data.initializeprobpmap = nothing |
| 338 | + fn = ODEFunction(rhs2; initialization_data = initdata) |
| 339 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 340 | + integ = init(prob; initializealg = NoInit()) |
| 341 | + |
| 342 | + u0, p, success = SciMLBase.get_initial_values( |
| 343 | + prob, integ, fn, SciMLBase.OverrideInit(), |
| 344 | + Val(false); nlsolve_alg = NewtonRaphson(), abstol, reltol) |
| 345 | + |
| 346 | + @test u0 ≈ [2.0, 2.0] |
| 347 | + @test p ≈ 0.0 |
| 348 | + @test success |
| 349 | + end |
| 350 | + |
| 351 | + @testset "Initialization status for `SCCNonlinearProblem`" begin |
| 352 | + initprob = SCCNonlinearProblem([initprob], [Returns(nothing)]) |
| 353 | + initialization_data = SciMLBase.OverrideInitData( |
| 354 | + initprob, nothing, nothing, nothing) |
| 355 | + fn = ODEFunction(rhs2; initialization_data) |
| 356 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 357 | + @test SciMLBase.initialization_status(prob) == SciMLBase.FULLY_DETERMINED |
| 358 | + end |
| 359 | + |
| 360 | + @testset "Trivial initialization" begin |
| 361 | + initprob = NonlinearProblem(Returns(nothing), nothing, [1.0]) |
| 362 | + update_initializeprob! = function (iprob, integ) |
| 363 | + # just to access the current time and use it as a number, so this errors |
| 364 | + # if run on a problem with `current_time(prob) === nothing` |
| 365 | + iprob.p[1] = current_time(integ) + 1 |
| 366 | + iprob.p[1] = state_values(integ)[1] |
| 367 | + end |
| 368 | + initprobmap = function (nlsol) |
| 369 | + u1 = parameter_values(nlsol)[1] |
| 370 | + return [u1, u1] |
| 371 | + end |
| 372 | + initprobpmap = function (_, nlsol) |
| 373 | + return 0.0 |
| 374 | + end |
| 375 | + initialization_data = SciMLBase.OverrideInitData( |
| 376 | + initprob, update_initializeprob!, initprobmap, initprobpmap) |
| 377 | + fn = ODEFunction(rhs2; initialization_data) |
| 378 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 379 | + @test SciMLBase.initialization_status(prob) == SciMLBase.FULLY_DETERMINED |
| 380 | + integ = init(prob; initializealg = NoInit()) |
| 381 | + |
| 382 | + u0, p, success = SciMLBase.get_initial_values( |
| 383 | + prob, integ, fn, SciMLBase.OverrideInit(), Val(false) |
| 384 | + ) |
| 385 | + @test u0 ≈ [2.0, 2.0] |
| 386 | + @test p ≈ 0.0 |
| 387 | + @test success |
| 388 | + |
| 389 | + @testset "Doesn't run in `remake` if `tspan == (nothing, nothing)`" begin |
| 390 | + prob = ODEProblem(fn, [2.0, 0.0], (nothing, nothing), 0.0) |
| 391 | + @test_nowarn remake(prob) |
| 392 | + end |
| 393 | + end |
| 394 | + |
| 395 | + @testset "Initialization status for overdetermined case" begin |
| 396 | + initfn = NonlinearFunction(; resid_prototype = ones(3)) do u, p |
| 397 | + return [u[1] - 1.0, u[2] - 1.0, u[1] * u[2] - 1.0] |
| 398 | + end |
| 399 | + initprob = NonlinearLeastSquaresProblem(initfn, ones(2)) |
| 400 | + initialization_data = SciMLBase.OverrideInitData( |
| 401 | + initprob, nothing, nothing, nothing) |
| 402 | + fn = ODEFunction(rhs2; initialization_data) |
| 403 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 404 | + @test SciMLBase.initialization_status(prob) == SciMLBase.OVERDETERMINED |
| 405 | + end |
| 406 | + |
| 407 | + @testset "Initialization status for underdetermined case" begin |
| 408 | + initfn = NonlinearFunction(; resid_prototype = ones(1)) do u, p |
| 409 | + return [u[1] - 1.0] |
| 410 | + end |
| 411 | + initprob = NonlinearLeastSquaresProblem(initfn, ones(2)) |
| 412 | + initialization_data = SciMLBase.OverrideInitData( |
| 413 | + initprob, nothing, nothing, nothing) |
| 414 | + fn = ODEFunction(rhs2; initialization_data) |
| 415 | + prob = ODEProblem(fn, [2.0, 0.0], (0.0, 1.0), 0.0) |
| 416 | + @test SciMLBase.initialization_status(prob) == SciMLBase.UNDERDETERMINED |
| 417 | + end |
| 418 | +end |
| 419 | + |
| 420 | +@testset "NoInit" begin |
| 421 | + prob = ODEProblem(ones(2), (0.0, 1.0), ones(2)) do u, p, t |
| 422 | + return u |
| 423 | + end |
| 424 | + u, p, success = SciMLBase.get_initial_values( |
| 425 | + prob, prob, prob.f, SciMLBase.NoInit(), Val(true)) |
| 426 | + @test u == ones(2) |
| 427 | + @test p == ones(2) |
| 428 | + @test success |
| 429 | +end |
0 commit comments