| Line | Exclusive | Inclusive | Code |
|---|---|---|---|
| 1 | """ | ||
| 2 | step!(integ::DEIntegrator [, dt [, stop_at_tdt]]) | ||
| 3 | |||
| 4 | Perform one (successful) step on the integrator. | ||
| 5 | |||
| 6 | Alternative, if a `dt` is given, then `step!` the integrator until | ||
| 7 | there is a temporal difference `≥ dt` in `integ.t`. When `true` is | ||
| 8 | passed to the optional third argument, the integrator advances exactly | ||
| 9 | `dt`. | ||
| 10 | """ | ||
| 11 | function step!(d::DEIntegrator) | ||
| 12 | error("Integrator stepping is not implemented") | ||
| 13 | end | ||
| 14 | |||
| 15 | """ | ||
| 16 | resize!(integrator::DEIntegrator,k::Int) | ||
| 17 | |||
| 18 | Resizes the DE to a size `k`. This chops off the end of the array, or adds blank values at the end, depending on whether | ||
| 19 | `k > length(integrator.u)`. | ||
| 20 | """ | ||
| 21 | function Base.resize!(i::DEIntegrator, ii::Int) | ||
| 22 | error("resize!: method has not been implemented for the integrator") | ||
| 23 | end | ||
| 24 | |||
| 25 | """ | ||
| 26 | deleteat!(integrator::DEIntegrator,idxs) | ||
| 27 | |||
| 28 | Shrinks the ODE by deleting the `idxs` components. | ||
| 29 | """ | ||
| 30 | function Base.deleteat!(i::DEIntegrator, ii) | ||
| 31 | error("deleteat!: method has not been implemented for the integrator") | ||
| 32 | end | ||
| 33 | |||
| 34 | """ | ||
| 35 | addat!(integrator::DEIntegrator,idxs,val) | ||
| 36 | |||
| 37 | Grows the ODE by adding the `idxs` components. Must be contiguous indices. | ||
| 38 | """ | ||
| 39 | function addat!(i::DEIntegrator, idxs, val = zeros(length(idxs))) | ||
| 40 | error("addat!: method has not been implemented for the integrator") | ||
| 41 | end | ||
| 42 | |||
| 43 | """ | ||
| 44 | get_tmp_cache(i::DEIntegrator) | ||
| 45 | |||
| 46 | Returns a tuple of internal cache vectors which are safe to use as temporary arrays. This should be used | ||
| 47 | for integrator interface and callbacks which need arrays to write into in order to be non-allocating. | ||
| 48 | The length of the tuple is dependent on the method. | ||
| 49 | """ | ||
| 50 | function get_tmp_cache(i::DEIntegrator) | ||
| 51 | error("get_tmp_cache!: method has not been implemented for the integrator") | ||
| 52 | end | ||
| 53 | function user_cache(i::DEIntegrator) | ||
| 54 | error("user_cache: method has not been implemented for the integrator") | ||
| 55 | end | ||
| 56 | function u_cache(i::DEIntegrator) | ||
| 57 | error("u_cache: method has not been implemented for the integrator") | ||
| 58 | end | ||
| 59 | function du_cache(i::DEIntegrator) | ||
| 60 | error("du_cache: method has not been implemented for the integrator") | ||
| 61 | end | ||
| 62 | ratenoise_cache(i::DEIntegrator) = () | ||
| 63 | rand_cache(i::DEIntegrator) = () | ||
| 64 | |||
| 65 | """ | ||
| 66 | full_cache(i::DEIntegrator) | ||
| 67 | |||
| 68 | Returns an iterator over the cache arrays of the method. This can be used to change internal values as needed. | ||
| 69 | """ | ||
| 70 | function full_cache(i::DEIntegrator) | ||
| 71 | error("full_cache: method has not been implemented for the integrator") | ||
| 72 | end | ||
| 73 | |||
| 74 | """ | ||
| 75 | resize_non_user_cache!(integrator::DEIntegrator,k::Int) | ||
| 76 | |||
| 77 | Resizes the non-user facing caches to be compatible with a DE of size `k`. This includes resizing Jacobian caches. | ||
| 78 | |||
| 79 | !!! note | ||
| 80 | |||
| 81 | In many cases, [`resize!`](@ref) simply resizes [`full_cache`](@ref) variables and then | ||
| 82 | calls this function. This finer control is required for some `AbstractArray` | ||
| 83 | operations. | ||
| 84 | """ | ||
| 85 | function resize_non_user_cache!(i::DEIntegrator, ii::Int) | ||
| 86 | error("resize_non_user_cache!: method has not been implemented for the integrator") | ||
| 87 | end | ||
| 88 | |||
| 89 | """ | ||
| 90 | deleteat_non_user_cache!(integrator::DEIntegrator,idxs) | ||
| 91 | |||
| 92 | [`deleteat!`](@ref)s the non-user facing caches at indices `idxs`. This includes resizing Jacobian caches. | ||
| 93 | |||
| 94 | !!! note | ||
| 95 | |||
| 96 | In many cases, `deleteat!` simply `deleteat!`s [`full_cache`](@ref) variables and then | ||
| 97 | calls this function. This finer control is required for some `AbstractArray` | ||
| 98 | operations. | ||
| 99 | """ | ||
| 100 | function deleteat_non_user_cache!(i::DEIntegrator, idxs) | ||
| 101 | error("deleteat_non_user_cache!: method has not been implemented for the integrator") | ||
| 102 | end | ||
| 103 | |||
| 104 | """ | ||
| 105 | addat_non_user_cache!(i::DEIntegrator,idxs) | ||
| 106 | |||
| 107 | [`addat!`](@ref)s the non-user facing caches at indices `idxs`. This includes resizing Jacobian caches. | ||
| 108 | |||
| 109 | !!! note | ||
| 110 | |||
| 111 | In many cases, `addat!` simply `addat!`s [`full_cache`](@ref) variables and then | ||
| 112 | calls this function. This finer control is required for some `AbstractArray` | ||
| 113 | operations. | ||
| 114 | """ | ||
| 115 | function addat_non_user_cache!(i::DEIntegrator, idxs) | ||
| 116 | error("addat_non_user_cache!: method has not been implemented for the integrator") | ||
| 117 | end | ||
| 118 | |||
| 119 | """ | ||
| 120 | terminate!(i::DEIntegrator[, retcode = :Terminated]) | ||
| 121 | |||
| 122 | Terminates the integrator by emptying `tstops`. This can be used in events and callbacks to immediately | ||
| 123 | end the solution process. Optionally, `retcode` may be specified (see: [Return Codes (RetCodes)](@ref retcodes)). | ||
| 124 | """ | ||
| 125 | function terminate!(i::DEIntegrator) | ||
| 126 | error("terminate!: method has not been implemented for the integrator") | ||
| 127 | end | ||
| 128 | |||
| 129 | """ | ||
| 130 | get_du(i::DEIntegrator) | ||
| 131 | |||
| 132 | Returns the derivative at `t`. | ||
| 133 | """ | ||
| 134 | function get_du(i::DEIntegrator) | ||
| 135 | error("get_du: method has not been implemented for the integrator") | ||
| 136 | end | ||
| 137 | |||
| 138 | """ | ||
| 139 | get_du!(out,i::DEIntegrator) | ||
| 140 | |||
| 141 | Write the current derivative at `t` into `out`. | ||
| 142 | """ | ||
| 143 | function get_du!(out, i::DEIntegrator) | ||
| 144 | error("get_du: method has not been implemented for the integrator") | ||
| 145 | end | ||
| 146 | function get_dt(i::DEIntegrator) | ||
| 147 | error("get_dt: method has not been implemented for the integrator") | ||
| 148 | end | ||
| 149 | |||
| 150 | """ | ||
| 151 | get_proposed_dt(i::DEIntegrator) | ||
| 152 | |||
| 153 | Gets the proposed `dt` for the next timestep. | ||
| 154 | """ | ||
| 155 | function get_proposed_dt(i::DEIntegrator) | ||
| 156 | error("get_proposed_dt: method has not been implemented for the integrator") | ||
| 157 | end | ||
| 158 | |||
| 159 | """ | ||
| 160 | set_proposed_dt(i::DEIntegrator,dt) | ||
| 161 | set_proposed_dt(i::DEIntegrator,i2::DEIntegrator) | ||
| 162 | |||
| 163 | Sets the proposed `dt` for the next timestep. If the second argument isa `DEIntegrator`, then it sets the timestepping of | ||
| 164 | the first argument to match that of the second one. Note that due to PI control and step acceleration, this is more than matching | ||
| 165 | the factors in most cases. | ||
| 166 | """ | ||
| 167 | function set_proposed_dt!(i::DEIntegrator, dt) | ||
| 168 | error("set_proposed_dt!: method has not been implemented for the integrator") | ||
| 169 | end | ||
| 170 | |||
| 171 | """ | ||
| 172 | savevalues!(integrator::DEIntegrator, | ||
| 173 | force_save=false) -> Tuple{Bool, Bool} | ||
| 174 | |||
| 175 | Try to save the state and time variables at the current time point, or the | ||
| 176 | `saveat` point by using interpolation when appropriate. It returns a tuple that | ||
| 177 | is `(saved, savedexactly)`. If `savevalues!` saved value, then `saved` is true, | ||
| 178 | and if `savevalues!` saved at the current time point, then `savedexactly` is | ||
| 179 | true. | ||
| 180 | |||
| 181 | The saving priority/order is as follows: | ||
| 182 | |||
| 183 | - `save_on` | ||
| 184 | |||
| 185 | + `saveat` | ||
| 186 | + `force_save` | ||
| 187 | + `save_everystep` | ||
| 188 | """ | ||
| 189 | function savevalues!(i::DEIntegrator) | ||
| 190 | error("savevalues!: method has not been implemented for the integrator") | ||
| 191 | end | ||
| 192 | |||
| 193 | """ | ||
| 194 | u_modified!(i::DEIntegrator,bool) | ||
| 195 | |||
| 196 | Sets `bool` which states whether a change to `u` occurred, allowing the solver to handle the discontinuity. By default, | ||
| 197 | this is assumed to be true if a callback is used. This will result in the re-calculation of the derivative at | ||
| 198 | `t+dt`, which is not necessary if the algorithm is FSAL and `u` does not experience a discontinuous change at the | ||
| 199 | end of the interval. Thus, if `u` is unmodified in a callback, a single call to the derivative calculation can be | ||
| 200 | eliminated by `u_modified!(integrator,false)`. | ||
| 201 | """ | ||
| 202 | function u_modified!(i::DEIntegrator, bool) | ||
| 203 | error("u_modified!: method has not been implemented for the integrator") | ||
| 204 | end | ||
| 205 | |||
| 206 | """ | ||
| 207 | add_tstop!(i::DEIntegrator,t) | ||
| 208 | |||
| 209 | Adds a `tstop` at time `t`. | ||
| 210 | """ | ||
| 211 | function add_tstop!(i::DEIntegrator, t) | ||
| 212 | error("add_tstop!: method has not been implemented for the integrator") | ||
| 213 | end | ||
| 214 | |||
| 215 | """ | ||
| 216 | has_tstop(i::DEIntegrator) | ||
| 217 | |||
| 218 | Checks if integrator has any stopping times defined. | ||
| 219 | """ | ||
| 220 | function has_tstop(i::DEIntegrator) | ||
| 221 | error("has_tstop: method has not been implemented for the integrator") | ||
| 222 | end | ||
| 223 | |||
| 224 | """ | ||
| 225 | first_tstop(i::DEIntegrator) | ||
| 226 | |||
| 227 | Gets the first stopping time of the integrator. | ||
| 228 | """ | ||
| 229 | function first_tstop(i::DEIntegrator) | ||
| 230 | error("first_tstop: method has not been implemented for the integrator") | ||
| 231 | end | ||
| 232 | |||
| 233 | """ | ||
| 234 | pop_tstop!(i::DEIntegrator) | ||
| 235 | |||
| 236 | Pops the last stopping time from the integrator. | ||
| 237 | """ | ||
| 238 | function pop_tstop!(i::DEIntegrator) | ||
| 239 | error("pop_tstop!: method has not been implemented for the integrator") | ||
| 240 | end | ||
| 241 | |||
| 242 | """ | ||
| 243 | add_saveat!(i::DEIntegrator,t) | ||
| 244 | |||
| 245 | Adds a `saveat` time point at `t`. | ||
| 246 | """ | ||
| 247 | function add_saveat!(i::DEIntegrator, t) | ||
| 248 | error("add_saveat!: method has not been implemented for the integrator") | ||
| 249 | end | ||
| 250 | |||
| 251 | function set_abstol!(i::DEIntegrator, t) | ||
| 252 | error("set_abstol!: method has not been implemented for the integrator") | ||
| 253 | end | ||
| 254 | function set_reltol!(i::DEIntegrator, t) | ||
| 255 | error("set_reltol!: method has not been implemented for the integrator") | ||
| 256 | end | ||
| 257 | |||
| 258 | """ | ||
| 259 | reinit!(integrator::DEIntegrator,args...; kwargs...) | ||
| 260 | |||
| 261 | The reinit function lets you restart the integration at a new value. | ||
| 262 | |||
| 263 | # Arguments | ||
| 264 | |||
| 265 | - `u0`: Value of `u` to start at. Default value is `integrator.sol.prob.u0` | ||
| 266 | |||
| 267 | # Keyword Arguments | ||
| 268 | |||
| 269 | - `t0`: Starting timepoint. Default value is `integrator.sol.prob.tspan[1]` | ||
| 270 | - `tf`: Ending timepoint. Default value is `integrator.sol.prob.tspan[2]` | ||
| 271 | - `erase_sol=true`: Whether to start with no other values in the solution, or keep the previous solution. | ||
| 272 | - `tstops`, `d_discontinuities`, & `saveat`: Cache where these are stored. Default is the original cache. | ||
| 273 | - `reset_dt`: Set whether to reset the current value of `dt` using the automatic `dt` determination algorithm. Default is | ||
| 274 | `(integrator.dtcache == zero(integrator.dt)) && integrator.opts.adaptive` | ||
| 275 | - `reinit_callbacks`: Set whether to run the callback initializations again (and `initialize_save` is for that). Default is `true`. | ||
| 276 | - `reinit_cache`: Set whether to re-run the cache initialization function (i.e. resetting FSAL, not allocating vectors) | ||
| 277 | which should usually be true for correctness. Default is `true`. | ||
| 278 | |||
| 279 | Additionally, once can access [`auto_dt_reset!`](@ref) which will run the auto `dt` initialization algorithm. | ||
| 280 | """ | ||
| 281 | function reinit!(integrator::DEIntegrator, args...; kwargs...) | ||
| 282 | error("reinit!: method has not been implemented for the integrator") | ||
| 283 | end | ||
| 284 | |||
| 285 | """ | ||
| 286 | initialize_dae!(integrator::DEIntegrator,initializealg = integrator.initializealg) | ||
| 287 | |||
| 288 | Runs the DAE initialization to find a consistent state vector. The optional | ||
| 289 | argument `initializealg` can be used to specify a different initialization | ||
| 290 | algorithm to use. | ||
| 291 | """ | ||
| 292 | function initialize_dae!(integrator::DEIntegrator) | ||
| 293 | error("initialize_dae!: method has not been implemented for the integrator") | ||
| 294 | end | ||
| 295 | function initialize_dae!(integrator::DEIntegrator, initializealg) | ||
| 296 | if !(initializealg isa NoInit) | ||
| 297 | error("initialize_dae!: $(typeof(initializealg)) method has not been implemented for the integrator") | ||
| 298 | end | ||
| 299 | end | ||
| 300 | |||
| 301 | """ | ||
| 302 | auto_dt_reset!(integrator::DEIntegrator) | ||
| 303 | |||
| 304 | Run the auto `dt` initialization algorithm. | ||
| 305 | """ | ||
| 306 | function auto_dt_reset!(integrator::DEIntegrator) | ||
| 307 | error("auto_dt_reset!: method has not been implemented for the integrator") | ||
| 308 | end | ||
| 309 | |||
| 310 | """ | ||
| 311 | change_t_via_interpolation!(integrator::DEIntegrator,t,modify_save_endpoint=Val{false}) | ||
| 312 | |||
| 313 | Modifies the current `t` and changes all of the corresponding values using the local interpolation. If the current solution | ||
| 314 | has already been saved, one can provide the optional value `modify_save_endpoint` to also modify the endpoint of `sol` in the | ||
| 315 | same manner. | ||
| 316 | """ | ||
| 317 | function change_t_via_interpolation!(i::DEIntegrator, args...) | ||
| 318 | error("change_t_via_interpolation!: method has not been implemented for the integrator") | ||
| 319 | end | ||
| 320 | |||
| 321 | addsteps!(i::DEIntegrator, args...) = nothing | ||
| 322 | |||
| 323 | """ | ||
| 324 | reeval_internals_due_to_modification!(integrator::DEIntegrator, continuous_modification=true) | ||
| 325 | |||
| 326 | Update DE integrator after changes by callbacks. | ||
| 327 | For DAEs (either implicit or semi-explicit), this requires re-solving alebraic variables. | ||
| 328 | If continuous_modification is true (or unspecified), this should also recalculate interpolation data. | ||
| 329 | Otherwise the integrator is allowed to skip recalculating the interpolation. | ||
| 330 | """ | ||
| 331 | function reeval_internals_due_to_modification!( | ||
| 332 | integrator::DEIntegrator, continuous_modification) | ||
| 333 | reeval_internals_due_to_modification!(integrator::DEIntegrator) | ||
| 334 | end | ||
| 335 | reeval_internals_due_to_modification!(integrator::DEIntegrator) = nothing | ||
| 336 | |||
| 337 | """ | ||
| 338 | set_t!(integrator::DEIntegrator, t) | ||
| 339 | |||
| 340 | Set current time point of the `integrator` to `t`. | ||
| 341 | """ | ||
| 342 | function set_t!(integrator::DEIntegrator, t) | ||
| 343 | error("set_t!: method has not been implemented for the integrator") | ||
| 344 | end | ||
| 345 | |||
| 346 | """ | ||
| 347 | set_u!(integrator::DEIntegrator, u) | ||
| 348 | set_u!(integrator::DEIntegrator, sym, val) | ||
| 349 | |||
| 350 | Set current state of the `integrator` to `u`. Alternatively, set the state of variable | ||
| 351 | `sym` to value `val`. | ||
| 352 | """ | ||
| 353 | function set_u! end | ||
| 354 | |||
| 355 | function set_u!(integrator::DEIntegrator, u) | ||
| 356 | error("set_u!: method has not been implemented for the integrator") | ||
| 357 | end | ||
| 358 | |||
| 359 | function set_u!(integrator::DEIntegrator, sym, val) | ||
| 360 | # So any error checking happens to ensure we actually _can_ set state | ||
| 361 | set_u!(integrator, integrator.u) | ||
| 362 | |||
| 363 | if symbolic_type(sym) == NotSymbolic() | ||
| 364 | error("sym must be a symbol") | ||
| 365 | end | ||
| 366 | i = variable_index(integrator, sym) | ||
| 367 | |||
| 368 | if isnothing(i) | ||
| 369 | error("$sym is not a state variable") | ||
| 370 | end | ||
| 371 | |||
| 372 | integrator.u[i] = val | ||
| 373 | u_modified!(integrator, true) | ||
| 374 | end | ||
| 375 | |||
| 376 | """ | ||
| 377 | set_ut!(integrator::DEIntegrator, u, t) | ||
| 378 | |||
| 379 | Set current state of the `integrator` to `u` and `t` | ||
| 380 | """ | ||
| 381 | function set_ut!(integrator::DEIntegrator, u, t) | ||
| 382 | set_u!(integrator, u) | ||
| 383 | set_t!(integrator, t) | ||
| 384 | end | ||
| 385 | |||
| 386 | ### Addat isn't a real thing. Let's make it a real thing Gretchen | ||
| 387 | |||
| 388 | function addat!(a::AbstractArray, idxs, val = nothing) | ||
| 389 | if val === nothing | ||
| 390 | resize!(a, length(a) + length(idxs)) | ||
| 391 | else | ||
| 392 | error("real addat! on arrays isn't supported yet") | ||
| 393 | #= | ||
| 394 | flip_range = last(idxs):-1:idxs.start | ||
| 395 | @show idxs,flip_range | ||
| 396 | splice!(a,flip_range,val) | ||
| 397 | =# | ||
| 398 | end | ||
| 399 | end | ||
| 400 | |||
| 401 | ### Indexing | ||
| 402 | function getsyms(integrator::DEIntegrator) | ||
| 403 | syms = variable_symbols(integrator) | ||
| 404 | if isempty(syms) | ||
| 405 | syms = keys(integrator.u) | ||
| 406 | end | ||
| 407 | return syms | ||
| 408 | end | ||
| 409 | |||
| 410 | function getindepsym(integrator::DEIntegrator) | ||
| 411 | syms = independent_variable_symbols(integrator) | ||
| 412 | if isempty(syms) | ||
| 413 | return nothing | ||
| 414 | end | ||
| 415 | return syms | ||
| 416 | end | ||
| 417 | |||
| 418 | function getparamsyms(integrator::DEIntegrator) | ||
| 419 | psyms = parameter_symbols(integrator) | ||
| 420 | if isempty(psyms) | ||
| 421 | return nothing | ||
| 422 | end | ||
| 423 | return psyms | ||
| 424 | end | ||
| 425 | |||
| 426 | function getobserved(integrator::DEIntegrator) | ||
| 427 | if has_observed(integrator.f) | ||
| 428 | return integrator.f.observed | ||
| 429 | else | ||
| 430 | return DEFAULT_OBSERVED | ||
| 431 | end | ||
| 432 | end | ||
| 433 | |||
| 434 | function sym_to_index(sym, integrator::DEIntegrator) | ||
| 435 | idx = variable_index(integrator, sym) | ||
| 436 | if idx === nothing | ||
| 437 | idx = findfirst(isequal(sym), keys(integrator.u)) | ||
| 438 | end | ||
| 439 | return idx | ||
| 440 | end | ||
| 441 | |||
| 442 | # SymbolicIndexingInterface | ||
| 443 | SymbolicIndexingInterface.symbolic_container(A::DEIntegrator) = A.f | ||
| 444 | SymbolicIndexingInterface.parameter_values(A::DEIntegrator) = A.p | ||
| 445 | SymbolicIndexingInterface.state_values(A::DEIntegrator) = A.u | ||
| 446 | SymbolicIndexingInterface.current_time(A::DEIntegrator) = A.t | ||
| 447 | function SymbolicIndexingInterface.set_state!(A::DEIntegrator, val, idx) | ||
| 448 | # So any error checking happens to ensure we actually _can_ set state | ||
| 449 | set_u!(A, A.u) | ||
| 450 | A.u[idx] = val | ||
| 451 | u_modified!(A, true) | ||
| 452 | end | ||
| 453 | |||
| 454 | SymbolicIndexingInterface.is_time_dependent(::DEIntegrator) = true | ||
| 455 | |||
| 456 | # TODO make this nontrivial once dynamic state selection works | ||
| 457 | SymbolicIndexingInterface.constant_structure(::DEIntegrator) = true | ||
| 458 | |||
| 459 | function Base.getproperty(A::DEIntegrator, sym::Symbol) | ||
| 460 | if sym === :destats && hasfield(typeof(A), :stats) | ||
| 461 | @warn "destats has been deprecated for stats" | ||
| 462 | getfield(A, :stats) | ||
| 463 | elseif sym === :ps | ||
| 464 | return ParameterIndexingProxy(A) | ||
| 465 | else | ||
| 466 | return getfield(A, sym) | ||
| 467 | end | ||
| 468 | end | ||
| 469 | |||
| 470 | Base.@propagate_inbounds function _getindex(A::DEIntegrator, | ||
| 471 | ::NotSymbolic, | ||
| 472 | I::Union{Int, AbstractArray{Int}, | ||
| 473 | CartesianIndex, Colon, BitArray, | ||
| 474 | AbstractArray{Bool}}...) | ||
| 475 | A.u[I...] | ||
| 476 | end | ||
| 477 | |||
| 478 | Base.@propagate_inbounds function _getindex(A::DEIntegrator, ::ScalarSymbolic, sym) | ||
| 479 | if is_variable(A, sym) | ||
| 480 | return A[variable_index(A, sym)] | ||
| 481 | elseif is_parameter(A, sym) | ||
| 482 | error("Indexing with parameters is deprecated. Use `getp(sys, $sym)(integrator)` for parameter indexing") | ||
| 483 | elseif is_independent_variable(A, sym) | ||
| 484 | return A.t | ||
| 485 | elseif is_observed(A, sym) | ||
| 486 | return SymbolicIndexingInterface.observed(A, sym)(A.u, A.p, A.t) | ||
| 487 | else | ||
| 488 | error("Tried to index integrator with a Symbol that was not found in the system.") | ||
| 489 | end | ||
| 490 | end | ||
| 491 | |||
| 492 | Base.@propagate_inbounds function _getindex(A::DEIntegrator, ::ArraySymbolic, sym) | ||
| 493 | return A[collect(sym)] | ||
| 494 | end | ||
| 495 | |||
| 496 | Base.@propagate_inbounds function _getindex( | ||
| 497 | A::DEIntegrator, ::ScalarSymbolic, sym::Union{Tuple, AbstractArray}) | ||
| 498 | return getindex.((A,), sym) | ||
| 499 | end | ||
| 500 | |||
| 501 | Base.@propagate_inbounds function Base.getindex(A::DEIntegrator, sym) | ||
| 502 | symtype = symbolic_type(sym) | ||
| 503 | elsymtype = symbolic_type(eltype(sym)) | ||
| 504 | |||
| 505 | if symtype != NotSymbolic() | ||
| 506 | return _getindex(A, symtype, sym) | ||
| 507 | else | ||
| 508 | return _getindex(A, elsymtype, sym) | ||
| 509 | end | ||
| 510 | end | ||
| 511 | |||
| 512 | Base.@propagate_inbounds function Base.getindex( | ||
| 513 | A::DEIntegrator, ::SymbolicIndexingInterface.SolvedVariables) | ||
| 514 | return getindex(A, variable_symbols(A)) | ||
| 515 | end | ||
| 516 | |||
| 517 | Base.@propagate_inbounds function Base.getindex( | ||
| 518 | A::DEIntegrator, ::SymbolicIndexingInterface.AllVariables) | ||
| 519 | return getindex(A, all_variable_symbols(A)) | ||
| 520 | end | ||
| 521 | |||
| 522 | function observed(A::DEIntegrator, sym) | ||
| 523 | getobserved(A)(sym, A.u, A.p, A.t) | ||
| 524 | end | ||
| 525 | |||
| 526 | function Base.setindex!(A::DEIntegrator, val, sym) | ||
| 527 | has_sys(A.f) || | ||
| 528 | error("Invalid indexing of integrator: Integrator does not support indexing without a system") | ||
| 529 | if symbolic_type(sym) == ScalarSymbolic() | ||
| 530 | if is_variable(A, sym) | ||
| 531 | A.u[variable_index(A, sym)] = val | ||
| 532 | u_modified!(A, true) | ||
| 533 | elseif is_parameter(A, sym) | ||
| 534 | error("Parameter indexing is deprecated. Use `setp(sys, $sym)(integrator, $val)` to set parameter value.") | ||
| 535 | else | ||
| 536 | error("Invalid indexing of integrator: $sym is not a state or parameter, it may be an observed variable.") | ||
| 537 | end | ||
| 538 | return A | ||
| 539 | elseif symbolic_type(sym) == ArraySymbolic() | ||
| 540 | setindex!.((A,), val, collect(sym)) | ||
| 541 | return A | ||
| 542 | else | ||
| 543 | sym isa AbstractArray || error("Invalid indexing of integrator") | ||
| 544 | setindex!.((A,), val, sym) | ||
| 545 | return A | ||
| 546 | end | ||
| 547 | end | ||
| 548 | |||
| 549 | ### Integrator traits | ||
| 550 | |||
| 551 | has_reinit(i::DEIntegrator) = false | ||
| 552 | |||
| 553 | ### Display | ||
| 554 | |||
| 555 | function Base.summary(io::IO, I::DEIntegrator) | ||
| 556 | type_color, no_color = get_colorizers(io) | ||
| 557 | print(io, | ||
| 558 | type_color, nameof(typeof(I)), | ||
| 559 | no_color, " with uType ", | ||
| 560 | type_color, typeof(I.u), | ||
| 561 | no_color, " and tType ", | ||
| 562 | type_color, typeof(I.t), | ||
| 563 | no_color) | ||
| 564 | end | ||
| 565 | function Base.show(io::IO, A::DEIntegrator) | ||
| 566 | println(io, string("t: ", A.t)) | ||
| 567 | print(io, "u: ") | ||
| 568 | show(io, A.u) | ||
| 569 | end | ||
| 570 | function Base.show(io::IO, m::MIME"text/plain", A::DEIntegrator) | ||
| 571 | println(io, string("t: ", A.t)) | ||
| 572 | print(io, "u: ") | ||
| 573 | show(io, m, A.u) | ||
| 574 | end | ||
| 575 | |||
| 576 | ### Error check (retcode) | ||
| 577 | |||
| 578 | last_step_failed(integrator::DEIntegrator) = false | ||
| 579 | |||
| 580 | """ | ||
| 581 | check_error(integrator) | ||
| 582 | |||
| 583 | Check state of `integrator` and return one of the | ||
| 584 | [Return Codes](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#retcodes) | ||
| 585 | """ | ||
| 586 | function check_error(integrator::DEIntegrator) | ||
| 587 | if integrator.sol.retcode ∉ (ReturnCode.Success, ReturnCode.Default) | ||
| 588 | return integrator.sol.retcode | ||
| 589 | end | ||
| 590 | # This implementation is intended to be used for ODEIntegrator and | ||
| 591 | # SDEIntegrator. | ||
| 592 | if isnan(integrator.dt) | ||
| 593 | if integrator.opts.verbose | ||
| 594 | @warn("NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.") | ||
| 595 | end | ||
| 596 | return ReturnCode.DtNaN | ||
| 597 | end | ||
| 598 | if integrator.iter > integrator.opts.maxiters | ||
| 599 | if integrator.opts.verbose | ||
| 600 | @warn("Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).") | ||
| 601 | end | ||
| 602 | return ReturnCode.MaxIters | ||
| 603 | end | ||
| 604 | |||
| 605 | # The last part: | ||
| 606 | # If you are close to the end, don't exit: let the user hit the end! | ||
| 607 | # However, if we try that and the step fails, exit instead of infinite loop | ||
| 608 | if !integrator.opts.force_dtmin && integrator.opts.adaptive && | ||
| 609 | abs(integrator.dt) <= abs(integrator.opts.dtmin) && | ||
| 610 | (((hasproperty(integrator, :opts) && hasproperty(integrator.opts, :tstops)) ? | ||
| 611 | integrator.t + integrator.dt < integrator.tdir * first(integrator.opts.tstops) : | ||
| 612 | true) || (hasproperty(integrator, :accept_step) && !integrator.accept_step)) | ||
| 613 | if integrator.opts.verbose | ||
| 614 | if isdefined(integrator, :EEst) | ||
| 615 | EEst = ", and step error estimate = $(integrator.EEst)" | ||
| 616 | else | ||
| 617 | EEst = "" | ||
| 618 | end | ||
| 619 | @warn("dt($(integrator.dt)) <= dtmin($(integrator.opts.dtmin)) at t=$(integrator.t)$EEst. Aborting. There is either an error in your model specification or the true solution is unstable.") | ||
| 620 | end | ||
| 621 | return ReturnCode.DtLessThanMin | ||
| 622 | end | ||
| 623 | if integrator.opts.unstable_check(integrator.dt, integrator.u, integrator.p, | ||
| 624 | integrator.t) | ||
| 625 | if integrator.opts.verbose | ||
| 626 | @warn("Instability detected. Aborting") | ||
| 627 | end | ||
| 628 | return ReturnCode.Unstable | ||
| 629 | end | ||
| 630 | if last_step_failed(integrator) | ||
| 631 | if integrator.opts.verbose | ||
| 632 | @warn("Newton steps could not converge and algorithm is not adaptive. Use a lower dt.") | ||
| 633 | end | ||
| 634 | return ReturnCode.ConvergenceFailure | ||
| 635 | end | ||
| 636 | return ReturnCode.Success | ||
| 637 | end | ||
| 638 | |||
| 639 | function postamble! end | ||
| 640 | |||
| 641 | """ | ||
| 642 | check_error!(integrator) | ||
| 643 | |||
| 644 | Same as `check_error` but also set solution's return code | ||
| 645 | (`integrator.sol.retcode`) and run `postamble!`. | ||
| 646 | """ | ||
| 647 | function check_error!(integrator::DEIntegrator) | ||
| 648 | code = check_error(integrator) | ||
| 649 | if code != ReturnCode.Success | ||
| 650 | integrator.sol = solution_new_retcode(integrator.sol, code) | ||
| 651 | postamble!(integrator) | ||
| 652 | end | ||
| 653 | return code | ||
| 654 | end | ||
| 655 | |||
| 656 | ### Default Iterator Interface | ||
| 657 | function done(integrator::DEIntegrator) | ||
| 658 | if !(integrator.sol.retcode in (ReturnCode.Default, ReturnCode.Success)) | ||
| 659 | return true | ||
| 660 | elseif isempty(integrator.opts.tstops) | ||
| 661 | postamble!(integrator) | ||
| 662 | return true | ||
| 663 | elseif integrator.just_hit_tstop | ||
| 664 | integrator.just_hit_tstop = false | ||
| 665 | if integrator.opts.stop_at_next_tstop | ||
| 666 | postamble!(integrator) | ||
| 667 | return true | ||
| 668 | end | ||
| 669 | end | ||
| 670 | false | ||
| 671 | end | ||
| 672 | function Base.iterate(integrator::DEIntegrator, state = 0) | ||
| 673 | done(integrator) && return nothing | ||
| 674 | state += 1 | ||
| 675 | step!(integrator) # Iter updated in the step! header | ||
| 676 | # Next is callbacks -> iterator -> top | ||
| 677 | return integrator, state | ||
| 678 | end | ||
| 679 | |||
| 680 | Base.eltype(::Type{T}) where {T <: DEIntegrator} = T | ||
| 681 | Base.IteratorSize(::Type{<:DEIntegrator}) = Base.SizeUnknown() | ||
| 682 | |||
| 683 | ### Other Iterators | ||
| 684 | |||
| 685 | struct IntegratorTuples{I} | ||
| 686 | integrator::I | ||
| 687 | end | ||
| 688 | |||
| 689 | function Base.iterate(tup::IntegratorTuples, state = 0) | ||
| 690 | done(tup.integrator) && return nothing | ||
| 691 | step!(tup.integrator) # Iter updated in the step! header | ||
| 692 | state += 1 | ||
| 693 | # Next is callbacks -> iterator -> top | ||
| 694 | return (tup.integrator.u, tup.integrator.t), state | ||
| 695 | end | ||
| 696 | |||
| 697 | function Base.eltype(::Type{ | ||
| 698 | IntegratorTuples{I}, | ||
| 699 | }) where {U, T, | ||
| 700 | I <: | ||
| 701 | DEIntegrator{<:Any, <:Any, U, T}} | ||
| 702 | Tuple{U, T} | ||
| 703 | end | ||
| 704 | Base.IteratorSize(::Type{<:IntegratorTuples}) = Base.SizeUnknown() | ||
| 705 | |||
| 706 | RecursiveArrayTools.tuples(integrator::DEIntegrator) = IntegratorTuples(integrator) | ||
| 707 | |||
| 708 | """ | ||
| 709 | $(TYPEDEF) | ||
| 710 | """ | ||
| 711 | struct IntegratorIntervals{I} | ||
| 712 | integrator::I | ||
| 713 | end | ||
| 714 | |||
| 715 | function Base.iterate(tup::IntegratorIntervals, state = 0) | ||
| 716 | done(tup.integrator) && return nothing | ||
| 717 | state += 1 | ||
| 718 | step!(tup.integrator) # Iter updated in the step! header | ||
| 719 | # Next is callbacks -> iterator -> top | ||
| 720 | return (tup.integrator.uprev, tup.integrator.tprev, tup.integrator.u, tup.integrator.t), | ||
| 721 | state | ||
| 722 | end | ||
| 723 | |||
| 724 | function Base.eltype(::Type{ | ||
| 725 | IntegratorIntervals{I}, | ||
| 726 | }) where {U, T, | ||
| 727 | I <: | ||
| 728 | DEIntegrator{<:Any, <:Any, U, T | ||
| 729 | }} | ||
| 730 | Tuple{U, T, U, T} | ||
| 731 | end | ||
| 732 | Base.IteratorSize(::Type{<:IntegratorIntervals}) = Base.SizeUnknown() | ||
| 733 | |||
| 734 | intervals(integrator::DEIntegrator) = IntegratorIntervals(integrator) | ||
| 735 | |||
| 736 | struct TimeChoiceIterator{T, T2} | ||
| 737 | integrator::T | ||
| 738 | ts::T2 | ||
| 739 | end | ||
| 740 | |||
| 741 | function Base.iterate(iter::TimeChoiceIterator, state = 1) | ||
| 742 | state > length(iter.ts) && return nothing | ||
| 743 | t = iter.ts[state] | ||
| 744 | integrator = iter.integrator | ||
| 745 | if isinplace(integrator.sol.prob) | ||
| 746 | tmp = first(get_tmp_cache(integrator)) | ||
| 747 | if t == integrator.t | ||
| 748 | tmp .= integrator.u | ||
| 749 | elseif t < integrator.t | ||
| 750 | integrator(tmp, t) | ||
| 751 | else | ||
| 752 | step!(integrator, t - integrator.t) | ||
| 753 | integrator(tmp, t) | ||
| 754 | end | ||
| 755 | return (tmp, t), state + 1 | ||
| 756 | else | ||
| 757 | if t == integrator.t | ||
| 758 | tmp = integrator.u | ||
| 759 | elseif t < integrator.t | ||
| 760 | tmp = integrator(t) | ||
| 761 | else | ||
| 762 | step!(integrator, t - integrator.t) | ||
| 763 | tmp = integrator(t) | ||
| 764 | end | ||
| 765 | return (tmp, t), state + 1 | ||
| 766 | end | ||
| 767 | end | ||
| 768 | |||
| 769 | Base.length(iter::TimeChoiceIterator) = length(iter.ts) | ||
| 770 | |||
| 771 | @recipe function f(integrator::DEIntegrator; | ||
| 772 | denseplot = (integrator.opts.calck || | ||
| 773 | integrator isa AbstractSDEIntegrator) && | ||
| 774 | integrator.iter > 0, | ||
| 775 | plotdensity = 10, | ||
| 776 | plot_analytic = false, vars = nothing, idxs = nothing) | ||
| 777 | if vars !== nothing | ||
| 778 | Base.depwarn( | ||
| 779 | "To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.", | ||
| 780 | :f; force = true) | ||
| 781 | (idxs !== nothing) && | ||
| 782 | error("Simultaneously using keywords vars and idxs is not supported. Please only use idxs.") | ||
| 783 | idxs = vars | ||
| 784 | end | ||
| 785 | |||
| 786 | int_vars = interpret_vars(idxs, integrator.sol) | ||
| 787 | |||
| 788 | if denseplot | ||
| 789 | # Generate the points from the plot from dense function | ||
| 790 | plott = collect(range(integrator.tprev, integrator.t; length = plotdensity)) | ||
| 791 | if plot_analytic | ||
| 792 | plot_analytic_timeseries = [integrator.sol.prob.f.analytic( | ||
| 793 | integrator.sol.prob.u0, | ||
| 794 | integrator.sol.prob.p, | ||
| 795 | t) for t in plott] | ||
| 796 | end | ||
| 797 | else | ||
| 798 | plott = nothing | ||
| 799 | end | ||
| 800 | |||
| 801 | dims = length(int_vars[1]) | ||
| 802 | for var in int_vars | ||
| 803 | @assert length(var) == dims | ||
| 804 | end | ||
| 805 | # Should check that all have the same dims! | ||
| 806 | |||
| 807 | plot_vecs = [] | ||
| 808 | for i in 2:dims | ||
| 809 | push!(plot_vecs, []) | ||
| 810 | end | ||
| 811 | |||
| 812 | labels = String[]# Array{String, 2}(1, length(int_vars)*(1+plot_analytic)) | ||
| 813 | strs = String[] | ||
| 814 | varsyms = variable_symbols(integrator) | ||
| 815 | @show plott | ||
| 816 | |||
| 817 | for x in int_vars | ||
| 818 | for j in 2:dims | ||
| 819 | if denseplot | ||
| 820 | if (x[j] isa Integer && x[j] == 0) || | ||
| 821 | isequal(x[j], getindepsym_defaultt(integrator)) | ||
| 822 | push!(plot_vecs[j - 1], plott) | ||
| 823 | else | ||
| 824 | push!(plot_vecs[j - 1], Vector(integrator(plott; idxs = x[j]))) | ||
| 825 | end | ||
| 826 | else # just get values | ||
| 827 | if x[j] == 0 | ||
| 828 | push!(plot_vecs[j - 1], integrator.t) | ||
| 829 | elseif x[j] == 1 && !(integrator.u isa AbstractArray) | ||
| 830 | push!(plot_vecs[j - 1], integrator.u) | ||
| 831 | else | ||
| 832 | push!(plot_vecs[j - 1], integrator.u[x[j]]) | ||
| 833 | end | ||
| 834 | end | ||
| 835 | |||
| 836 | if !isempty(varsyms) && x[j] isa Integer | ||
| 837 | push!(strs, String(getname(varsyms[x[j]]))) | ||
| 838 | elseif hasname(x[j]) | ||
| 839 | push!(strs, String(getname(x[j]))) | ||
| 840 | else | ||
| 841 | push!(strs, "u[$(x[j])]") | ||
| 842 | end | ||
| 843 | end | ||
| 844 | add_labels!(labels, x, dims, integrator.sol, strs) | ||
| 845 | end | ||
| 846 | |||
| 847 | if plot_analytic | ||
| 848 | for x in int_vars | ||
| 849 | for j in 1:dims | ||
| 850 | if denseplot | ||
| 851 | push!(plot_vecs[j], | ||
| 852 | u_n(plot_timeseries, x[j], sol, plott, plot_timeseries)) | ||
| 853 | else # Just get values | ||
| 854 | if x[j] == 0 | ||
| 855 | push!(plot_vecs[j], integrator.t) | ||
| 856 | elseif x[j] == 1 && !(integrator.u isa AbstractArray) | ||
| 857 | push!(plot_vecs[j], | ||
| 858 | integrator.sol.prob.f(Val{:analytic}, integrator.t, | ||
| 859 | integrator.sol[1])) | ||
| 860 | else | ||
| 861 | push!(plot_vecs[j], | ||
| 862 | integrator.sol.prob.f(Val{:analytic}, integrator.t, | ||
| 863 | integrator.sol[1])[x[j]]) | ||
| 864 | end | ||
| 865 | end | ||
| 866 | end | ||
| 867 | add_labels!(labels, x, dims, integrator.sol, strs) | ||
| 868 | end | ||
| 869 | end | ||
| 870 | |||
| 871 | xflip --> integrator.tdir < 0 | ||
| 872 | |||
| 873 | if denseplot | ||
| 874 | seriestype --> :path | ||
| 875 | else | ||
| 876 | seriestype --> :scatter | ||
| 877 | end | ||
| 878 | |||
| 879 | # Special case labels when idxs = (:x,:y,:z) or (:x) or [:x,:y] ... | ||
| 880 | if idxs isa Tuple && (typeof(idxs[1]) == Symbol && typeof(idxs[2]) == Symbol) | ||
| 881 | xlabel --> idxs[1] | ||
| 882 | ylabel --> idxs[2] | ||
| 883 | if length(idxs) > 2 | ||
| 884 | zlabel --> idxs[3] | ||
| 885 | end | ||
| 886 | end | ||
| 887 | if getindex.(int_vars, 1) == zeros(length(int_vars)) || | ||
| 888 | getindex.(int_vars, 2) == zeros(length(int_vars)) | ||
| 889 | xlabel --> "t" | ||
| 890 | end | ||
| 891 | |||
| 892 | linewidth --> 3 | ||
| 893 | #xtickfont --> font(11) | ||
| 894 | #ytickfont --> font(11) | ||
| 895 | #legendfont --> font(11) | ||
| 896 | #guidefont --> font(11) | ||
| 897 | label --> reshape(labels, 1, length(labels)) | ||
| 898 | (plot_vecs...,) | ||
| 899 | end | ||
| 900 | |||
| 901 | function step!(integ::DEIntegrator, dt, stop_at_tdt = false) | ||
| 902 | (dt * integ.tdir) < 0 * oneunit(dt) && error("Cannot step backward.") | ||
| 903 | t = integ.t | ||
| 904 | next_t = t + dt | ||
| 905 | stop_at_tdt && add_tstop!(integ, next_t) | ||
| 906 | while integ.t * integ.tdir < next_t * integ.tdir | ||
| 907 | 58 (100 %) |
58 (100 %)
samples spent calling
step!
step!(integ)
|
|
| 908 | integ.sol.retcode in (ReturnCode.Default, ReturnCode.Success) || break | ||
| 909 | end | ||
| 910 | end | ||
| 911 | |||
| 912 | has_stats(i::DEIntegrator) = false | ||
| 913 | |||
| 914 | """ | ||
| 915 | is_integrator_adaptive(i::DEIntegrator) | ||
| 916 | |||
| 917 | Checks if the integrator is adaptive | ||
| 918 | """ | ||
| 919 | function isadaptive(integrator::DEIntegrator) | ||
| 920 | isdefined(integrator.opts, :adaptive) ? integrator.opts.adaptive : false | ||
| 921 | end |