diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b32b533a9..38a8099f4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -17,7 +17,7 @@ jobs: strategy: fail-fast: false matrix: - version: ['1.0', '1', 'nightly'] + version: ['1.0', '1.6', '1', 'nightly'] os: - ubuntu-latest - macos-latest diff --git a/.gitignore b/.gitignore index 23aa59dca..067fbd84d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ .DS_Store +Manifest.toml examples/.ipynb_checkpoints examples/JuliaCon2017/.ipynb_checkpoints/ docs/build/ diff --git a/Project.toml b/Project.toml index 7d2fbab3e..c4f242380 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorIntegration" uuid = "92b13dbe-c966-51a2-8445-caca9f8a7d42" repo = "https://github.com/PerezHz/TaylorIntegration.jl.git" -version = "0.8.11" +version = "0.9.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" @@ -9,6 +9,7 @@ Elliptic = "b305315f-e792-5b7a-8f41-49f472929428" Espresso = "6912e4f1-e036-58b0-9138-08d1e6358ea9" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" @@ -34,11 +35,11 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" Elliptic = "b305315f-e792-5b7a-8f41-49f472929428" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["TaylorSeries", "LinearAlgebra", "Test", "Elliptic", "InteractiveUtils", - "Pkg", "DiffEqBase", "StaticArrays"] +test = ["DiffEqBase", "Elliptic", "InteractiveUtils", "LinearAlgebra", "Logging", "Pkg", "StaticArrays", "TaylorSeries", "Test"] diff --git a/docs/src/taylorize.md b/docs/src/taylorize.md index 6df77a451..2705d944b 100644 --- a/docs/src/taylorize.md +++ b/docs/src/taylorize.md @@ -16,17 +16,20 @@ the function where the ODEs of the problem are defined, in order to compute the recurrence relations that are used to construct the Taylor expansion of the solution. This is done for each order of the series in [`TaylorIntegration.jetcoeffs!`](@ref). These computations are -not optimized: they waste memory due to allocations of some temporary +not optimized: they waste memory due to repeated allocations of some temporary arrays, and perform some operations whose result has been previously computed. Here we describe one way to optimize this: The idea is to replace the -default method of [`TaylorIntegration.jetcoeffs!`](@ref) by another (with the -same name) which is called by dispatch, that in principle performs better. -The new method is constructed specifically for the actual function -defining the equations of motion by parsing its expression; the new -function performs in principle *exactly* the same operations, but avoids -the extra allocations and the repetition of some operations. +default method of [`TaylorIntegration.jetcoeffs!`](@ref) by another +method (same function name) which is called by dispatch, and that in principle +performs better. +The new method is constructed specifically for the specific function +defining the equations of motion by parsing its expression. This new +method performs in principle *exactly* the same operations, but avoids +the extra allocations; to do the latter, the macro also creates an *internal* +function `TaylorIntegration._allocate_jetcoeffs!`, which allocates all temporary +`Taylor1` objects as well as the declared `Array{Taylor1,1}`s. ## An example @@ -46,13 +49,14 @@ end # Initial time (t0), final time (tf) and initial condition (q0) t0 = 0.0 -tf = 100.0 +tf = 10000.0 q0 = [1.3, 0.0] # The actual integration t1, x1 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); # warm-up run e1 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); -e1 +all1 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); +e1, all1 ``` We note that the initial number of methods defined for @@ -60,15 +64,23 @@ We note that the initial number of methods defined for ```@example taylorize length(methods(TaylorIntegration.jetcoeffs!)) == 2 # initial value ``` -Using `@taylorize` will increase this number by creating a new method. +Similarly, the number of methods for `TaylorIntegration._allocate_jetcoeffs!` is +also 2. +```@example taylorize +length(methods(TaylorIntegration._allocate_jetcoeffs!)) == 2 # initial value +``` +Using `@taylorize` will increase this number by creating a new method for these +two functions. The macro [`@taylorize`](@ref) is intended to be used in front of the function that implements the equations of motion. The macro does the following: it -first parses the actual function as it is, so the integration can be computed +first parses the function as it is, so the integration can be computed using [`taylorinteg`](@ref) as above, by explicitly using the keyword -argument `parse_eqs=false`. It then creates and evaluates a new method of +argument `parse_eqs=false`; this also declares the function of the ODEs, whose name +is used for parsing. It then creates and evaluates a new method of [`TaylorIntegration.jetcoeffs!`](@ref), which is the specialized method -(through `Val`) on the specific function passed to the macro. +(through `Val`) on the specific function passed to the macro as well as a specialized +`TaylorIntegration._allocate_jetcoeffs!`. ```@example taylorize @taylorize function pendulum!(dx, x, p, t) @@ -85,27 +97,32 @@ println(methods(TaylorIntegration.jetcoeffs!)) ``` We see that there is only one method of `pendulum!`, and -there is a new method of `TaylorIntegration.jetcoeffs!`, whose signature appears -in this documentation as `Val{Main.ex-taylorize.pendulum!}`; it is -an specialized version for the function `pendulum!` (with some extra information +there is a *new* method (three in total) of `TaylorIntegration.jetcoeffs!`, +whose signature appears +in this documentation as `Val{Main.ex-taylorize.pendulum!}`. It is +a specialized version for the function `pendulum!` (with some extra information about the module where the function was created). This method is selected internally if it exists (default), exploiting dispatch, when calling [`taylorinteg`](@ref) or [`lyap_taylorinteg`](@ref); to integrate -using the hard-coded method of [`TaylorIntegration.jetcoeffs!`](@ref) of the +using the hard-coded (default) method of [`TaylorIntegration.jetcoeffs!`](@ref) of the integration above, the keyword argument `parse_eqs` has to be set to `false`. +Similarly, one can check that there exists a new method of +`TaylorIntegration._allocate_jetcoeffs!` Now we carry out the integration using the specialized method; note that we -use the same instruction as above. +use the same instruction as above; the default value for the keyword argument `parse_eqs` +is `true`, so we may omit it. ```@example taylorize t2, x2 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); # warm-up run e2 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); -e2 +all2 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500); +e2, all2 ``` -We note the difference in the performance: +We note the difference in the performance and allocations: ```@example taylorize -e1/e2 +e1/e2, all1/all2 ``` We can check that both integrations yield the same results. @@ -121,8 +138,8 @@ setting it to `false` causes the standard method to be used. taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500, parse_eqs=false); # warm-up run e3 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500, parse_eqs=false); - -e1/e3 +all3 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=1500, parse_eqs=false); +e1/e3, all1/all3 ``` We now illustrate the possibility of exploiting the macro @@ -144,13 +161,13 @@ The speed-up obtained comes from the design of the new (specialized) method of allocations and some repeated computations. This is achieved by knowing the specific AST of the function of the ODEs integrated, which is walked through and *translated* into the actual implementation, where some -required auxiliary arrays are created and the low-level functions defined in +required auxiliary arrays are created, and the low-level functions defined in `TaylorSeries.jl` are used. For this, we heavily rely on [`Espresso.jl`](https://github.com/dfdx/Espresso.jl) and some metaprogramming; we thank Andrei Zhabinski for his help and comments. -The new `jetcoeffs!` method can be inspected by constructing the expression -corresponding to the function, and using +The new `TaylorIntegration.jetcoeffs!` and `TaylorIntegration._allocate_jetcoeffs!` method can be inspected by +constructing the expression corresponding to the function, and using [`TaylorIntegration._make_parsed_jetcoeffs`](@ref): ```@example taylorize @@ -160,12 +177,13 @@ ex = :(function pendulum!(dx, x, p, t) return dx end) -new_ex = TaylorIntegration._make_parsed_jetcoeffs(ex) +new_ex1, new_ex2 = TaylorIntegration._make_parsed_jetcoeffs(ex) ``` -This function has a similar structure as the hard-coded method of +The first function has a similar structure as the hard-coded method of `TaylorIntegration.jetcoeffs!`, but uses low-level functions in `TaylorSeries` -(e.g., `sincos!` above) and explicitly allocates the necessary temporary arrays. +(e.g., `sincos!` above); all temporary `Taylor1` objects as well as declared +arrays are allocated once by `TaylorIntegration._allocate_jetcoeffs!`. More complex functions quickly become very difficult to read. Note that, if necessary, one can further optimize `new_ex` manually. @@ -174,10 +192,10 @@ if necessary, one can further optimize `new_ex` manually. The construction of the internal function obtained by using [`@taylorize`](@ref) is somewhat complicated and limited. Here we -list some limitations and advice. +list some limitations and provide some advice. - It is best to have expressions which involve two arguments at most, which - requires parentheses be appropriately used: For example, `res = a+b+c` should + requires using parentheses appropriately: For example, `res = a+b+c` should be written as `res = (a+b)+c`. - Updating operators such as `+=`, `*=`, etc., are not supported. For @@ -200,13 +218,13 @@ list some limitations and advice. `DifferentialEquations` may not be able to exploit it. - `if-else` blocks are recognized in their long form, but short-circuit - conditional operators (`&&` and `||`) are not. When comparands are subject to + conditional operators (`&&` and `||`) are not. When comparing to a Taylor expansion, use operators such as `iszero` for `if-else` tests rather than comparing against numeric literals. - Input and output lengths should be determined at the time of `@taylorize` application, not at runtime. Do not use the length of the input as an - implicit indicator of whether or not to write all elements of the output. If + implicit indicator of whether to write all elements of the output. If conditional output of auxiliary equations is desired, use explicit methods, such as through parameters or by setting auxiliary `t0` vector elements to zero, and assigning unwanted auxiliary outputs zero. @@ -216,16 +234,16 @@ list some limitations and advice. heuristics used, especially for vectors, may not work for all cases. - Use `local` for internal parameters (simple constant values); this improves - performance. Do not use it if the variable is Taylor expanded. + performance. Do not use it if the variable is Taylor expanded during the integration. - `@taylorize` supports multi-threading via `Threads.@threads`. **WARNING**: this feature is experimental. Since thread-safety depends on the definition of each ODE, we cannot guarantee the resulting code to be thread-safe in advance. The user should check the resulting code to ensure that it is indeed thread-safe. For more information about multi-threading, the reader is - referred to the [Julia documentation](https://docs.julialang.org/en/v1/manual/parallel-computing#man-multithreading-1). + referred to the [Julia documentation](https://docs.julialang.org/en/v1/manual/multi-threading/#man-multithreading). It is recommended to skim `test/taylorize.jl`, which implements different -cases. +cases and highlights cases where the macro doesn't work and how to solve the problem. Please report any problems you may encounter. diff --git a/src/common.jl b/src/common.jl index 7a6ce84a7..6558450a8 100644 --- a/src/common.jl +++ b/src/common.jl @@ -50,25 +50,32 @@ struct TaylorMethodCache{uType, rateType, tTType, uTType} <: OrdinaryDiffEqMutab duT::uTType uauxT::uTType parse_eqs::Ref{Bool} + tmpTaylor::uTType + arrTaylor::Vector{uTType} end full_cache(c::TaylorMethodCache) = begin - tuple(c.u, c.uprev, c.tmp, c.k, c.fsalfirst, c.tT, c.uT, c.duT, c.uauxT, c.parse_eqs) + tuple(c.u, c.uprev, c.tmp, c.k, c.fsalfirst, c.tT, c.uT, c.duT, c.uauxT, c.parse_eqs, + c.tmpTaylor, c.arrTaylor) end struct TaylorMethodConstantCache{uTType} <: OrdinaryDiffEqConstantCache uT::uTType parse_eqs::Ref{Bool} + tmpTaylor::Vector{uTType} + arrTaylor::Vector{Vector{uTType}} end function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck,::Val{true}) - tT = Taylor1(typeof(t), alg.order) + order = alg.order + tT = Taylor1(typeof(t), order) tT[0] = t - uT = Taylor1.(u, tT.order) - duT = zero.(Taylor1.(u, tT.order)) + uT = Taylor1.(u, order) + duT = zero.(Taylor1.(u, order)) uauxT = similar(uT) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p) TaylorMethodCache( u, uprev, @@ -79,21 +86,60 @@ function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits, uT, duT, uauxT, - Ref(alg.parse_eqs) + Ref(parse_eqs), + tmpTaylor, + arrTaylor ) end -alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits, - uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck, - ::Val{false}) = TaylorMethodConstantCache(Taylor1(u, alg.order), Ref(alg.parse_eqs)) +# This method is used for DynamicalODEFunction's (`parse_eqs=false`): tmpT1 and arrT1 +# must have the proper type to build `TaylorMethodCache` +function alg_cache(alg::TaylorMethod, u::ArrayPartition, rate_prototype, uEltypeNoUnits, + uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, + calck, ::Val{true}) + order = alg.order + tT = Taylor1(typeof(t), order) + tT[0] = t + uT = Taylor1.(u, order) + duT = zero.(Taylor1.(u, order)) + uauxT = similar(uT) + parse_eqs, _, _ = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p) + ## `tmpT1` must have the same type than `uT` and `arrT1` is a `[tmpT1]` + tmpT1 = similar(uT) + arrT1 = similar([tmpT1]) + TaylorMethodCache( + u, + uprev, + similar(u), + zero(rate_prototype), + zero(rate_prototype), + tT, + uT, + duT, + uauxT, + Ref(parse_eqs), + tmpT1, + arrT1 + ) +end + +function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits, + uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{false}) + order = alg.order + tT = Taylor1(typeof(t), order) + tT[0] = t + uT = Taylor1(u, order) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(alg.parse_eqs, f, tT, uT, p) + TaylorMethodConstantCache(Taylor1(u, alg.order), Ref(parse_eqs), tmpTaylor, arrTaylor) +end function initialize!(integrator, c::TaylorMethodConstantCache) @unpack u, t, f, p = integrator tT = Taylor1(typeof(t), integrator.alg.order) tT[0] = t c.uT .= Taylor1(u, tT.order) - c.parse_eqs.x = _determine_parsing!(c.parse_eqs.x, f, tT, c.uT, p) - __jetcoeffs!(Val(c.parse_eqs.x), f, tT, c.uT, p) + __jetcoeffs!(Val(c.parse_eqs.x), f, tT, c.uT, p, c.tmpTaylor, c.arrTaylor) # FSAL stuff integrator.kshortsize = 2 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) @@ -111,7 +157,7 @@ function perform_step!(integrator,cache::TaylorMethodConstantCache) tT[0] = t+dt u = evaluate(cache.uT, dt) cache.uT[0] = u - __jetcoeffs!(Val(cache.parse_eqs.x), f, tT, cache.uT, p) + __jetcoeffs!(Val(cache.parse_eqs.x), f, tT, cache.uT, p, cache.tmpTaylor, cache.arrTaylor) k = f(u, p, t+dt) # For the interpolation, needs k at the updated point integrator.destats.nf += 1 integrator.fsallast = k @@ -123,8 +169,7 @@ end function initialize!(integrator, cache::TaylorMethodCache) @unpack u, t, f, p = integrator @unpack k, fsalfirst, tT, uT, duT, uauxT, parse_eqs = cache - parse_eqs.x = _determine_parsing!(parse_eqs.x, f, tT, uT, duT, p) - __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p) + __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, cache.tmpTaylor, cache.arrTaylor) # FSAL for interpolation integrator.fsalfirst = fsalfirst integrator.fsallast = k @@ -138,19 +183,20 @@ end function perform_step!(integrator, cache::TaylorMethodCache) @unpack t, dt, u, f, p = integrator - @unpack k, tT, uT, duT, uauxT, parse_eqs = cache + @unpack k, tT, uT, duT, uauxT, parse_eqs, tmpTaylor, arrTaylor = cache evaluate!(uT, dt, u) tT[0] = t+dt for i in eachindex(u) @inbounds uT[i][0] = u[i] duT[i].coeffs .= zero(duT[i][0]) end - __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p) + __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, tmpTaylor, arrTaylor) k = constant_term.(duT) # For the interpolation, needs k at the updated point integrator.destats.nf += 1 end -stepsize_controller!(integrator,alg::TaylorMethod) = stepsize(integrator.cache.uT, integrator.opts.abstol) +stepsize_controller!(integrator,alg::TaylorMethod) = + stepsize(integrator.cache.uT, integrator.opts.abstol) step_accept_controller!(integrator, alg::TaylorMethod, q) = q function DiffEqBase.solve( @@ -166,7 +212,9 @@ function DiffEqBase.solve( end f = prob.f + parse_eqs = haskey(kwargs, :parse_eqs) ? kwargs[:parse_eqs] : true # `true` is the default if !isinplace && typeof(prob.u0) <: AbstractArray + ### TODO: allow `parse_eqs=true` for DynamicalODEFunction if prob.f isa DynamicalODEFunction f1! = (dv, v, u, p, t) -> (dv .= prob.f.f1(v, u, p, t); 0) f2! = (du, v, u, p, t) -> (du .= prob.f.f2(v, u, p, t); 0) @@ -181,14 +229,11 @@ function DiffEqBase.solve( _prob = DynamicalODEProblem(f1!, f2!, _u0.x[1], _u0.x[2], prob.tspan, prob.p; prob.kwargs...) else f! = (du, u, p, t) -> (du .= f(u, p, t); 0) - _alg = TaylorMethod(alg.order, parse_eqs = false) + _alg = TaylorMethod(alg.order, parse_eqs = parse_eqs) _prob = ODEProblem(f!, prob.u0, prob.tspan, prob.p; prob.kwargs...) end - elseif haskey(kwargs, :parse_eqs) - _alg = TaylorMethod(alg.order, parse_eqs = kwargs[:parse_eqs]) - _prob = prob else - _alg = TaylorMethod(alg.order) + _alg = TaylorMethod(alg.order, parse_eqs = parse_eqs) _prob = prob end @@ -201,12 +246,12 @@ end # used in continuous callbacks and related methods to update Taylor expansions cache function update_jetcoeffs_cache!(u,f,p,cache::TaylorMethodCache) - @unpack tT, uT, duT, uauxT, parse_eqs = cache + @unpack tT, uT, duT, uauxT, parse_eqs, tmpTaylor, arrTaylor = cache for i in eachindex(u) @inbounds uT[i][0] = u[i] duT[i].coeffs .= zero(duT[i][0]) end - __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p) + __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, tmpTaylor, arrTaylor) return nothing end @@ -231,14 +276,17 @@ function DiffEqBase.addsteps!(k, t, uprev, u, dt, f, p, cache::TaylorMethodCache nothing end -@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x, params) = - jetcoeffs!(f.f, t, x, params) -@inline __jetcoeffs!(::Val{true}, f::ODEFunction, t, x, params) = - jetcoeffs!(Val(f.f), t, x, params) -@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x, dx, xaux, params) = - jetcoeffs!(f.f, t, x, dx, xaux, params) -@inline __jetcoeffs!(::Val{true}, f::ODEFunction, t, x, dx, xaux, params) = - jetcoeffs!(Val(f.f), t, x, dx, params) - -_determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, params) = _determine_parsing!(parse_eqs::Bool, f.f, t, x, params) -_determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, dx, params) = _determine_parsing!(parse_eqs::Bool, f.f, t, x, dx, params) +@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x, params, tmpTaylor, arrTaylor) = + __jetcoeffs!(Val(false), f.f, t, x, params, tmpTaylor, arrTaylor) +@inline __jetcoeffs!(::Val{true}, f::ODEFunction, t, x, params, tmpTaylor, arrTaylor) = + __jetcoeffs!(Val(true), f.f, t, x, params, tmpTaylor, arrTaylor) +@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x, dx, xaux, params, tmpTaylor, arrTaylor) = + __jetcoeffs!(Val(false), f.f, t, x, dx, xaux, params, tmpTaylor, arrTaylor) +@inline __jetcoeffs!(::Val{true}, f::ODEFunction, t, x, dx, xaux, params, tmpTaylor, arrTaylor) = + __jetcoeffs!(Val(true), f.f, t, x, dx, xaux, params, tmpTaylor, arrTaylor) + + +_determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, params) = + _determine_parsing!(parse_eqs, f.f, t, x, params) +_determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, dx, params) = + _determine_parsing!(parse_eqs, f.f, t, x, dx, params) diff --git a/src/explicitode.jl b/src/explicitode.jl index 0342efd51..9c774bd0e 100644 --- a/src/explicitode.jl +++ b/src/explicitode.jl @@ -88,20 +88,28 @@ end # __jetcoeffs """ - __jetcoeffs!(::Val{bool}, f, t, x, params) - __jetcoeffs!(::Val{bool}, f, t, x, dx, xaux, params) + __jetcoeffs!(::Val{bool}, f, t, x, params, tmpTaylor, arrTaylor) + __jetcoeffs!(::Val{bool}, f, t, x, dx, xaux, params, tmpTaylor, arrTaylor) -Chooses a method of [`jetcoeffs!`](@ref) (hard-coded or generated by +Chooses a method of [`jetcoeffs!`](@ref) (hard-coded) or the generated by [`@taylorize`](@ref)) depending on `Val{bool}` (`bool::Bool`). """ -@inline __jetcoeffs!(::Val{false}, f, t, x, params) = +@inline __jetcoeffs!(::Val{false}, f, t, x, params, tmpTaylor, arrTaylor) = jetcoeffs!(f, t, x, params) -@inline __jetcoeffs!(::Val{true}, f, t, x, params) = - jetcoeffs!(Val(f), t, x, params) -@inline __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params) = +@inline __jetcoeffs!(::Val{true}, f, t, x, params, tmpTaylor, arrTaylor) = + jetcoeffs!(Val(f), t, x, params, tmpTaylor, arrTaylor) +@inline __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params, tmpTaylor, arrTaylor) = jetcoeffs!(f, t, x, dx, xaux, params) -@inline __jetcoeffs!(::Val{true}, f, t, x, dx, xaux, params) = - jetcoeffs!(Val(f), t, x, dx, params) +@inline __jetcoeffs!(::Val{true}, f, t, x, dx, xaux, params, tmpTaylor, arrTaylor) = + jetcoeffs!(Val(f), t, x, dx, params, tmpTaylor, arrTaylor) + + +# _allocate_jetcoeffs! fallbacks; needed to define the corresponding parsed method +_allocate_jetcoeffs!(::Taylor1{T}, x::Taylor1{S}, params) where {T,S} = + (typeof(x)[], Vector{typeof(x)}[]) +_allocate_jetcoeffs!(::Taylor1{T}, x::AbstractArray{Taylor1{S}, N}, + ::AbstractArray{Taylor1{S}, N}, params) where {T,S,N} = + (eltype(x)[], Vector{eltype(x)}[]) # stepsize @@ -190,8 +198,8 @@ end #taylorstep @doc doc""" - taylorstep!(f, t, x, t0, order, abstol, params, parse_eqs=true) -> δt - taylorstep!(f!, t, x, dx, xaux, t0, order, abstol, params, parse_eqs=true) -> δt + taylorstep!(f, t, x, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt + taylorstep!(f!, t, x, dx, xaux, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt One-step Taylor integration for the one-dependent variable ODE ``\dot{x}=dx/dt=f(x, p, t)`` with initial conditions ``x(t_0)=x_0``. @@ -215,10 +223,10 @@ created with [`@taylorize`](@ref); the default is `true` (parse the equations). """ function taylorstep!(f, t::Taylor1{T}, x::Taylor1{U}, abstol::T, params, - parse_eqs::Bool=true) where {T<:Real, U<:Number} + tmpTaylor, arrTaylor, parse_eqs::Bool=true) where {T<:Real, U<:Number} # Compute the Taylor coefficients - __jetcoeffs!(Val(parse_eqs), f, t, x, params) + __jetcoeffs!(Val(parse_eqs), f, t, x, params, tmpTaylor, arrTaylor) # Compute the step-size of the integration using `abstol` δt = stepsize(x, abstol) @@ -231,10 +239,10 @@ end function taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, abstol::T, params, - parse_eqs::Bool=true) where {T<:Real, U<:Number} + tmpTaylor, arrTaylor, parse_eqs::Bool=true) where {T<:Real, U<:Number} # Compute the Taylor coefficients - __jetcoeffs!(Val(parse_eqs), f!, t, x, dx, xaux, params) + __jetcoeffs!(Val(parse_eqs), f!, t, x, dx, xaux, params, tmpTaylor, arrTaylor) # Compute the step-size of the integration using `abstol` δt = stepsize(x, abstol) @@ -252,34 +260,44 @@ Check if the parsed method of `jetcoeffs!` exists and check it runs without error. """ function _determine_parsing!(parse_eqs::Bool, f, t, x, params) - parse_eqs = parse_eqs && - !isempty(methodswith(Val{f}, TaylorIntegration.jetcoeffs!)) + parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) && + !isempty(methodswith(Val{f}, _allocate_jetcoeffs!)) + + tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, params) + if parse_eqs try - jetcoeffs!(Val(f), t, x, params) + tmpTaylor, arrTaylor = _allocate_jetcoeffs!(Val(f), t, x, params) + jetcoeffs!(Val(f), t, x, params, tmpTaylor, arrTaylor) catch - @warn("""Unable to use the parsed method of `jetcoeffs!` + @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`, despite of having `parse_eqs=true`, due to some internal error. - Using `parse_eqs = false`""") + Using `parse_eqs = false`.""") parse_eqs = false end end - return parse_eqs + + return parse_eqs, tmpTaylor, arrTaylor end function _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params) - parse_eqs = parse_eqs && - !isempty(methodswith(Val{f}, TaylorIntegration.jetcoeffs!)) + parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) && + !isempty(methodswith(Val{f}, _allocate_jetcoeffs!)) + + tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, dx, params) + if parse_eqs try - jetcoeffs!(Val(f), t, x, dx, params) + tmpTaylor, arrTaylor = _allocate_jetcoeffs!(Val(f), t, x, dx, params) + jetcoeffs!(Val(f), t, x, dx, params, tmpTaylor, arrTaylor) catch - @warn("""Unable to use the parsed method of `jetcoeffs!` + @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`, despite of having `parse_eqs=true`, due to some internal error. - Using `parse_eqs = false`""") + Using `parse_eqs = false`.""") parse_eqs = false end end - return parse_eqs + + return parse_eqs, tmpTaylor, arrTaylor end @@ -369,11 +387,11 @@ function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T, sign_tstep = copysign(1, tmax-t0) # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f, t, x, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params) # Integration while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f, t, x, abstol, params, parse_eqs) # δt is positive! + δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) x0 = evaluate(x, δt) # new initial condition @@ -417,11 +435,12 @@ function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T, sign_tstep = copysign(1, tmax-t0) # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f, t, x, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params) + # Integration while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f, t, x, abstol, params, parse_eqs) # δt is positive! + δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) x0 = evaluate(x, δt) # new initial condition @@ -470,12 +489,12 @@ function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T, sign_tstep = copysign(1, tmax-t0) # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params) # Integration nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive! + δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) evaluate!(x, δt, x0) # new initial condition @@ -530,12 +549,12 @@ function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T, sign_tstep = copysign(1, tmax-t0) # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params) # Integration nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive! + δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) evaluate!(x, δt, x0) # new initial condition @@ -609,13 +628,13 @@ function taylorinteg(f, x0::U, trange::AbstractVector{T}, @inbounds xv[1] = x0 # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f, t, x, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params) # Integration iter = 2 nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f, t, x, abstol, params, parse_eqs)# δt is positive! + δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor, parse_eqs)# δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) x0 = evaluate(x, δt) # new initial condition @@ -683,13 +702,13 @@ function taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, @inbounds xv[:,1] .= q0 # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params) # Integration iter = 2 nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive! + δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) evaluate!(x, δt, x0) # new initial condition diff --git a/src/lyapunovspectrum.jl b/src/lyapunovspectrum.jl index c613b5112..d116ca361 100644 --- a/src/lyapunovspectrum.jl +++ b/src/lyapunovspectrum.jl @@ -15,7 +15,7 @@ Optionally, the user may provide a Jacobian function `jacobianfunc!` to compute function stabilitymatrix!(eqsdiff!, t::Taylor1{T}, x::Vector{Taylor1{U}}, δx::Vector{TaylorN{Taylor1{U}}}, dδx::Vector{TaylorN{Taylor1{U}}}, jac::Matrix{Taylor1{U}}, _δv::Vector{TaylorN{Taylor1{U}}}, params, - jacobianfunc! =nothing) where {T<:Real, U<:Number} + jacobianfunc! = nothing) where {T<:Real, U<:Number} if isa(jacobianfunc!, Nothing) # Set δx equal to current value of x plus 1st-order variations @@ -170,7 +170,7 @@ function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, δx::Array{TaylorN{Taylor1{U}},1}, dδx::Array{TaylorN{Taylor1{U}},1}, jac::Array{Taylor1{U},2}, abstol::T, _δv::Vector{TaylorN{Taylor1{U}}}, - varsaux::Array{Taylor1{U},3}, params, parse_eqs::Bool=true, + varsaux::Array{Taylor1{U},3}, params, tmpTaylor, arrTaylor, parse_eqs::Bool=true, jacobianfunc! =nothing) where {T<:Real, U<:Number} # Dimensions of phase-space: dof @@ -178,7 +178,8 @@ function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dof = length(δx) # Compute the Taylor coefficients associated to trajectory - __jetcoeffs!(Val(parse_eqs), f!, t, view(x, 1:dof), view(dx, 1:dof), view(xaux, 1:dof), params) + __jetcoeffs!(Val(parse_eqs), f!, t, view(x, 1:dof), view(dx, 1:dof), + view(xaux, 1:dof), params, tmpTaylor, arrTaylor) # Compute stability matrix stabilitymatrix!(f!, t, x, δx, dδx, jac, _δv, params, jacobianfunc!) @@ -259,13 +260,14 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, vⱼ = similar(aⱼ) # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f!, t, view(x, 1:dof), view(dx, 1:dof), params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, + view(x, 1:dof), view(dx, 1:dof), params) # Integration nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, - abstol, _δv, varsaux, params, parse_eqs, jacobianfunc!) # δt is positive! + δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, + varsaux, params, tmpTaylor, arrTaylor, parse_eqs, jacobianfunc!) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) evaluate!(x, δt, x0) # Update x0 @@ -355,14 +357,15 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, vⱼ = similar(aⱼ) # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f!, t, view(x, 1:dof), view(dx, 1:dof), params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!( + parse_eqs, f!, t, view(x, 1:dof), view(dx, 1:dof), params) # Integration iter = 2 nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, - abstol, _δv, varsaux, params, parse_eqs, jacobianfunc!) # δt is positive! + δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, + varsaux, params, tmpTaylor, arrTaylor, parse_eqs, jacobianfunc!) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) evaluate!(x, δt, x0) # Update x0 diff --git a/src/parse_eqs.jl b/src/parse_eqs.jl index 086cd016d..a5a6b1d7a 100644 --- a/src/parse_eqs.jl +++ b/src/parse_eqs.jl @@ -5,93 +5,190 @@ using Espresso: subs, simplify, ExGraph, ExH, to_expr, sanitize, genname, find_vars, findex, find_indices, isindexed -# Define some constants to create the new (parsed) functions -# The (irrelevant) `nothing` below is there to have a :block Expr; deleted later -const _HEAD_PARSEDFN_SCALAR = sanitize(:( -function TaylorIntegration.jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T}, __x::Taylor1{_S}, __params) where - {_T<:Real, _S<:Number} +""" + BookKeeping + +Mutable struct that contains all the bookkeeping vectors/dictionaries used within +`_make_parsed_jetcoeffs`: + - `d_indx` : Dictionary mapping new variables (symbols) to old (perhaps indexed) symbols + - `d_assign` : Dictionary with the numeric assignments (that are substituted) + - `d_decl` : Dictionary declared arrays + - `v_newvars` : Symbols of auxiliary indexed vars + - `v_arraydecl`: Symbols which are explicitly declared as Array or Vector + - `v_preamb` : Symbols or Expr used in the preamble (declarations, etc) + - `retvar` : *Guessed* returned variable, which defines the LHS of the ODEs + +""" +mutable struct BookKeeping + d_indx :: Dict{Symbol, Expr} + d_assign :: Dict{Union{Symbol,Expr}, Number} + d_decl :: Dict{Symbol, Expr} + v_newvars :: Vector{Symbol} + v_arraydecl :: Vector{Symbol} + v_preamb :: Vector{Union{Symbol,Expr}} + retvar :: Symbol + + function BookKeeping() + return new(Dict{Symbol, Expr}(), Dict{Union{Symbol,Expr}, Number}(), + Dict{Symbol, Expr}(), Symbol[], Symbol[], Union{Symbol,Expr}[], :nothing) + end +end + +""" +`inbookkeeping(v, bkkeep::BookKeeping)` + +Checks if `v` is declared in `bkkeep`, considering the `d_indx`, `v_newvars` and +`v_arraydecl` fields. + +""" +@inline inbookkeeping(v, bkkeep::BookKeeping) = (v ∈ keys(bkkeep.d_indx) || + v ∈ bkkeep.v_newvars || v ∈ bkkeep.v_arraydecl) - order = __tT.order - nothing -end) + +# Constants to create the structure of the new jetcoeffs! and _allocate_jetcoeffs! methods. +# The (irrelevant) `nothing` below is there to have a `:block` Expr; it is deleted later +const _HEAD_PARSEDFN_SCALAR = sanitize(:( + function TaylorIntegration.jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T}, __x::Taylor1{_S}, + __params, __tmpTaylor, __arrTaylor) where {_T<:Real, _S<:Number} + order = __tT.order + nothing + end) ); const _HEAD_PARSEDFN_VECTOR = sanitize(:( -function TaylorIntegration.jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T}, __x::AbstractArray{Taylor1{_S}, _N}, - __dx::AbstractArray{Taylor1{_S}, _N}, __params) where {_T<:Real, _S<:Number, _N} + function TaylorIntegration.jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T}, + __x::AbstractArray{Taylor1{_S}, _N}, __dx::AbstractArray{Taylor1{_S}, _N}, + __params, __tmpTaylor, __arrTaylor) where {_T<:Real, _S<:Number, _N} + order = __tT.order + nothing + end) +); - order = __tT.order - nothing -end) +const _HEAD_ALLOC_TAYLOR1_SCALAR = sanitize(:( + function TaylorIntegration._allocate_jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T}, + __x::Taylor1{_S}, __params) where {_T<:Real, _S<:Number} + order = __tT.order + nothing + end) +); + +const _HEAD_ALLOC_TAYLOR1_VECTOR = sanitize(:( + function TaylorIntegration._allocate_jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T}, + __x::AbstractArray{Taylor1{_S}, _N}, __dx::AbstractArray{Taylor1{_S}, _N}, + __params) where {_T<:Real, _S<:Number, _N} + order = __tT.order + nothing + end) ); # Constants for the initial declaration and initialization of arrays const _DECL_ARRAY = sanitize( Expr(:block, :(__var1 = Array{Taylor1{_S}}(undef, __var2)), - :(__var1 .= Taylor1( zero(_S), order )) - ) ) + :(__var1 .= Taylor1( zero(_S), order ))) +); """ `_make_parsed_jetcoeffs( ex, debug=false )` -This function constructs a new function, equivalent to the -differential equations, which exploits the mutating functions -of TaylorSeries.jl. +This function constructs the expressions of two new methods, the first equivalent to the +differential equations (jetcoeffs!), which exploits the mutating functions of TaylorSeries.jl, +and the second one (_allocate_jetcoeffs) preallocates any auxiliary `Taylor1` or +`Vector{Taylor1{T}}` needed. """ function _make_parsed_jetcoeffs(ex::Expr, debug=false) - # Extract the name, args and body of the function + # `fn` name of the function having the ODEs + # `fnargs` arguments of the function + # `fnbody` future transformed function body fn, fnargs, fnbody = _extract_parts(ex) debug && (println("****** _extract_parts ******"); @show(fn, fnargs, fnbody); println()) - # Set up new function - newfunction = _newhead(fn, fnargs) - - # for-loop block (recursion); it will contain the parsed body - forloopblock = Expr(:for, :(ord = 1:order-1), - Expr(:block, :(ordnext = ord + 1)) ) - - #= Transform graph representation of the body of the function - - `defspreamble` includes the definitions used for zeroth order (preamble) - - `fnbody` is the transformed function body, using mutating functions - from TaylorSeries, and is used within the for-loop block - - `retvar` is the (guessed) return variable, which defines the LHS - of the ODEs - =# - defspreamble, fnbody, retvar = _preamble_body(fnbody, fnargs, debug) - debug && (println("****** _preamble_body ******"); - @show(defspreamble); println(); @show(fnbody); println();) + # Set up the Expr for the new functions + new_jetcoeffs, new_allocjetcoeffs = _newhead(fn, fnargs) + + # Expr for the for-loop block for the recursion (of the `x` variable) + forloopblock = Expr(:for, :(ord = 1:order-1), Expr(:block, :(ordnext = ord + 1)) ) - #= Create body of recursion loop; temporary assignements may be needed. - - rec_preamb - - rec_fnbody - =# - debug && println("****** _recursionloop ******") - rec_preamb, rec_fnbody = _recursionloop(fnargs, retvar) + # Transform the graph representation of the body of the functions: + # defspreamble: inicializations used for the zeroth order (preamble) + # defsprealloc: definitions (declarations) of auxiliary Taylor1's + # fnbody: transformed function body, using mutating functions from TaylorSeries; + # used later within the recursion loop + # bkkeep: book-keeping structure having info of the variables + defspreamble, defsprealloc, fnbody, bkkeep = _preamble_body(fnbody, fnargs, debug) + debug && (println("****** _preamble_body ******"); + @show(defspreamble); println(); @show(defsprealloc); println(); + @show(fnbody); println(); @show(bkkeep); println()) - # Add preamble to newfunction - push!(newfunction.args[2].args, defspreamble..., rec_preamb) + # Create body of recursion loop; temporary assignements may be needed. + # rec_preamb: recursion loop for the preamble (first order correction) + # rec_fnbody: recursion loop for the body-function (recursion loop for higher orders) + rec_preamb, rec_fnbody = _recursionloop(fnargs, bkkeep) + debug && (println("****** _recursionloop ******"); + @show(rec_preamb); println(); @show(rec_fnbody); println()) - # Add parsed fnbody to forloopblock + # Add rec_fnbody to forloopblock push!(forloopblock.args[2].args, fnbody.args[1].args..., rec_fnbody) - # Push preamble and forloopblock to newfunction - push!(newfunction.args[2].args, forloopblock, Meta.parse("return nothing")) + # Add preamble and recursion body to `new_allocjetcoeffs` + push!(new_jetcoeffs.args[2].args, defspreamble..., rec_preamb) + + # Push preamble and forloopblock to `new_allocjetcoeffs` and return line + push!(new_jetcoeffs.args[2].args, forloopblock, Meta.parse("return nothing")) + + # Define the expressions of the returned vectors in `new_allocjetcoeffs` + push!(new_allocjetcoeffs.args[2].args, defsprealloc...) + if isempty(bkkeep.v_newvars) + ret_defsprealloc = :(Taylor1{_S}[]) + else + ret_defsprealloc = :([$(bkkeep.v_newvars...),]) + end + if isempty(bkkeep.v_arraydecl) + ret_tmparrays = :([Vector{Taylor1{_S}}(undef, 0),]) + else + ret_tmparrays = :([$(bkkeep.v_arraydecl...),]) + end + ret_ret = :(return $(ret_defsprealloc), $(ret_tmparrays)) + + # Add return line to `new_allocjetcoeffs` + push!(new_allocjetcoeffs.args[2].args, ret_ret) - # Rename variables of the body of the new function + # Rename variables in the calling form of the new methods if length(fnargs) == 3 - newfunction = subs(newfunction, Dict(:__tT => fnargs[3], + new_jetcoeffs = subs(new_jetcoeffs, Dict(:__tT => fnargs[3], :__params => fnargs[2], :(__x) => fnargs[1], :(__fn) => fn )) + new_allocjetcoeffs = subs(new_allocjetcoeffs, + Dict(:__tT => fnargs[3], :__params => fnargs[2], + :(__x) => fnargs[1], :(__fn) => fn )) + elseif length(fnargs) == 4 - newfunction = subs(newfunction, Dict(:__tT => fnargs[4], + new_jetcoeffs = subs(new_jetcoeffs, Dict(:__tT => fnargs[4], :__params => fnargs[3], :(__x) => fnargs[2], :(__dx) => fnargs[1], :(__fn) => fn )) + new_allocjetcoeffs = subs(new_allocjetcoeffs, + Dict(:__tT => fnargs[4], :__params => fnargs[3], + :(__x) => fnargs[2], :(__dx) => fnargs[1], :(__fn) => fn )) + + else + # A priori this is not needed + throw(ArgumentError("Wrong number of arguments in `fnargs`")) + end + + # Bring back old variable-names to `new_jecoeffs`; includes proper renaming + # of the temporal arrays + dd = Dict{Symbol, Expr}() + @inbounds for (ind, vnew) in enumerate(bkkeep.v_newvars) + merge!(dd, Dict(vnew => Expr(:ref, :__tmpTaylor, ind))) end + @inbounds for (ind, vnew) in enumerate(bkkeep.v_arraydecl) + merge!(dd, Dict(vnew => Expr(:ref, :__arrTaylor, ind))) + end + new_jetcoeffs = subs(new_jetcoeffs, dd) - return newfunction + return new_jetcoeffs, new_allocjetcoeffs end @@ -106,7 +203,6 @@ do not work). """ function _extract_parts(ex::Expr) - # Capture name, args and body # @capture( shortdef(ex), ffn_(ffnargs__) = fnbody_ ) || # throw(ArgumentError("It must be a function call:\n $ex")) @@ -154,27 +250,28 @@ end """ -`_capture_fn_args_body!(ex, vout::Dict{Symbol, Any})` +`_capture_fn_args_body!(ex, dd::Dict{Symbol, Any})` Captures the name of a function, arguments, body and other properties, returning them as the values of the dictionary `dd`, which is updated in place. """ -function _capture_fn_args_body!(ex::Expr, vout::Dict{Symbol, Any} = Dict()) +function _capture_fn_args_body!(ex::Expr, dd::Dict{Symbol, Any} = Dict()) exh = ExH(ex) if exh.head == :(=) || exh.head == :function - _capture_fn_args_body!.(ex.args, Ref(vout)) + _capture_fn_args_body!.(ex.args, Ref(dd)) elseif exh.head == :call - push!(vout, Pair(:fname, exh.args[1])) - push!(vout, Pair(:fnargs, exh.args[2:end])) + push!(dd, Pair(:fname, exh.args[1])) + push!(dd, Pair(:fnargs, exh.args[2:end])) elseif exh.head == :block - push!(vout, Pair(:fnbody, ex)) + push!(dd, Pair(:fnbody, ex)) elseif exh.head == :where - push!(vout, Pair(:where, ex.args[2:end])) - _capture_fn_args_body!(ex.args[1], vout) + push!(dd, Pair(:where, ex.args[2:end])) + _capture_fn_args_body!(ex.args[1], dd) end - vout + + return nothing end @@ -182,18 +279,21 @@ end """ `_newhead(fn, fnargs)` -Creates the head of the new method of `jetcoeffs!`. Here, -`fn` is the name of the passed function and `fnargs` is a -vector with its arguments (which are two or three). +Creates the head of the new method of `jetcoeffs!` and `_allocate_jetcoeffs`. +`fn` is the name of the passed function and `fnargs` is a vector with its arguments +defning the function (which are either three or four). """ function _newhead(fn, fnargs) - # Construct common elements of the new expression if length(fnargs) == 3 newfunction = copy(_HEAD_PARSEDFN_SCALAR) + newdeclfunc = copy(_HEAD_ALLOC_TAYLOR1_SCALAR) + elseif length(fnargs) == 4 newfunction = copy(_HEAD_PARSEDFN_VECTOR) + newdeclfunc = copy(_HEAD_ALLOC_TAYLOR1_VECTOR) + else throw(ArgumentError( "Wrong number of arguments in the definition of the function $fn")) @@ -201,8 +301,9 @@ function _newhead(fn, fnargs) # Delete the irrelevant `nothing` pop!(newfunction.args[2].args) + pop!(newdeclfunc.args[2].args) - return newfunction + return newfunction, newdeclfunc end @@ -210,72 +311,61 @@ end """ `_preamble_body(fnbody, fnargs, debug=false)` -Returns the preamble, the body and a guessed return -variable, which will be used to build the parsed -function. `fnbody` is the expression with the body of -the original function, `fnargs` is a vector of symbols +Returns expressions for the preamble, the declaration of +arrays, the body and the bookkeeping struct, which will be used to build +the new functions. `fnbody` is the expression with the body of +the original function (already adapted), `fnargs` is a vector of symbols of the original diferential equations function. """ function _preamble_body(fnbody, fnargs, debug=false) + # Inicialize BookKeeping struct + bkkeep = BookKeeping() - # Bookkeeping: - # v_vars: List of symbols with created variables - # v_assign: Numeric assignments (to be deleted) - v_vars = Union{Symbol,Expr}[] - v_assign = Dict{Union{Symbol,Expr}, Number}() - - #= Rename vars to have the body in non-indexed form - `d_indx` is a dictionary mapping new variables (symbols) to old - (perhaps indexed) symbols - =# - fnbody, d_indx = _rename_indexedvars(fnbody) + # Rename vars to have the body in non-indexed form; bkkeep has different entries + # for bookkeeping variables/symbolds, including indexed ones + fnbody, bkkeep.d_indx = _rename_indexedvars(fnbody) debug && (println("------ _rename_indexedvars ------"); - @show(fnbody); println(); @show(d_indx); println()) - - #= Create newfnbody - - v_newindx: symbols of auxiliary indexed vars - - v_arraydecl: symbols which are explicitly declared as Array or Vector - - `newfnbody` corresponds to `fnbody`, cleaned (without irrelevant comments) - and with all new variables in place - =# - newfnbody, v_newindx, v_arraydecl = _newfnbody(fnbody, fnargs, d_indx) + @show(fnbody); println(); @show(bkkeep.d_indx); println()) + + # Create `newfnbody` which corresponds to `fnbody`, cleaned (without irrelevant comments) + # and with all new variables in place; bkkeep.d_indx is updated + newfnbody = _newfnbody!(fnbody, fnargs, bkkeep) debug && (println("------ _newfnbody ------"); - @show(v_newindx); println(); @show(v_arraydecl); println(); - @show(newfnbody); println()) + @show(newfnbody); println(); @show(bkkeep); println()) - # Parse `newfnbody` and create `prepreamble`, updating the bookkeeping vectors. - # The returned `newfnbody` and `prepreamble` use the mutating functions - # of TaylorSeries. - prepreamble = Expr(:block,) - _parse_newfnbody!(newfnbody, prepreamble, v_vars, v_assign, d_indx, - v_newindx, v_arraydecl) - preamble = prepreamble.args[1] + # Parse `newfnbody` and create `prepreamble` and `prealloc`, updating `bkkeep`. + # These objects use the mutating functions from TaylorSeries. + preamble = Expr(:block,) + prealloc = Expr(:block,) + _parse_newfnbody!(newfnbody, preamble, prealloc, bkkeep, false) + # Get rid of the initial `:block` + preamble = preamble.args[1] + prealloc = prealloc.args[1] # Substitute numeric assignments in new function's body - newfnbody = subs(newfnbody, Dict(v_assign)) - preamble = subs(preamble, Dict(v_assign)) + newfnbody = subs(newfnbody, bkkeep.d_assign) + preamble = subs(preamble, bkkeep.d_assign) + prealloc = subs(prealloc, bkkeep.d_assign) debug && (println("------ _parse_newfnbody! ------"); @show(newfnbody); println(); @show(preamble); println(); - @show(v_vars); println(); @show(d_indx); println(); - @show(v_assign); println(); @show(v_newindx); println()) + @show(prealloc); println(); @show(bkkeep); println()) # Include the assignement of indexed auxiliary variables - defspreamble = _defs_preamble!(preamble, fnargs, - d_indx, v_newindx, v_arraydecl, Union{Symbol,Expr}[], empty(d_indx)) - # preamble = subs(preamble, d_indx) - + defsprealloc = _defs_allocs!(prealloc, fnargs, bkkeep, false) + preamble = subs(preamble, bkkeep.d_indx) + defspreamble = Expr[preamble.args...] # Bring back substitutions - newfnbody = subs(newfnbody, d_indx) + newfnbody = subs(newfnbody, bkkeep.d_indx) - # Define retvar; for scalar eqs is the last entry included in v_vars - retvar = length(fnargs) == 3 ? subs(v_vars[end], d_indx) : fnargs[1] + # Define retvar; for scalar eqs is the last entry included in v_newvars + bkkeep.retvar = length(fnargs) == 3 ? subs(bkkeep.v_newvars[end], bkkeep.d_indx) : fnargs[1] - debug && (println("------ _defs_preamble! ------"); - @show(d_indx); println(); @show(v_newindx); println(); - @show(newfnbody); println(); @show(retvar); println()) + debug && (println("------ _defs_allocs! ------"); + @show(defsprealloc); println(); @show(preamble); println(); + @show(newfnbody); println(); @show(bkkeep); println()) - return defspreamble, newfnbody, retvar + return defspreamble, defsprealloc, newfnbody, bkkeep end @@ -295,90 +385,104 @@ function _rename_indexedvars(fnbody) indexed_vars = findex(:(_X[_i...]), fnbody) st = Dict(ivar => genname() for ivar in indexed_vars) new_fnbody = subs(fnbody, st) - return new_fnbody, Dict(v => k for (k, v) in st) + return new_fnbody, Dict{Symbol, Expr}(v => k for (k, v) in st) end """ -`_newfnbody(fnbody, fnargs, d_indx)` +`_newfnbody!(fnbody, fnargs, bkkeep)` Returns a new (modified) body of the function, a priori unfolding -the expression graph (AST) as unary and binary calls, and a vector -(`v_newindx`) with the name of auxiliary indexed variables. +the expression graph (AST) as unary and binary calls, and updates +the bookkeeping structure bkkeep. """ -function _newfnbody(fnbody, fnargs, d_indx) - # `fnbody` is assumed to be a `:block` `Expr` +function _newfnbody!(fnbody, fnargs, bkkeep::BookKeeping) + # `fnbody` is assumed to be a `:block`-Expr newfnbody = Expr(:block,) - v_newindx = Symbol[] - v_arraydecl = Symbol[] - # The magic happens HERE!! + # Local definition for possible args for `@threads` call + local v_threads_args = (Expr(:., :Threads, QuoteNode(Symbol("@threads"))), + Symbol("@threads")) + + # Magic happens HERE!! # Each line of fnbody (fnbody.args) is parsed separately for (i, ex) in enumerate(fnbody.args) if isa(ex, Expr) ex_head = ex.head - # Ignore the following cases - (ex_head == :return) && continue + # The `return` statement is treated separately + if ex.head == :return + @assert length(ex.args) == 1 # one sole returned value + + # Do nothing if it is `nothing` or if its in fnargs + (ex.args[1] == :nothing || ex.args[1] in fnargs) && continue + + # Otherwise, create a new assignment, and process it + new_aux_var = genname() + ex = :($new_aux_var = $(ex.args[1])) + ex_head = ex.head + end # Treat `for` loops, `Threads.@threads for` and `if` blocks separately if ex_head == :block - newblock, tmp_newindx, tmp_arraydecl = _newfnbody(ex, fnargs, d_indx) + newblock = _newfnbody!(ex, fnargs, bkkeep) push!(newfnbody.args, newblock ) - append!(v_newindx, tmp_newindx) - append!(v_arraydecl, tmp_arraydecl) + elseif ex_head == :for push!(newfnbody.args, Expr(:for, ex.args[1])) - loopbody, tmp_newindx, tmp_arraydecl = _newfnbody( ex.args[2], fnargs, d_indx ) + loopbody = _newfnbody!( ex.args[2], fnargs, bkkeep ) push!(newfnbody.args[end].args, loopbody) - append!(v_newindx, tmp_newindx) - append!(v_arraydecl, tmp_arraydecl) + elseif ex_head == :macrocall + # TODO: Deal with `@inbounds` # Deal with `Threads.@threads` and `@threads` cases - if ex.args[1] in [Expr(:., :Threads, QuoteNode(Symbol("@threads"))), Symbol("@threads")] + if ex.args[1] in v_threads_args push!(newfnbody.args, Expr(:macrocall, ex.args[1])) # Although here `ex.args[2]` is a `LineNumberNode`, - # we add it to `newfnbody` because `@threads` call expressions require 3 args + # we add it to `newfnbody` because `@threads` call expressions + # require 3 args push!(newfnbody.args[end].args, ex.args[2]) # Since `@threads` is called before a `for` loop, we deal # with `ex.args[3]` as a `for` loop push!(newfnbody.args[end].args, Expr(:for, ex.args[3].args[1])) - atthreadsbody, tmp_newindx, tmp_arraydecl = _newfnbody( ex.args[3].args[2], fnargs, d_indx ) + + atthreadsbody = _newfnbody!( ex.args[3].args[2], fnargs, bkkeep ) push!(newfnbody.args[end].args[end].args, atthreadsbody) - append!(v_newindx, tmp_newindx) - append!(v_arraydecl, tmp_arraydecl) + else - # If macro not implemented, throw an `ArgumentError` + # If macro is not implemented, throw an `ArgumentError` throw(ArgumentError("Macro $(ex.args[1]) is not yet implemented")) end + elseif ex_head == :if # The first argument of an `if` expression is the condition, the # second one is the block the condition is true, and the third # is the `else`-block or an `ifelse`-block. push!(newfnbody.args, Expr(:if, ex.args[1])) + for exx in ex.args[2:end] if exx.head == :elseif exxx = Expr(:block, Expr(:if, exx.args[1].args[2], exx.args[2:end]...)) - else #if exx.head == :block # `then` or `else` blocks + else exxx = exx end - ifbody, tmp_newindx, tmp_arraydecl = _newfnbody(exxx, fnargs, d_indx) + + ifbody = _newfnbody!(exxx, fnargs, bkkeep) push!(newfnbody.args[end].args, ifbody) - append!(v_newindx, tmp_newindx) - append!(v_arraydecl, tmp_arraydecl) end - elseif ex_head == :(=) || ex_head == :call # assignements or function calls + elseif ex_head == :(=) || ex_head == :call + # Assignements or function calls ex_lhs = ex.args[1] ex_rhs = ex.args[2] + # TODO: Include cases where the definition uses `similar` # Case of explicit declaration of Array or Vector - # TO-DO: Include cases where the definition uses `similar` if !isempty(findex(:(_Array{_TT...}), ex)) || !isempty(findex(:(Vector{_TT...}), ex)) - push!(v_arraydecl, ex_lhs) + inbookkeeping(ex_lhs, bkkeep) || push!(bkkeep.v_arraydecl, ex_lhs) push!(newfnbody.args, ex) continue end @@ -390,24 +494,26 @@ function _newfnbody(fnbody, fnargs, d_indx) # Construct new expressions if ex_rhs == fnargs[3] # params for i in eachindex(ex_lhs.args) - ex_indx = Expr(:local, Expr(:(=), ex_lhs.args[i], Expr(:ref, ex_rhs, i))) + ex_indx = Expr(:local, Expr(:(=), ex_lhs.args[i], + Expr(:ref, ex_rhs, i))) push!(newfnbody.args, ex_indx) end + elseif ex_rhs == fnargs[2] # `x` for i in eachindex(ex_lhs.args) ex_indx = Expr(:(=), ex_lhs.args[i], Expr(:ref, ex_rhs, i)) push!(newfnbody.args, ex_indx) end + else throw(ArgumentError("Assignment no implemented in $(ex)")) end - # end continue end # Unfold AST graph - nex = deepcopy(ex) + # nex = deepcopy(ex) # try # nex = to_expr(ExGraph(simplify(ex))) # catch @@ -419,34 +525,35 @@ function _newfnbody(fnbody, fnargs, d_indx) push!(newfnbody.args, nex.args[2:end]...) # Bookkeeping of indexed vars, to define assignements - isindx_lhs = haskey(d_indx, ex_lhs) + isindx_lhs = haskey(bkkeep.d_indx, ex_lhs) for nexargs in nex.args[2:end] - # (nexargs.head == :line || - # haskey(d_indx, nexargs.args[1])) && continue - haskey(d_indx, nexargs.args[1]) && continue + haskey(bkkeep.d_indx, nexargs.args[1]) && continue vars_nex = find_vars(nexargs) - # any(haskey.(d_indx, vars_nex[:])) && - # !in(vars_nex[1], v_newindx) && - (isindx_lhs || vars_nex[1] != ex_lhs) && push!(v_newindx, vars_nex[1]) + # any(haskey.(bkkeep.d_indx, Ref(vars_nex))) && + # !in(vars_nex[1], bkkeep.v_newvars) && + (isindx_lhs || vars_nex[1] != ex_lhs) && + !inbookkeeping(vars_nex[1], bkkeep) && + push!(bkkeep.v_newvars, vars_nex[1]) end + elseif (ex_head == :local) || (ex_head == :continue) || (ex_head == :break) - # If declared as `local` or `continue`, copy `ex` as it is. - # In some cases this, using `local` helps performance. Very + # If declared as `local`, or it is a `continue` or `break`, copy `ex` as it is. + # In some cases using `local` helps performance. Very # useful for including (numeric) constants push!(newfnbody.args, ex) - # + else # If not implemented, throw an `ArgumentError` throw(ArgumentError("$(ex_head) is not yet implemented; $(typeof(ex))")) end - # + elseif isa(ex, LineNumberNode) continue # Ignore `LineNumberNode`s + else # In any other case, throw an `ArgumentError` throw(ArgumentError("$ex is not an `Expr`; $(typeof(ex))")) - # end end @@ -455,61 +562,72 @@ function _newfnbody(fnbody, fnargs, d_indx) newfnbody = Expr(:block, newfnbody) end - return newfnbody, v_newindx, v_arraydecl + return newfnbody end """ -`_parse_newfnbody!(ex, preex, v_vars, v_assign, d_indx, v_newindx, v_arraydecl, - [inloop=false])` +_parse_newfnbody!(ex::Expr, preex::Expr, prealloc::Expr, bkkeep::BookKeeping, inloop::Bool) Parses `ex` (the new body of the function) replacing the expressions to use the mutating functions of TaylorSeries, and building the preamble -`preex`. This is done by traversing recursively the args of `ex`, updating -the bookkeeping vectors `v_vars` and `v_assign`. `d_indx` is -used to substitute back the explicit indexed variables. +`preex` and `prealloc` expressions. This is done by traversing recursively (again) +the args of `ex`, updating the bookkeeping struct `bkkeep`, in particular +the fieldnames `v_newvars` and `d_assign`. """ -function _parse_newfnbody!(ex::Expr, preex::Expr, - v_vars, v_assign, d_indx, v_newindx, v_arraydecl, inloop=false) - +function _parse_newfnbody!(ex::Expr, preex::Expr, prealloc::Expr, bkkeep::BookKeeping, + inloop::Bool) # Numeric assignements to be deleted indx_rm = Int[] + # Definition for possible args for `@threads` call + local v_threads_args = (Expr(:., :Threads, QuoteNode(Symbol("@threads"))), Symbol("@threads")) + # Magic happens HERE - # Each line of ex (ex.args) is parsed separately + # Each line of `ex` (i.e., `ex.args`) is parsed separately for (i, aa) in enumerate(ex.args) - - # Treat for loops, @threads for loops, blocks and if blocks separately + # Treat for-loops, @threads macro, blocks and if-blocks separately + # (through `inloop=true`) if (aa.head == :for) push!(preex.args, Expr(:for, aa.args[1])) - _parse_newfnbody!(aa, preex.args[end], - v_vars, v_assign, d_indx, v_newindx, v_arraydecl, true) - # - elseif (aa.head == :macrocall && aa.args[1] in [Expr(:., :Threads, QuoteNode(Symbol("@threads"))), Symbol("@threads")]) + push!(prealloc.args, Expr(:for, aa.args[1])) + + _parse_newfnbody!(aa, preex.args[end], prealloc.args[end], bkkeep, true) + + elseif (aa.head == :macrocall && aa.args[1] in v_threads_args) push!(preex.args, Expr(:macrocall, aa.args[1])) push!(preex.args[end].args, aa.args[2]) push!(preex.args[end].args, Expr(:for, aa.args[3].args[1])) - _parse_newfnbody!(aa.args[3], preex.args[end].args[end], - v_vars, v_assign, d_indx, v_newindx, v_arraydecl, true) - # - elseif (aa.head == :block) + push!(prealloc.args, Expr(:macrocall, aa.args[1])) + push!(prealloc.args[end].args, aa.args[2]) + push!(prealloc.args[end].args, Expr(:for, aa.args[3].args[1])) + _parse_newfnbody!(aa.args[3], preex.args[end].args[end], prealloc.args[end].args[end], + bkkeep, true) + + elseif (aa.head == :block) push!(preex.args, Expr(:block)) - _parse_newfnbody!(aa, preex.args[end], - v_vars, v_assign, d_indx, v_newindx, v_arraydecl) - # + push!(prealloc.args, Expr(:block)) + + # `inloop` needs to be `false` here + _parse_newfnbody!(aa, preex.args[end], prealloc.args[end], bkkeep, false) + elseif (aa.head == :if) push!(preex.args, Expr(:if, aa.args[1])) + push!(prealloc.args, Expr(:if, aa.args[1])) + for exx in aa.args[2:end] push!(preex.args[end].args, Expr(:block)) + push!(prealloc.args[end].args, Expr(:block)) _parse_newfnbody!(exx, preex.args[end].args[end], - v_vars, v_assign, d_indx, v_newindx, v_arraydecl) + prealloc.args[end].args[end], bkkeep, inloop) end - # - elseif (aa.head == :(=)) + elseif (aa.head == :(=)) + # Skip parsing the first for-loop expression (`inloop=true`), which is of + # the form `:(i = ...)` i == 1 && inloop && continue # Main case @@ -519,39 +637,43 @@ function _parse_newfnbody!(ex::Expr, preex::Expr, # Replace expressions when needed, and bookkeeping if isa(aa_rhs, Expr) (aa_rhs.args[1] == :eachindex || aa_rhs.head == :(:)) && continue + # Replace expression - _replace_expr!(ex, preex, i, aa_lhs, aa_rhs, v_vars, d_indx, v_newindx) + _replace_expr!(ex, preex, prealloc, i, aa_lhs, aa_rhs, bkkeep) # Remove new aa_lhs declaration from `ex`; it is already declared - in(aa_lhs, v_arraydecl) && push!(indx_rm, i) - # - elseif isa(aa_rhs, Symbol) # occurs when there is a simple assignment + in(aa_lhs, bkkeep.v_arraydecl) && push!(indx_rm, i) + elseif isa(aa_rhs, Symbol) # occurs when there is a simple assignment bb = subs(aa, Dict(aa_rhs => :(identity($aa_rhs)))) bb_lhs = bb.args[1] bb_rhs = bb.args[2] + # Replace expression - _replace_expr!(ex, preex, i, bb_lhs, bb_rhs, v_vars, d_indx, v_newindx) - # + _replace_expr!(ex, preex, prealloc, i, bb_lhs, bb_rhs, bkkeep) + elseif isa(aa_rhs, Number) # case of numeric values - push!(v_assign, aa_lhs => aa_rhs) + push!(bkkeep.d_assign, aa_lhs => aa_rhs) push!(indx_rm, i) - # + if aa_lhs ∈ bkkeep.v_newvars + iis = findall(x->x==aa_lhs, bkkeep.v_newvars) + deleteat!(bkkeep.v_newvars, iis) + end + else # needed? - error("Either $aa or $typeof(aa_rhs[2]) are different from `Expr`, `Symbol` or `Number`") - # + error("Either $aa or $typeof(aa_rhs[2]) are not an `Expr`, `Symbol` or `Number`") end - # + elseif aa.head == :local - # If declared as `local`, copy `ex` as it is, and delete it from - # the recursion body + # If declared as `local`, copy `ex` as it is, and delete it from the recursion body push!(preex.args, aa) + push!(prealloc.args, aa) push!(indx_rm, i) # delete associated expr in body function - # - else + else # Unrecognized head; pass the expression as it is push!(preex.args, aa) + push!(prealloc.args, aa) end end @@ -564,285 +686,370 @@ end """ -`_replace_expr!(ex, preex, i, aalhs, aarhs, v_vars, d_indx, v_newindx)` +`_replace_expr!(ex::Expr, preex::Expr, , prealloc::Expr, i::Int, aalhs, aarhs, + bkkeep::BookKeeping)` -Replaces the calls in `ex.args[i]`, and updates `preex` with the definitions -of the expressions, based on the the LHS (`aalhs`) and RHS (`aarhs`) of the -base assignment. The bookkeeping vectors (`v_vars`, `v_newindx`) -are updated. `d_indx` is used to bring back the indexed variables. +Replaces the calls in `ex.args[i]`, and updates `preex` and `prealloc` with the +appropriate expressions, based on the the LHS (`aalhs`) and RHS (`aarhs`) of the +base assignment. The bookkeeping struct is updated (`v_newvars`) within `_replacecalls!`. +`d_indx` is used to bring back the indexed variables. """ -function _replace_expr!(ex::Expr, preex::Expr, i::Int, - aalhs, aarhs, v_vars, d_indx, v_newindx) - +function _replace_expr!(ex::Expr, preex::Expr, prealloc::Expr, i::Int, + aalhs, aarhs, bkkeep::BookKeeping) # Replace calls - fnexpr, def_fnexpr, auxfnexpr = _replacecalls!(aarhs, aalhs, v_vars) + fnexpr, def_fnexpr, auxfnexpr, def_alloc, aux_alloc = _replacecalls!(bkkeep, aarhs, aalhs) # Bring back indexed variables - fnexpr = subs(fnexpr, d_indx) - def_fnexpr = subs(def_fnexpr, d_indx) - auxfnexpr = subs(auxfnexpr, d_indx) + fnexpr = subs(fnexpr, bkkeep.d_indx) + def_fnexpr = subs(def_fnexpr, bkkeep.d_indx) + auxfnexpr = subs(auxfnexpr, bkkeep.d_indx) + def_alloc = subs(def_alloc, bkkeep.d_indx) + aux_alloc = subs(aux_alloc, bkkeep.d_indx) + + # Update `ex`, `preex` and `prealloc` + if (fnexpr != def_fnexpr) + ex.args[i] = fnexpr + end - # Update `ex` and `preex` - push!(preex.args, def_fnexpr) - ex.args[i] = fnexpr + if def_fnexpr.head != :nothing + if def_fnexpr.head == :block + append!(preex.args, def_fnexpr.args) + else + push!(preex.args, def_fnexpr) + end + end + + if def_alloc.head != :nothing + if def_alloc.head == :block + append!(prealloc.args, def_alloc.args) + else + push!(prealloc.args, def_alloc) + end + end # Same for the auxiliary expressions if auxfnexpr.head != :nothing - push!(preex.args, auxfnexpr) - # in(aalhs, v_newindx) && push!(v_newindx, auxfnexpr.args[1]) - push!(v_newindx, auxfnexpr.args[1]) + if auxfnexpr.head == :block + append!(preex.args, auxfnexpr.args) + else + push!(preex.args, auxfnexpr) + end + aux_alloc.head == :nothing || push!(prealloc.args, aux_alloc) end - !in(aalhs, v_vars) && push!(v_vars, subs(aalhs, d_indx)) return nothing end """ -`_replacecalls!(fnold, newvar, v_vars)` +`_replacecalls!(bkkeep, fnold, newvar)` Replaces the symbols of unary and binary calls of the expression `fnold`, which defines `newvar`, by the mutating functions in TaylorSeries.jl. -The vector `v_vars` is updated if new auxiliary -variables are introduced. +The vector `bkkeep.v_vars` is updated if new auxiliary +variables are introduced (bookkeeping). """ -function _replacecalls!(fnold::Expr, newvar::Symbol, v_vars) - +function _replacecalls!(bkkeep::BookKeeping, fnold::Expr, newvar::Symbol) ll = length(fnold.args) dcall = fnold.args[1] newarg1 = fnold.args[2] - # If call is not in mutating functions dictionaries, copy original one + # If call is not in the mutating functions dictionaries, copy original one + # to def_fnexpr and def_alloc, except if it is an Array/Vector definition if !( in(dcall, keys(TaylorSeries._dict_unary_calls)) || - in(dcall, keys(TaylorSeries._dict_binary_calls)) ) + in(dcall, keys(TaylorSeries._dict_binary_calls)) ) fnexpr = :($newvar = $fnold) - def_fnexpr = fnexpr - aux_fnexpr = Expr(:nothing) - return fnexpr, def_fnexpr, aux_fnexpr + if isa(dcall, Expr) && (dcall.args[1] == :Array || dcall.args[1] == :Vector) + def_fnexpr = Expr(:nothing) + inbookkeeping(newvar, bkkeep) || push!(bkkeep.v_arraydecl, newvar) + else + def_fnexpr = fnexpr + inbookkeeping(newvar, bkkeep) || push!(bkkeep.v_newvars, newvar) + end + def_alloc = fnexpr + + return fnexpr, def_fnexpr, Expr(:nothing), def_alloc, Expr(:nothing) end - if ll == 2 + # Bookkeeping + inbookkeeping(newvar, bkkeep) || push!(bkkeep.v_newvars, newvar) + + # Initializing + def_alloc = Expr(:nothing) + aux_alloc = Expr(:nothing) + if ll == 2 # Unary call # Replacements fnexpr, def_fnexpr, aux_fnexpr = TaylorSeries._dict_unary_calls[dcall] fnexpr = subs(fnexpr, Dict(:_res => newvar, :_arg1 => newarg1, :_k => :ord)) - def_fnexpr = :( _res = Taylor1($(def_fnexpr.args[2]), order) ) + + def_alloc = :( _res = Taylor1($(def_fnexpr.args[2]), order) ) + def_alloc = subs(def_alloc, + Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), :_k => :ord)) + + def_fnexpr = Expr(:block, + :(_res.coeffs[1] = $(def_fnexpr.args[2])), + :(_res.coeffs[2:order+1] .= zero(_res.coeffs[1])) ) def_fnexpr = subs(def_fnexpr, - Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), - :_k => :ord)) + Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), :_k => :ord)) + # def_fnexpr = Expr(:block, + # :(_res[0] = $(def_fnexpr.args[2])), + # :(_res[1:order] .= zero(_res[0])) ) + # def_fnexpr = subs(def_fnexpr, + # Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), :_k => :ord)) # Auxiliary expression if aux_fnexpr.head != :nothing newaux = genname() - aux_fnexpr = :( _res = Taylor1($(aux_fnexpr.args[2]), order) ) + + aux_alloc = :( _res = Taylor1($(aux_fnexpr.args[2]), order) ) + aux_alloc = subs(aux_alloc, + Dict(:_res => newaux, :_arg1 => :(constant_term($(newarg1))), :_aux => newaux)) + + aux_fnexpr = Expr(:block, + :(_res.coeffs[1] = $(aux_fnexpr.args[2])), + :(_res.coeffs[2:order+1] .= zero(_res.coeffs[1])) ) aux_fnexpr = subs(aux_fnexpr, - Dict(:_res => newaux, :_arg1 => :(constant_term($(newarg1))), - :_aux => newaux)) + Dict(:_res => newaux, :_arg1 => :(constant_term($(newarg1))), :_aux => newaux)) + fnexpr = subs(fnexpr, Dict(:_aux => newaux)) - !in(newaux, v_vars) && push!(v_vars, newaux) + if newvar ∈ bkkeep.v_arraydecl + push!(bkkeep.v_arraydecl, newaux) + else + push!(bkkeep.v_newvars, newaux) + end end - elseif ll == 3 + elseif ll == 3 # Binary call; no auxiliary expressions needed newarg2 = fnold.args[3] # Replacements fnexpr, def_fnexpr, aux_fnexpr = TaylorSeries._dict_binary_calls[dcall] fnexpr = subs(fnexpr, - Dict(:_res => newvar, :_arg1 => newarg1, :_arg2 => newarg2, - :_k => :ord)) + Dict(:_res => newvar, :_arg1 => newarg1, :_arg2 => newarg2, :_k => :ord)) - def_fnexpr = :(_res = Taylor1($(def_fnexpr.args[2]), order) ) + def_alloc = :(_res = Taylor1($(def_fnexpr.args[2]), order) ) + def_alloc = subs(def_alloc, Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), + :_arg2 => :(constant_term($(newarg2))), :_k => :ord) ) + + def_fnexpr = Expr(:block, + :(_res.coeffs[1] = $(def_fnexpr.args[2])), + :(_res.coeffs[2:order+1] .= zero(_res.coeffs[1])) ) def_fnexpr = subs(def_fnexpr, - Dict(:_res => newvar, - :_arg1 => :(constant_term($(newarg1))), - :_arg2 => :(constant_term($(newarg2))), - :_k => :ord)) + Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), + :_arg2 => :(constant_term($(newarg2))), :_k => :ord)) + # def_fnexpr = Expr(:block, + # :(_res[0] = $(def_fnexpr.args[2])), + # :(_res[1:order] .= zero(_res[0])) ) + # def_fnexpr = subs(def_fnexpr, + # Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), + # :_arg2 => :(constant_term($(newarg2))), :_k => :ord)) + else - # Recognized call, but not as a unary or binary call + # Recognized call, but not a unary or binary call; copy expression fnexpr = :($newvar = $fnold) def_fnexpr = fnexpr aux_fnexpr = Expr(:nothing) end - return fnexpr, def_fnexpr, aux_fnexpr + return fnexpr, def_fnexpr, aux_fnexpr, def_alloc, aux_alloc end """ -`_defs_preamble!(preamble, fnargs, d_indx, v_newindx, v_arraydecl, v_preamb, - d_decl, [inloop=false])` +`_defs_allocs!(preamble, fnargs, bkkeep, [inloop=false, ex_aux::Expr(:block,)])` Returns a vector with expressions defining the auxiliary variables -in the preamble; it may modify `d_indx` if new variables are introduced. -`v_preamb` is for bookkeeping the introduced variables. +in the preamble, and the declaration of the arrays. This function +may modify `bkkeep.d_indx` if new variables are introduced. +`bkkeep.v_preamb` is for bookkeeping the introduced variables. """ -function _defs_preamble!(preamble::Expr, fnargs, - d_indx, v_newindx, v_arraydecl, v_preamb, d_decl, inloop::Bool=false, - ex_aux::Expr = Expr(:block,)) +function _defs_allocs!(preamble::Expr, fnargs, bkkeep::BookKeeping, + inloop::Bool, ex_aux::Expr = Expr(:block,)) # Initializations defspreamble = Expr[] + defs_alloc = Expr[] # declaration of arrays - for (i, ex) in enumerate(preamble.args) + # Local definition for possible args for `@threads` call + local v_threads_args = (Expr(:., :Threads, QuoteNode(Symbol("@threads"))), Symbol("@threads")) - if isa(ex, Expr) + for ex in preamble.args + isa(ex, Expr) || throw(ArgumentError("$ex is not an `Expr`; $(typeof(ex))")) + + # Treat block, for loops, @threads for loops, if separately + if (ex.head == :block) + newdefspr = _defs_allocs!(ex, fnargs, bkkeep, inloop, ex_aux) + + append!(defspreamble, newdefspr) + + elseif (ex.head == :for) + push!(ex_aux.args, ex.args[1]) + + newdefspr = _defs_allocs!(ex.args[2], fnargs, bkkeep, true, ex_aux) + + append!(defspreamble, newdefspr) + pop!(ex_aux.args) + + elseif (ex.head == :macrocall && ex.args[1] in v_threads_args) + push!(ex_aux.args, ex.args[3].args[1]) + + newdefspr = _defs_allocs!(ex.args[3].args[2], fnargs, bkkeep, true, ex_aux) + + append!(defspreamble, newdefspr) + pop!(ex_aux.args) + + elseif (ex.head == :if) + for exx in ex.args[2:end] + newdefspr = _defs_allocs!(exx, fnargs, bkkeep, inloop, ex_aux) - # Treat block, for loops, @threads for loops, if separately - if (ex.head == :block) - newdefspr = _defs_preamble!(ex, fnargs, d_indx, - v_newindx, v_arraydecl, v_preamb, d_decl, inloop, ex_aux) - append!(defspreamble, newdefspr) - elseif (ex.head == :for) - push!(ex_aux.args, ex.args[1]) - newdefspr = _defs_preamble!(ex.args[2], fnargs, d_indx, - v_newindx, v_arraydecl, v_preamb, d_decl, true, ex_aux) - append!(defspreamble, newdefspr) - pop!(ex_aux.args) - elseif (ex.head == :macrocall && ex.args[1] in [Expr(:., :Threads, QuoteNode(Symbol("@threads"))), Symbol("@threads")]) - push!(ex_aux.args, ex.args[3].args[1]) - newdefspr = _defs_preamble!(ex.args[3].args[2], fnargs, d_indx, - v_newindx, v_arraydecl, v_preamb, d_decl, true, ex_aux) append!(defspreamble, newdefspr) - pop!(ex_aux.args) - elseif (ex.head == :if) - for exx in ex.args[2:end] - newdefspr = _defs_preamble!(exx, fnargs, d_indx, - v_newindx, v_arraydecl, v_preamb, d_decl, inloop, ex_aux) - append!(defspreamble, newdefspr) + end + continue + + elseif (ex.head == :local) + push!(defspreamble, subs(ex, bkkeep.d_indx)) + continue + + elseif (ex.head == :(=)) + # `ex.head` is a :(=) of some kind + alhs = ex.args[1] + arhs = ex.args[2] + arhs = subs(ex.args[2], bkkeep.d_indx) # substitute updated vars in rhs + + # Outside of a loop + if !inloop + # `alhs` is declared as an array + if in(alhs, bkkeep.v_arraydecl) || !in(alhs, bkkeep.v_preamb) + push!(defspreamble, subs(ex, bkkeep.d_decl)) + push!(bkkeep.v_preamb, alhs) end - continue - elseif (ex.head == :(=)) - # `ex.head` is a :(=) of some kind - alhs = ex.args[1] - arhs = subs(ex.args[2], d_indx) # substitute updated vars in rhs + continue + end - # Outside of a loop - if !inloop - in(alhs, v_preamb) && continue - push!(defspreamble, subs(ex, d_decl)) - push!(v_preamb, alhs) - continue + # Inside a loop + if isindexed(alhs) + # `var1` may be a vector or a matrix; declaring it is subtle + var1 = alhs.args[1] + while isa(var1, Expr) + var1 = var1.args[1] end + (in(var1, fnargs) || in(var1, bkkeep.v_arraydecl) || + in(alhs, bkkeep.v_preamb)) && continue + + # indices of var1 + d_subs = Dict(veexx.args[1] => veexx.args[2] for veexx in ex_aux.args) + indx1 = alhs.args[2:end] + ex_tuple = :( [$(indx1...)] ) + ex_tuple = subs(ex_tuple, d_subs) + ex_tuple = :( size( $(Expr(:tuple, ex_tuple.args...)) )) + exx = subs(copy(_DECL_ARRAY), + Dict( :__var1 => :($var1), :__var2 => :($ex_tuple)) ) + + # Bookkeeping + push!(defspreamble, exx.args...) + push!(bkkeep.v_preamb, var1) + continue - # Inside a loop - if isindexed(alhs) - - # `var1` may be a vector or a matrix; declaring it is subtle - var1 = alhs.args[1] - (in(var1, fnargs) || in(var1, v_arraydecl) || - in(alhs, v_preamb)) && continue - - # indices of var1 - d_subs = Dict(veexx.args[1] => veexx.args[2] - for veexx in ex_aux.args) - indx1 = alhs.args[2:end] - ex_tuple = :( [$(indx1...)] ) - ex_tuple = subs(ex_tuple, d_subs) - ex_tuple = :( size( $(Expr(:tuple, ex_tuple.args...)) )) - exx = subs(_DECL_ARRAY, Dict(:__var1 => :($var1), - :__var2 => :($ex_tuple)) ) - push!(defspreamble, exx.args...) - push!(v_preamb, var1) - continue - # - elseif in(alhs, v_newindx) || isindexed(arhs) - # `alhs` is an aux indexed var, so something in `arhs` - # is indexed. - - in(alhs, v_preamb) && continue - - vars_indexed = findex(:(_X[_i...]), arhs) - - # NOTE: Uses the size of the var with more indices - # to define the declaration of the new array. - iimax, ii_indx = findmax( - [length(find_indices(aa)[1]) for aa in vars_indexed] ) - var1 = vars_indexed[ii_indx].args[1] - indx1 = vars_indexed[ii_indx].args[2:end] - - exx_indx = ones(Int, length(indx1)) - push!(d_indx, alhs => :($alhs[$(indx1...)]) ) - push!(d_decl, alhs => :($alhs[$exx_indx...])) - exx = subs(_DECL_ARRAY, - Dict(:__var1 => :($alhs), :__var2 => :(size($var1))) ) - push!(defspreamble, exx.args...) - push!(v_preamb, alhs) - continue - # - else + elseif isindexed(arhs) + # `alhs` is an aux indexed var, so something in `arhs` is indexed. + in(alhs, bkkeep.v_preamb) && continue - # `alhs` is not indexed nor an aux indexed var - in(alhs, v_preamb) && continue - vars_indexed = findex(:(_X[_i...]), arhs) - if !isempty(vars_indexed) - ex = subs(ex, Dict(vv => :(one(S)) for vv in vars_indexed)) - end - ex = subs(ex, d_decl) - push!(defspreamble, ex) - push!(v_preamb, alhs) - continue - # + vars_indexed = findex(:(_X[_i...]), arhs) + + # NOTE: Uses the size of the var with more indices + # to define the declaration of the new array. + ii_indx = argmax( [length(find_indices(aa)[1]) for aa in vars_indexed] ) + var1 = vars_indexed[ii_indx].args[1] + indx1 = vars_indexed[ii_indx].args[2:end] + + exx_indx = ones(Int, length(indx1)) + exx = subs(copy(_DECL_ARRAY), + Dict(:__var1 => :($alhs), :__var2 => :(size($var1))) ) + + # Bookkeeping + inbookkeeping(alhs, bkkeep) && alhs ∈ bkkeep.v_newvars && + deleteat!(bkkeep.v_newvars, findall(x->x==alhs, bkkeep.v_newvars)) + + push!(bkkeep.d_indx, alhs => :($alhs[$(indx1...)]) ) + push!(bkkeep.d_decl, alhs => :($alhs[$(exx_indx...)])) + push!(bkkeep.v_arraydecl, alhs) + push!(bkkeep.v_preamb, alhs) + + append!(defs_alloc, exx.args) + continue + + else + # `alhs` is not indexed nor an aux indexed var + in(alhs, bkkeep.v_preamb) && continue + vars_indexed = findex(:(_X[_i...]), arhs) + if !isempty(vars_indexed) + ex = subs(ex, Dict(vv => :(one(S)) for vv in vars_indexed)) end - # + ex = subs(ex, bkkeep.d_decl) + push!(defspreamble, ex) + push!(bkkeep.v_preamb, alhs) + continue end - else - throw(ArgumentError("$ex is not an `Expr`; $(typeof(ex))")) - # - end - exx = subs(ex, d_indx) - !inloop && push!(defspreamble, exx) + # # CHECK: Is the following needed? + # else + # # If this is reached, other `ex.head` was encountered (e.g., `.=`); + # # we include this in the `defs_tmparrays` + # inloop || (push!(defs_tmparrays, subs(ex, bkkeep.d_decl)); continue) + end + exx = subs(ex, bkkeep.d_indx) + inloop || push!(defspreamble, exx) end + + # Include allocations first in `defspreamble` + pushfirst!(defspreamble, defs_alloc...) + return defspreamble end - """ -`_recursionloop(fnargs, retvar)` +`_recursionloop(fnargs, bkkeep)` Build the expression for the recursion-loop. """ -function _recursionloop(fnargs, retvar) - +function _recursionloop(fnargs, bkkeep::BookKeeping) ll = length(fnargs) + if ll == 3 + rec_preamb = sanitize( :( $(fnargs[1]).coeffs[2] = $(bkkeep.retvar).coeffs[1] ) ) + rec_fnbody = sanitize( :( $(fnargs[1]).coeffs[ordnext+1] = $(bkkeep.retvar).coeffs[ordnext]/ordnext ) ) - rec_preamb = sanitize( :( $(fnargs[1])[1] = $(retvar)[0] ) ) - rec_fnbody = sanitize( :( $(fnargs[1])[ordnext] = $(retvar)[ord]/ordnext ) ) - # elseif ll == 4 - - retvar = fnargs[1] + bkkeep.retvar = fnargs[1] rec_preamb = sanitize(:( for __idx in eachindex($(fnargs[2])) - $(fnargs[2])[__idx].coeffs[2] = $(retvar)[__idx].coeffs[1] - end - )) + $(fnargs[2])[__idx].coeffs[2] = $(bkkeep.retvar)[__idx].coeffs[1] + end)) rec_fnbody = sanitize(:( for __idx in eachindex($(fnargs[2])) $(fnargs[2])[__idx].coeffs[ordnext+1] = - $(retvar)[__idx].coeffs[ordnext]/ordnext - end - )) - # - else + $(bkkeep.retvar)[__idx].coeffs[ordnext]/ordnext + end)) + else throw(ArgumentError( "Wrong number of arguments in the definition of the function $fn")) end + return rec_preamb, rec_fnbody end @@ -862,9 +1069,10 @@ See the [documentation](@ref taylorize) for more details and limitations. """ macro taylorize( ex ) - nex = _make_parsed_jetcoeffs(ex) + nex1, nex2 = _make_parsed_jetcoeffs(ex) esc(quote - $ex # evals to calling scope the passed function - $nex # evals the new method of `TaylorIntegration.jetcoeffs!` + $(ex) # evals to calling scope the passed function + $(nex1) # evals the new method of `jetcoeffs!` + $(nex2) # evals the new method of `_allocate_jetcoeffs` end) end diff --git a/src/rootfinding.jl b/src/rootfinding.jl index abd67bf3e..495fd0619 100644 --- a/src/rootfinding.jl +++ b/src/rootfinding.jl @@ -185,7 +185,7 @@ function taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T, sign_tstep = copysign(1, tmax-t0) # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params) # Some auxiliary arrays for root-finding/event detection/Poincaré surface of section evaluation g_tupl = g(dx, x, params, t) @@ -215,7 +215,7 @@ function taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T, nevents = 1 #number of detected events while sign_tstep*t0 < sign_tstep*tmax δt_old = δt - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive! + δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) evaluate!(x, δt, x0) # new initial condition @@ -288,7 +288,7 @@ function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T}, @inbounds xv[:,1] .= q0 # Determine if specialized jetcoeffs! method exists - parse_eqs = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params) # Some auxiliary arrays for root-finding/event detection/Poincaré surface of section evaluation g_tupl = g(dx, x, params, t) @@ -318,7 +318,7 @@ function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T}, nevents = 1 #number of detected events while sign_tstep*t0 < sign_tstep*tmax δt_old = δt - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive! + δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! # Below, δt has the proper sign according to the direction of the integration δt = sign_tstep * min(δt, sign_tstep*(tmax-t0)) evaluate!(x, δt, x0) # new initial condition diff --git a/test/TestPkg/src/TestPkg.jl b/test/TestPkg/src/TestPkg.jl index cc444b89c..4f97058ac 100644 --- a/test/TestPkg/src/TestPkg.jl +++ b/test/TestPkg/src/TestPkg.jl @@ -129,30 +129,30 @@ function TaylorIntegration.jetcoeffs!(::Val{f2!}, t::Taylor1{_T}, q::AbstractVec end ex1 = :(function f3!(dq, q, params, t) -local N = Int(length(q)/2) -local _eltype_q_ = eltype(q) -local μ = params -X = Array{_eltype_q_}(undef, N, N) -accX = Array{_eltype_q_}(undef, N) #acceleration -for j in 1:N - accX[j] = zero(q[1]) - dq[j] = q[N+j] -end -#compute accelerations -for j in 1:N + local N = Int(length(q)/2) + local _eltype_q_ = eltype(q) + local μ = params + X = Array{_eltype_q_}(undef, N, N) + accX = Array{_eltype_q_}(undef, N) #acceleration + for j in 1:N + accX[j] = zero(q[1]) + dq[j] = q[N+j] + end + #compute accelerations + for j in 1:N + for i in 1:N + if i == j + else + X[i,j] = q[i]-q[j] + temp_001 = accX[j] + (μ[i]*X[i,j]) + accX[j] = temp_001 + end #if i != j + end #for, i + end #for, j for i in 1:N - if i == j - else - X[i,j] = q[i]-q[j] - temp_001 = accX[j] + (μ[i]*X[i,j]) - accX[j] = temp_001 - end #if i != j - end #for, i -end #for, j -for i in 1:N - dq[N+i] = accX[i] -end -nothing + dq[N+i] = accX[i] + end + nothing end) ex2 = :(xdot2(x, p, t) = b1-x^2) @@ -165,9 +165,9 @@ ex3 = :(function harm_osc!(dx, x, p, t) return nothing end) -nex1 = TaylorIntegration._make_parsed_jetcoeffs(ex1) -nex2 = TaylorIntegration._make_parsed_jetcoeffs(ex2) -nex3 = TaylorIntegration._make_parsed_jetcoeffs(ex3) +nex1, narr1 = TaylorIntegration._make_parsed_jetcoeffs(ex1) +nex2, narr2 = TaylorIntegration._make_parsed_jetcoeffs(ex2) +nex3, narr3 = TaylorIntegration._make_parsed_jetcoeffs(ex3) greet(f, parse_eqs) = begin t = Taylor1(order) diff --git a/test/bigfloats.jl b/test/bigfloats.jl index 5ab494533..15bd30a53 100644 --- a/test/bigfloats.jl +++ b/test/bigfloats.jl @@ -1,9 +1,13 @@ using TaylorIntegration, Elliptic using Test using LinearAlgebra: norm +using Logging +import Logging: Warn @testset "Testing `bigfloats.jl`" begin + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + local _order = 90 local _abstol = 1.0E-77 @@ -18,7 +22,8 @@ using LinearAlgebra: norm # T is the pendulum's librational period == 4Elliptic.K(sin(q0[1]/2)^2) # we will evaluate the elliptic integral K using TaylorIntegration.jl: G(x, p, t) = 1/sqrt(1-((sin(big"1.3"/2))^2)*(sin(t)^2)) # K elliptic integral kernel - tvk, xvk = taylorinteg(G, 0.0, 0.0, BigFloat(π)/2, _order, _abstol) + tvk, xvk = (@test_logs min_level=Logging.Warn taylorinteg( + G, 0.0, 0.0, BigFloat(π)/2, _order, _abstol)) @test eltype(tvk) == BigFloat @test eltype(xvk) == BigFloat T = 4xvk[end] # T = 4Elliptic.K(sin(q0[1]/2)^2) @@ -27,7 +32,8 @@ using LinearAlgebra: norm t0 = 0.0 #the initial time - tv, xv = taylorinteg(pendulum!, q0, t0, T, _order, _abstol; maxsteps=1) + tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, q0, t0, T, _order, _abstol; maxsteps=1)) @test eltype(tv) == BigFloat @test eltype(xv) == BigFloat @test length(tv) == 2 @@ -35,7 +41,8 @@ using LinearAlgebra: norm @test length(xv[:,2]) == 2 #note that T is a BigFloat - tv, xv = taylorinteg(pendulum!, q0, t0, T, _order, _abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0, t0, T, _order, _abstol)) @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 diff --git a/test/common.jl b/test/common.jl index d567dcbe2..696af3b7e 100644 --- a/test/common.jl +++ b/test/common.jl @@ -1,57 +1,63 @@ using TaylorIntegration, Test, DiffEqBase using LinearAlgebra: norm using StaticArrays +using Logging +import Logging: Warn @testset "Testing `common.jl`" begin + local TI = TaylorIntegration + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + f(u,p,t) = u g(u,p,t) = cos(t) + f!(du, u, p, t) = (du .= u) + @testset "Test integration of ODE with numbers in common interface" begin u0 = 0.5 tspan = (0.0, 1.0) prob = ODEProblem(f, u0, tspan) - sol = solve(prob, TaylorMethod(50), abstol=1e-20) - @test TaylorIntegration.alg_order(TaylorMethod(50)) == 50 - @test TaylorIntegration.alg_order(sol.alg) == TaylorIntegration.alg_order(TaylorMethod(50)) + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(50), abstol=1e-20)) + @test TI.alg_order(TaylorMethod(50)) == 50 + @test TI.alg_order(sol.alg) == TI.alg_order(TaylorMethod(50)) @test abs(sol[end] - u0*exp(1)) < 1e-12 u0 = 0.0 tspan = (0.0, 11pi) prob = ODEProblem(g, u0, tspan) - sol = solve(prob, TaylorMethod(50), abstol=1e-20) + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(50), abstol=1e-20)) @test abs(sol[end] - sin(sol.t[end])) < 1e-12 end - f!(du, u, p, t) = (du .= u) - @testset "Test integration of ODE with abstract arrays in common interface" begin - u0 = rand(4, 2) - tspan = (0.0, 1.0) - prob = ODEProblem(f!, u0, tspan) - sol = solve(prob, TaylorMethod(50), abstol=1e-20) - @test norm(sol[end] - u0.*exp(1)) < 1e-12 - f_oop(u, p, t) = u - prob_oop = ODEProblem(f_oop, u0, tspan) - sol_oop = solve(prob_oop, TaylorMethod(50), abstol=1e-20) - @test norm(sol_oop[end] - u0.*exp(1)) < 1e-12 - @test sol.t == sol_oop.t - @test sol.u == sol_oop.u - end - - tspan = (0.0,5.0) - saveat_inputs = ([], 0:1:(tspan[2]+5), 0:1:tspan[2], 3:1:tspan[2], collect(0:1:tspan[2])) + # @testset "Test integration of ODE with abstract arrays in common interface" begin + # u0 = rand(4, 2) + # tspan = (0.0, 1.0) + # prob = ODEProblem(f!, u0, tspan) + # sol = solve(prob, TaylorMethod(50), abstol=1e-20) + # @test norm(sol[end] - u0.*exp(1)) < 1e-12 + # f_oop(u, p, t) = u + # prob_oop = ODEProblem(f_oop, u0, tspan) + # sol_oop = solve(prob_oop, TaylorMethod(50), abstol=1e-20) + # @test norm(sol_oop[end] - u0.*exp(1)) < 1e-12 + # @test sol.t == sol_oop.t + # @test sol.u == sol_oop.u + # end + + local tspan = (0.0,5.0) + local saveat_inputs = ([], 0:1:(tspan[2]+5), 0:1:tspan[2], 3:1:tspan[2], collect(0:1:tspan[2])) @testset "Test saveat behavior with numbers in common interface" begin u0 = 1.0 prob = ODEProblem(f, u0, tspan) - sol = solve(prob, TaylorMethod(20), abstol=1e-20) + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20)) s = saveat_inputs[1] - sol_taylor = solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s) + sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s)) @test all(sol.t .== sol_taylor.t) @test all(sol_taylor.u .== sol.u) @test length(sol_taylor.t) == length(sol_taylor.u) s = saveat_inputs[2] - sol_taylor = solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s) + sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s)) @test sol_taylor.t == saveat_inputs[3] for s ∈ saveat_inputs[3:end] - sol_taylor = solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s) + sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s)) @test all(s .== sol_taylor.t) @test length(sol_taylor.t) == length(sol_taylor.u) end @@ -60,17 +66,17 @@ using StaticArrays @testset "Test saveat behavior with abstract arrays in common interface" begin u0 = rand(2) prob = ODEProblem(f!, u0, tspan) - sol = solve(prob, TaylorMethod(20), abstol=1e-20) + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20)) s = saveat_inputs[1] - sol_taylor = solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s) + sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s)) @test all(sol.t .== sol_taylor.t) @test all(sol_taylor.u .== sol.u) @test length(sol_taylor.t) == length(sol_taylor.u) s = saveat_inputs[2] - sol_taylor = solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s) + sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s)) @test sol_taylor.t == saveat_inputs[3] for s ∈ saveat_inputs[3:end] - sol_taylor = solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s) + sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s)) @test all(s .== sol_taylor.t) @test length(sol_taylor.t) == length(sol_taylor.u) end @@ -81,14 +87,15 @@ using StaticArrays dx[2] = - x[1] return nothing end - tspan = (0.0, 10pi) - abstol = 1e-20 - order = 25 # Taylor expansion order wrt time - u0 = [1.0; 0.0] - prob = ODEProblem(harmosc!, u0, tspan) @testset "Test consistency with taylorinteg" begin - sol = solve(prob, TaylorMethod(order), abstol=abstol) - tv1, xv1 = taylorinteg(harmosc!, u0, tspan[1], tspan[2], order, abstol) + tspan = (0.0, 10pi) + abstol = 1e-20 + order = 25 # Taylor expansion order wrt time + u0 = [1.0; 0.0] + prob = ODEProblem(harmosc!, u0, tspan) + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(order), abstol=abstol)) + tv1, xv1 = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, u0, tspan[1], tspan[2], order, abstol)) @test sol.t == tv1 @test xv1[end,:] == sol[end] end @@ -96,11 +103,16 @@ using StaticArrays @testset "Test discrete callback in common interface" begin # discrete callback: example taken from DifferentialEquations.jl docs: # https://diffeq.sciml.ai/dev/features/callback_functions/#Using-Callbacks + tspan = (0.0, 10pi) + abstol = 1e-20 + order = 25 # Taylor expansion order wrt time + u0 = [1.0; 0.0] + prob = ODEProblem(harmosc!, u0, tspan) t_cb = 1.0pi condition(u,t,integrator) = t == t_cb affect!(integrator) = integrator.u[1] += 0.1 cb = DiscreteCallback(condition,affect!) - sol = solve(prob, TaylorMethod(order), abstol=abstol, tstops=[t_cb], callback=cb) + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(order), abstol=abstol, tstops=[t_cb], callback=cb)) @test sol.t[4] == t_cb @test sol.t[4] == sol.t[5] @test sol[4][1] + 0.1 == sol[5][1] @@ -126,7 +138,12 @@ using StaticArrays tspan = (0.0,15.0) p = 9.8 prob = ODEProblem(f,u0,tspan,p) - sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb) + # Avoid log-checks for Julia versions before v1.6 + if VERSION < v"1.6" + sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb) + else + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb)) + end tb = sqrt(2*50/9.8) # bounce time @test abs(tb - sol.t[9]) < 1e-14 @test sol.t[9] == sol.t[10] @@ -141,7 +158,7 @@ using StaticArrays @testset "Test vector continuous callback in common interface" begin # vector continuous callback: example taken from DifferentialEquations.jl docs: # https://diffeq.sciml.ai/dev/features/callback_functions/#VectorContinuousCallback-Example - @taylorize function f(du,u,p,t) + @taylorize function ff(du,u,p,t) local g_acc = p du[1] = u[2] du[2] = -g_acc + zero(u[2]) @@ -167,8 +184,13 @@ using StaticArrays u0 = [50.0, 0.0, 0.0, 2.0] tspan = (0.0, 15.0) p = 9.8 - prob = ODEProblem(f, u0, tspan, p) - sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb) + prob = ODEProblem(ff, u0, tspan, p) + # Avoid log-checks for Julia versions before v1.6 + if VERSION < v"1.6" + sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb) + else + sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb)) + end tb = sqrt(2*50/9.8) # bounce time @test abs(tb - sol.t[8]) < 1e-14 @test sol.t[8] == sol.t[9] @@ -185,14 +207,18 @@ using StaticArrays @testset "Test parsed jetcoeffs! method in common interface" begin b1 = 3.0 + order = 25 # Taylor expansion order wrt time + abstol = 1.0e-20 @taylorize xdot1(x, p, t) = b1-x^2 @test (@isdefined xdot1) x0 = 1.0 tspan = (0.0, 1000.0) prob = ODEProblem(xdot1, x0, tspan) - sol1 = solve(prob, TaylorMethod(order), abstol=abstol, parse_eqs=false) + sol1 = (@test_logs min_level=Logging.Warn solve( + prob, TaylorMethod(order), abstol=abstol, parse_eqs=false)) @test sol1.alg.parse_eqs == false - sol2 = solve(prob, TaylorMethod(order), abstol=abstol) + sol2 = (@test_logs min_level=Logging.Warn solve( + prob, TaylorMethod(order), abstol=abstol)) @test sol2.alg.parse_eqs == true @test sol1.t == sol2.t @test sol1.u == sol2.u @@ -207,13 +233,16 @@ using StaticArrays x0 = [0.0, 1.0] tspan = (0.0, 11pi) prob = ODEProblem(integ_vec, x0, tspan, [1.0]) - sol1 = solve(prob, TaylorMethod(order), abstol=abstol, parse_eqs=false) - sol2 = solve(prob, TaylorMethod(order), abstol=abstol) # parse_eqs=true + sol1 = (@test_logs min_level=Logging.Warn solve( + prob, TaylorMethod(order), abstol=abstol, parse_eqs=false)) + sol2 = (@test_logs min_level=Logging.Warn solve( + prob, TaylorMethod(order), abstol=abstol)) # parse_eqs=true @test sol2.alg.parse_eqs == true @test length(sol1.t) == length(sol2.t) @test sol1.t == sol2.t @test sol1.u == sol2.u - tv, xv = taylorinteg(integ_vec, x0, tspan[1], tspan[2], order, abstol, [1.0]) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + integ_vec, x0, tspan[1], tspan[2], order, abstol, [1.0])) @test sol1.t == tv @test sol1[1,:] == xv[:,1] @test sol1[2,:] == xv[:,2] @@ -246,9 +275,11 @@ using StaticArrays return nothing end prob = ODEProblem(pcr3bp!, q0, tspan, p) - sol1 = solve(prob, TaylorMethod(order), abstol=abstol) + sol1 = (@test_logs min_level=Logging.Warn solve( + prob, TaylorMethod(order), abstol=abstol)) @test sol1.alg.parse_eqs == true - sol2 = solve(prob, TaylorMethod(order), abstol=abstol, parse_eqs=false) + sol2 = (@test_logs min_level=Logging.Warn solve( + prob, TaylorMethod(order), abstol=abstol, parse_eqs=false)) @test sol2.alg.parse_eqs == false @test norm( H_pcr3bp(sol1.u[end]) - J0 ) < 1e-10 @test norm( H_pcr3bp(sol2.u[end]) - J0 ) < 1e-10 @@ -257,10 +288,11 @@ using StaticArrays end @testset "Test throwing errors in common interface" begin - u0 = rand(4, 2) + # NOTE: The original test defines u0 as a matrix (rand(4,2)) but that doesn't work + u0 = rand(4)#rand(4, 2) tspan = (0.0, 1.0) prob1 = ODEProblem(f!, u0, tspan) - sol = solve(prob1, TaylorMethod(50), abstol=1e-20) + sol = (@test_logs min_level=Logging.Warn solve(prob1, TaylorMethod(50), abstol=1e-20)) # `order` is not specified @test_throws ErrorException solve(prob, TaylorMethod(), abstol=1e-20) @@ -306,14 +338,14 @@ using StaticArrays iip_prob = DynamicalODEProblem(iip_ṗ, iip_q̇, iip_p0, iip_q0, (0., 100.)) oop_prob = DynamicalODEProblem(oop_ṗ, oop_q̇, oop_p0, oop_q0, (0., 100.)) - sol1 = solve(iip_prob, TaylorMethod(50), abstol=1e-20) + sol1 = (@test_logs min_level=Logging.Warn solve(iip_prob, TaylorMethod(50), abstol=1e-20)) @test energy_err(sol1) < 1e-10 - sol2 = solve(oop_prob, TaylorMethod(50), abstol=1e-20) + sol2 = (@test_logs min_level=Logging.Warn solve(oop_prob, TaylorMethod(50), abstol=1e-20)) @test energy_err(sol2) < 1e-10 @test sol1.t == sol2.t - @test sol1.u == sol2.u + @test sol1.u[:] == sol2.u[:] end end diff --git a/test/complex.jl b/test/complex.jl index e3a1f4339..60bef45a3 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -2,12 +2,16 @@ using TaylorIntegration using Test +using Logging +import Logging: Warn @testset "Testing `complex.jl`" begin local _order = 28 local _abstol = 1.0E-20 + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + eqs1(z, p, t) = -z eqs2(z, p, t) = im*z function eqs3!(Dz, z, p, t) @@ -20,28 +24,35 @@ using Test z0 = complex(0.0, 1.0) tr = 0.0:pi/8:2pi ts = 0.0:pi:2pi - zsol1 = taylorinteg(eqs1, z0, ts, _order, _abstol, maxsteps=1) + zsol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs1, z0, ts, _order, _abstol, maxsteps=1)) @test size(zsol1) == (length(ts),) - zsol1 = taylorinteg(eqs1, z0, tr, _order, _abstol, maxsteps=1, nothing) + zsol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs1, z0, tr, _order, _abstol, maxsteps=1, nothing)) @test length(zsol1) == length(tr) ta = vec(tr) - zsol1 = taylorinteg(eqs1, z0, ta, _order, _abstol, maxsteps=1) + zsol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs1, z0, ta, _order, _abstol, maxsteps=1)) @test length(zsol1) == length(ta) - zsol1 = taylorinteg(eqs1, z0, tr, _order, _abstol) + zsol1 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs1, z0, tr, _order, _abstol)) @test zsol1[1] == z0 @test isapprox( zsol1[2] , z0*exp(-tr[2]) ) @test isapprox( zsol1[6], z0*exp(-tr[6]) ) @test isapprox( zsol1[end], z0*exp(-tr[end]) ) - zsol11 = taylorinteg(eqs1, z0, ta, _order, _abstol, nothing) + zsol11 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs1, z0, ta, _order, _abstol, nothing)) @test zsol11 == zsol1 @test zsol11[1] == z0 @test isapprox( zsol11[2] , z0*exp(-tr[2]) ) @test isapprox( zsol11[6], z0*exp(-tr[6]) ) @test isapprox( zsol11[end], z0*exp(-tr[end]) ) - tt, zsol1t = taylorinteg(eqs1, z0, 0.0, 2pi, _order, _abstol, maxsteps=1) + tt, zsol1t = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs1, z0, 0.0, 2pi, _order, _abstol, maxsteps=1)) @test length(tt) == 2 @test length(zsol1t) == 2 - tt, zsol1t = taylorinteg(eqs1, z0, 0.0, 2pi, _order, _abstol) + tt, zsol1t = (@test_logs min_level=Logging.Warn taylorinteg( + eqs1, z0, 0.0, 2pi, _order, _abstol)) @test zsol1t[1] == z0 @test isapprox( zsol1t[2] , z0*exp(-tt[2]) ) @test isapprox( zsol1t[end], z0*exp(-tt[end]) ) @@ -52,26 +63,33 @@ using Test z0 = complex(0.0, 1.0) tr = 0.0:pi/8:2pi ts = 0.0:pi:2pi - zsol2 = taylorinteg(eqs2, z0, ts, _order, _abstol, maxsteps=1, nothing) + zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs2, z0, ts, _order, _abstol, maxsteps=1, nothing)) @test size(zsol2) == (length(ts),) - zsol2 = taylorinteg(eqs2, z0, tr, _order, _abstol, maxsteps=1) + zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs2, z0, tr, _order, _abstol, maxsteps=1)) @test length(zsol2) == length(tr) ta = vec(tr) - zsol2 = taylorinteg(eqs2, z0, ta, _order, _abstol, maxsteps=1) + zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs2, z0, ta, _order, _abstol, maxsteps=1)) @test length(zsol2) == length(ta) - zsol2 = taylorinteg(eqs2, z0, tr, _order, _abstol) + zsol2 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs2, z0, tr, _order, _abstol)) @test zsol2[1] == z0 @test isapprox( zsol2[3], z0*exp( complex(0.0, tr[3])) ) @test isapprox( zsol2[5], z0*exp( complex(0.0, tr[5])) ) - zsol22 = taylorinteg(eqs2, z0, ta, _order, _abstol, nothing) + zsol22 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs2, z0, ta, _order, _abstol, nothing)) @test zsol22 == zsol2 @test zsol22[1] == z0 @test isapprox( zsol22[3], z0*exp( complex(0.0, tr[3])) ) @test isapprox( zsol22[5], z0*exp( complex(0.0, tr[5])) ) - tt, zsol2 = taylorinteg(eqs2, z0, 0.0, 2pi, _order, _abstol, maxsteps=1) + tt, zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs2, z0, 0.0, 2pi, _order, _abstol, maxsteps=1)) @test length(tt) == 2 @test length(zsol2) == 2 - tt, zsol2 = taylorinteg(eqs2, z0, 0.0, 2pi, _order, _abstol) + tt, zsol2 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs2, z0, 0.0, 2pi, _order, _abstol)) @test zsol2[1] == z0 @test isapprox( zsol2[2], z0*exp( complex(0.0, tt[2])) ) @test isapprox( zsol2[end], z0*exp( complex(0.0, tt[end])) ) @@ -82,28 +100,35 @@ using Test zz0 = [z0, z0] tr = 0.0:pi/8:2pi ts = 0.0:pi:2pi - zsol3 = taylorinteg(eqs3!, zz0, ts, _order, _abstol, maxsteps=1) + zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs3!, zz0, ts, _order, _abstol, maxsteps=1)) @test size(zsol3) == (length(ts), length(zz0)) - zsol3 = taylorinteg(eqs3!, zz0, tr, _order, _abstol, nothing, maxsteps=1) + zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs3!, zz0, tr, _order, _abstol, nothing, maxsteps=1)) @test size(zsol3) == (length(tr), length(zz0)) ta = vec(tr) - zsol3 = taylorinteg(eqs3!, zz0, ta, _order, _abstol, maxsteps=1) + zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs3!, zz0, ta, _order, _abstol, maxsteps=1)) @test size(zsol3) == (length(ta), length(zz0)) - zsol3 = taylorinteg(eqs3!, [z0, z0], tr, _order, _abstol) + zsol3 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs3!, [z0, z0], tr, _order, _abstol)) @test isapprox( zsol3[4,1], z0*exp(-tr[4]) ) @test isapprox( zsol3[7,1], z0*exp(-tr[7]) ) @test isapprox( zsol3[4,2], z0*exp( complex(0.0, tr[4])) ) @test isapprox( zsol3[7,2], z0*exp( complex(0.0, tr[7])) ) - zsol33 = taylorinteg(eqs3!, [z0, z0], ta, _order, _abstol, nothing) + zsol33 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs3!, [z0, z0], ta, _order, _abstol, nothing)) @test zsol33 == zsol3 @test isapprox( zsol33[4,1], z0*exp(-tr[4]) ) @test isapprox( zsol33[7,1], z0*exp(-tr[7]) ) @test isapprox( zsol33[4,2], z0*exp( complex(0.0, tr[4])) ) @test isapprox( zsol33[7,2], z0*exp( complex(0.0, tr[7])) ) - tt, zsol3 = taylorinteg(eqs3!, zz0, 0.0, 2pi, _order, _abstol, maxsteps=1) + tt, zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs3!, zz0, 0.0, 2pi, _order, _abstol, maxsteps=1)) @test length(tt) == 2 @test size(zsol3) == (length(tt), length(zz0)) - tt, zsol3 = taylorinteg(eqs3!, zz0, 0.0, 2pi, _order, _abstol) + tt, zsol3 = (@test_logs min_level=Logging.Warn taylorinteg( + eqs3!, zz0, 0.0, 2pi, _order, _abstol)) @test zsol3[1,1] == z0 @test zsol3[1,2] == z0 @test isapprox( zsol3[2,1], z0*exp( -tt[2]) ) diff --git a/test/jettransport.jl b/test/jettransport.jl index e9c479fde..670562f21 100644 --- a/test/jettransport.jl +++ b/test/jettransport.jl @@ -3,12 +3,16 @@ using TaylorIntegration, Elliptic using LinearAlgebra: norm using Test +using Logging +import Logging: Warn @testset "Testing `jettransport.jl`" begin local _order = 28 local _abstol = 1.0E-20 + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + f(x, p, t) = x^2 g(x, p, t) = 0.3x @@ -18,11 +22,14 @@ using Test x0TN = x0 + p[1] #jet transport initial condition t0=0.0 tmax=0.3 - tvTN, xvTN = taylorinteg(f, x0TN, t0, tmax, _order, _abstol, maxsteps=1) + tvTN, xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + f, x0TN, t0, tmax, _order, _abstol, maxsteps=1)) @test size(tvTN) == (2,) @test size(xvTN) == (2,) - tvTN, xvTN = taylorinteg(f, x0TN, t0, tmax, _order, _abstol) - tv, xv = taylorinteg(f, x0, t0, tmax, _order, _abstol) + tvTN, xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0TN, t0, tmax, _order, _abstol)) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0, t0, tmax, _order, _abstol)) exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution δsol = exactsol(tvTN[end], x0TN, t0)-xvTN[end] δcoeffs = map(y->y[1], map(x->x.coeffs, δsol.coeffs)) @@ -36,7 +43,8 @@ using Test disp = 0.001*rand() #a small, random displacement x0_disp = x0+disp dv = map(x->[disp], tvTN) #a vector of identical displacements - tv_disp, xv_disp = taylorinteg(f, x0_disp, t0, tmax, _order, _abstol) + tv_disp, xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0_disp, t0, tmax, _order, _abstol)) xvTN_disp = evaluate.(xvTN, dv) @test norm(exactsol.(tvTN, x0_disp, t0)-xvTN_disp, Inf) < 1E-12 #analytical vs jet transport @test norm(x0_disp-evaluate(xvTN[1], [disp]), Inf) < 1E-12 @@ -47,11 +55,14 @@ using Test y0 = 1.0 #"nominal" initial condition u0 = 0.0 #initial time y0TN = y0 + p[1] #jet transport initial condition - uvTN, yvTN = taylorinteg(g, y0TN, u0, 10/0.3, _order, _abstol, maxsteps=1) + uvTN, yvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + g, y0TN, u0, 10/0.3, _order, _abstol, maxsteps=1)) @test size(uvTN) == (2,) @test size(yvTN) == (2,) - uvTN, yvTN = taylorinteg(g, y0TN, u0, 10/0.3, _order, _abstol) - uv, yv = taylorinteg(g, y0, u0, 10/0.3, _order, _abstol) + uvTN, yvTN = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0TN, u0, 10/0.3, _order, _abstol)) + uv, yv = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0, u0, 10/0.3, _order, _abstol)) exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) δsol_g = exactsol_g(uvTN[end], y0TN, u0)-yvTN[end] δcoeffs_g = map(y->y[1], map(x->x.coeffs, δsol_g.coeffs)) @@ -65,7 +76,8 @@ using Test disp = 0.001*rand() #a small, random displacement y0_disp = y0+disp dv = map(x->[disp], uvTN) #a vector of identical displacements - uv_disp, yv_disp = taylorinteg(g, y0_disp, u0, 10/0.3, _order, _abstol) + uv_disp, yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0_disp, u0, 10/0.3, _order, _abstol)) yvTN_disp = evaluate.(yvTN, dv) @test norm(exactsol_g.(uvTN, y0_disp, u0)-yvTN_disp, Inf) < 1E-9 #analytical vs jet transport @test norm(y0_disp-yvTN[1]([disp]), Inf) < 1E-9 @@ -80,11 +92,14 @@ using Test x0T1 = x0 + p #jet transport initial condition t0 = 0.0 tmax = 0.3 - tvT1, xvT1 = taylorinteg(f, x0T1, t0, tmax, _order, _abstol, maxsteps=1) + tvT1, xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + f, x0T1, t0, tmax, _order, _abstol, maxsteps=1)) @test size(tvT1) == (2,) @test size(xvT1) == (2,) - tvT1, xvT1 = taylorinteg(f, x0T1, t0, tmax, _order, _abstol) - tv, xv = taylorinteg(f, x0, t0, tmax, _order, _abstol) + tvT1, xvT1 = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0T1, t0, tmax, _order, _abstol)) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0, t0, tmax, _order, _abstol)) exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution δsol = exactsol(tvT1[end], x0T1, t0)-xvT1[end] #analytical vs jet transport diff at end of integration δcoeffs = δsol.coeffs @@ -97,7 +112,8 @@ using Test for i in 1:5 disp = 0.001*rand() #a small, random displacement x0_disp = x0+disp - tv_disp, xv_disp = taylorinteg(f, x0_disp, t0, tmax, _order, _abstol) + tv_disp, xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0_disp, t0, tmax, _order, _abstol)) xvT1_disp = evaluate.(xvT1, disp) @test norm(exactsol.(tvT1, x0_disp, t0)-xvT1_disp, Inf) < 1E-12 #analytical vs jet transport @test norm(x0_disp-evaluate(xvT1[1], disp), Inf) < 1E-12 @@ -108,12 +124,15 @@ using Test y0 = 1.0 #"nominal" initial condition u0 = 0.0 #initial time y0T1 = y0 + p #jet transport initial condition - uvT1, yvT1 = taylorinteg(g, y0T1, t0, 10/0.3, _order, _abstol, maxsteps=1) #warmup lap + uvT1, yvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + g, y0T1, t0, 10/0.3, _order, _abstol, maxsteps=1)) #warmup lap # test maxsteps break @test size(uvT1) == (2,) @test size(yvT1) == (2,) - uvT1, yvT1 = taylorinteg(g, y0T1, t0, 10/0.3, _order, _abstol) #Taylor1 jet transport integration - uv, yv = taylorinteg(g, y0, t0, 10/0.3, _order, _abstol) #reference integration + uvT1, yvT1 = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0T1, t0, 10/0.3, _order, _abstol)) #Taylor1 jet transport integration + uv, yv = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0, t0, 10/0.3, _order, _abstol)) #reference integration exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) #the analytical solution δsol_g = exactsol_g(uvT1[end], y0T1, t0)-yvT1[end] #analytical vs jet transport diff at end of integration δcoeffs_g = δsol_g.coeffs @@ -126,7 +145,8 @@ using Test for i in 1:5 disp = 0.001*rand() #a small, random displacement y0_disp = y0+disp - uv_disp, yv_disp = taylorinteg(g, y0_disp, u0, 10/0.3, _order, _abstol) + uv_disp, yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0_disp, u0, 10/0.3, _order, _abstol)) yvT1_disp = evaluate.(yvT1, disp) @test norm(exactsol_g.(uvT1, y0_disp, u0)-yvT1_disp, Inf) < 1E-9 #analytical vs jet transport @test norm(y0_disp-evaluate(yvT1[1], disp), Inf) < 1E-9 @@ -140,15 +160,20 @@ using Test x0 = 3.0 #"nominal" initial condition x0T1 = x0 + p #jet transport initial condition tv = 0.0:0.05:0.33 - xvT1 = taylorinteg(f, x0T1, tv, _order, _abstol, maxsteps=1) + xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + f, x0T1, tv, _order, _abstol, maxsteps=1)) @test size(xvT1) == (7,) ta = vec(tv) - xvT1 = taylorinteg(f, x0T1, ta, _order, _abstol, maxsteps=1) + xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + f, x0T1, ta, _order, _abstol, maxsteps=1)) @test size(xvT1) == (7,) - xvT1 = taylorinteg(f, x0T1, tv, _order, _abstol) - xvT11 = taylorinteg(f, x0T1, ta, _order, _abstol) + xvT1 = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0T1, tv, _order, _abstol)) + xvT11 = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0T1, ta, _order, _abstol)) @test xvT11 == xvT1 - xv = taylorinteg(f, x0, tv, _order, _abstol) + xv = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0, tv, _order, _abstol)) exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution δsol = exactsol(tv[end], x0T1, tv[1])-xvT1[end] δcoeffs = δsol.coeffs @@ -159,8 +184,10 @@ using Test @test isapprox(xv_analytical, xvT1_0, atol=1e-10, rtol=0) for i in 1:5 disp = 0.001*rand() #a small, random displacement - xv_disp = taylorinteg(f, x0+disp, tv, _order, _abstol) - xv_disp2 = taylorinteg(f, x0+disp, ta, _order, _abstol) + xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0+disp, tv, _order, _abstol)) + xv_disp2 = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0+disp, ta, _order, _abstol)) @test xv_disp == xv_disp2 xvT1_disp = xvT1.(disp) @test norm(exactsol.(tv, x0+disp, tv[1])-xvT1_disp, Inf) < 1E-12 #analytical vs jet transport @@ -170,10 +197,13 @@ using Test y0 = 1.0 #"nominal" initial condition y0T1 = y0 + p #jet transport initial condition uv = 0.0:1/0.3:10/0.3 - yvT1 = taylorinteg(g, y0T1, uv, _order, _abstol, maxsteps=1) + yvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + g, y0T1, uv, _order, _abstol, maxsteps=1)) @test size(yvT1) == (11,) - yvT1 = taylorinteg(g, y0T1, uv, _order, _abstol) - yv = taylorinteg(g, y0, uv, _order, _abstol) + yvT1 = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0T1, uv, _order, _abstol)) + yv = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0, uv, _order, _abstol)) exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) δsol_g = exactsol_g(uv[end], y0T1, uv[1])-yvT1[end] δcoeffs_g = δsol_g.coeffs @@ -184,7 +214,8 @@ using Test @test isapprox(yv_analytical, yvT1_0, atol=1e-10, rtol=0) for i in 1:5 disp = 0.001*rand() #a small, random displacement - yv_disp = taylorinteg(g, y0+disp, uv, _order, _abstol) + yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0+disp, uv, _order, _abstol)) yvT1_disp = evaluate.(yvT1, disp) @test norm(exactsol_g.(uv, y0+disp, uv[1])-yvT1_disp, Inf) < 1E-9 #analytical vs jet transport @test norm(yv_disp-yvT1_disp, Inf) < 1E-9 # integration vs jet transport @@ -196,15 +227,20 @@ using Test x0 = 3.0 #"nominal" initial condition x0TN = x0 + p[1] #jet transport initial condition tv = 0.0:0.05:0.33 - xvTN = taylorinteg(f, x0TN, tv, _order, _abstol, maxsteps=1) + xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + f, x0TN, tv, _order, _abstol, maxsteps=1)) @test size(xvTN) == (7,) ta = vec(tv) - xvTN = taylorinteg(f, x0TN, ta, _order, _abstol, maxsteps=1) + xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + f, x0TN, ta, _order, _abstol, maxsteps=1)) @test size(xvTN) == (7,) - xvTN = taylorinteg(f, x0TN, tv, _order, _abstol) - xvTN_ = taylorinteg(f, x0TN, ta, _order, _abstol) + xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0TN, tv, _order, _abstol)) + xvTN_ = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0TN, ta, _order, _abstol)) @test xvTN == xvTN_ - xv = taylorinteg(f, x0, tv, _order, _abstol) + xv = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0, tv, _order, _abstol)) exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution δsol = exactsol(tv[end], x0TN, tv[1])-xvTN[end] δcoeffs = map(y->y[1], map(x->x.coeffs, δsol.coeffs)) @@ -216,8 +252,10 @@ using Test for i in 1:5 disp = 0.001*rand() #a small, random displacement dv = map(x->[disp], tv) #a vector of identical displacements - xv_disp = taylorinteg(f, x0+disp, tv, _order, _abstol) - xv_disp_ = taylorinteg(f, x0+disp, ta, _order, _abstol) + xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0+disp, tv, _order, _abstol)) + xv_disp_ = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0+disp, ta, _order, _abstol)) @test xv_disp == xv_disp_ xvTN_disp = evaluate.(xvTN, dv) @test norm(exactsol.(tv, x0+disp, tv[1])-xvTN_disp, Inf) < 1E-12 #analytical vs jet transport @@ -227,10 +265,13 @@ using Test y0 = 1.0 #"nominal" initial condition y0TN = y0 + p[1] #jet transport initial condition uv = 0.0:1/0.3:10/0.3 - yvTN = taylorinteg(g, y0TN, uv, _order, _abstol, maxsteps=1) + yvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + g, y0TN, uv, _order, _abstol, maxsteps=1)) @test size(yvTN) == (11,) - yvTN = taylorinteg(g, y0TN, uv, _order, _abstol) - yv = taylorinteg(g, y0, uv, _order, _abstol) + yvTN = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0TN, uv, _order, _abstol)) + yv = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0, uv, _order, _abstol)) exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) δsol_g = exactsol_g(uv[end], y0TN, uv[1])-yvTN[end] δcoeffs_g = map(y->y[1], map(x->x.coeffs, δsol_g.coeffs)) @@ -242,112 +283,129 @@ using Test for i in 1:5 disp = 0.001*rand() #a small, random displacement dv = map(x->[disp], uv) #a vector of identical displacements - yv_disp = taylorinteg(g, y0+disp, uv, _order, _abstol) + yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + g, y0+disp, uv, _order, _abstol)) yvTN_disp = evaluate.(yvTN, dv) @test norm(exactsol_g.(uv, y0+disp, uv[1])-yvTN_disp, Inf) < 1E-9 #analytical vs jet transport @test norm(yv_disp-yvTN_disp, Inf) < 1E-9 # integration vs jet transport end end - function harmosc!(dx, x, p, t) #the harmonic oscillator ODE - dx[1] = x[2] - dx[2] = -x[1]*x[3]^2 - dx[3] = zero(x[1]) - nothing - end + let + function harmosc!(dx, x, p, t) #the harmonic oscillator ODE + dx[1] = x[2] + dx[2] = -x[1]*x[3]^2 + dx[3] = zero(x[1]) + nothing + end - @testset "Test Taylor1 jet transport (t0,tmax): harmonic oscillator" begin - t = Taylor1([0.0, 1.0], 10) - ω0 = 1.0 - x0 = [0.0,ω0,ω0] - x0T1 = x0+[0t,t,t] - tv1, xv1 = taylorinteg(harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=1) - @test size(tv1) == (2,) - @test size(xv1) == (2, 3) - tv1, xv1 = taylorinteg(harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=2000) - y0 = evaluate.(xv1) - x1(t,δω) = sin((ω0+δω)*t) - x2(t,δω) = (ω0+δω)*cos((ω0+δω)*t) - @test norm(y0[:,1]-x1.(tv1,0.0), Inf) < 1E-11 - @test norm(y0[:,2]-x2.(tv1,0.0), Inf) < 1E-11 - for i in 1:5 - δω=0.001*rand() - x0_disp = x0+[0.0,δω,δω] - tv, xv = taylorinteg(harmosc!, x0_disp, 0.0, 100pi, _order, _abstol, maxsteps=2000) - y1 = evaluate.(xv1, δω) - @test norm(y1[:,1]-x1.(tv1,δω),Inf) < 1E-11 - @test norm(y1[:,2]-x2.(tv1,δω),Inf) < 1E-11 - @test norm(x0_disp-evaluate.(xv1[1,:], δω), Inf) < 1E-11 - @test norm(xv[1,:]-evaluate.(xv1[1,:], δω), Inf) < 1E-11 - @test norm(xv[end,:]-evaluate.(xv1[end,:], δω), Inf) < 1E-11 + @testset "Test Taylor1 jet transport (t0,tmax): harmonic oscillator" begin + t = Taylor1([0.0, 1.0], 10) + ω0 = 1.0 + x0 = [0.0,ω0,ω0] + x0T1 = x0+[0t,t,t] + tv1, xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=1)) + @test size(tv1) == (2,) + @test size(xv1) == (2, 3) + tv1, xv1 = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=2000)) + y0 = evaluate.(xv1) + x1(t,δω) = sin((ω0+δω)*t) + x2(t,δω) = (ω0+δω)*cos((ω0+δω)*t) + @test norm(y0[:,1]-x1.(tv1,0.0), Inf) < 1E-11 + @test norm(y0[:,2]-x2.(tv1,0.0), Inf) < 1E-11 + for i in 1:5 + δω=0.001*rand() + x0_disp = x0+[0.0,δω,δω] + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0_disp, 0.0, 100pi, _order, _abstol, maxsteps=2000)) + y1 = evaluate.(xv1, δω) + @test norm(y1[:,1]-x1.(tv1,δω),Inf) < 1E-11 + @test norm(y1[:,2]-x2.(tv1,δω),Inf) < 1E-11 + @test norm(x0_disp-evaluate.(xv1[1,:], δω), Inf) < 1E-11 + @test norm(xv[1,:]-evaluate.(xv1[1,:], δω), Inf) < 1E-11 + @test norm(xv[end,:]-evaluate.(xv1[end,:], δω), Inf) < 1E-11 + end end - end - @testset "Test Taylor1 jet transport (trange): harmonic oscillator" begin - t = Taylor1([0.0, 1.0], 10) - ω0 = 1.0 - x0 = [0.0,ω0,ω0] - x0T1 = x0+[0t,t,t] - tv = 0.0:0.25*(2pi):100pi - xv1 = taylorinteg(harmosc!, x0T1, tv, _order, _abstol, maxsteps=1) - @test length(tv) == 201 - @test size(xv1) == (201, 3) - ta = vec(tv) - xv1 = taylorinteg(harmosc!, x0T1, ta, _order, _abstol, maxsteps=1) - @test length(ta) == 201 - @test size(xv1) == (201, 3) - xv1 = taylorinteg(harmosc!, x0T1, tv, _order, _abstol, maxsteps=2000) - xv1_ = taylorinteg(harmosc!, x0T1, ta, _order, _abstol, maxsteps=2000) - @test xv1 == xv1_ - y0 = evaluate.(xv1) - x1(t,δω) = sin((ω0+δω)*t) - x2(t,δω) = (ω0+δω)*cos((ω0+δω)*t) - @test norm(y0[:,1]-x1.(tv,0.0), Inf) < 1E-11 - @test norm(y0[:,2]-x2.(tv,0.0), Inf) < 1E-11 - for i in 1:5 - δω=0.001*rand() - x0_disp = x0+[0.0,δω,δω] - xv = taylorinteg(harmosc!, x0_disp, tv, _order, _abstol, maxsteps=2000) - xv_ = taylorinteg(harmosc!, x0_disp, ta, _order, _abstol, maxsteps=2000) - @test xv == xv_ - y1 = evaluate.(xv1, δω) - @test norm(y1[:,1]-x1.(tv,δω), Inf) < 1E-11 - @test norm(y1[:,2]-x2.(tv,δω), Inf) < 1E-11 - @test norm(y1-xv, Inf) < 1E-11 + @testset "Test Taylor1 jet transport (trange): harmonic oscillator" begin + t = Taylor1([0.0, 1.0], 10) + ω0 = 1.0 + x0 = [0.0,ω0,ω0] + x0T1 = x0+[0t,t,t] + tv = 0.0:0.25*(2pi):100pi + xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + harmosc!, x0T1, tv, _order, _abstol, maxsteps=1)) + @test length(tv) == 201 + @test size(xv1) == (201, 3) + ta = vec(tv) + xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + harmosc!, x0T1, ta, _order, _abstol, maxsteps=1)) + @test length(ta) == 201 + @test size(xv1) == (201, 3) + xv1 = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0T1, tv, _order, _abstol, maxsteps=2000)) + xv1_ = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0T1, ta, _order, _abstol, maxsteps=2000)) + @test xv1 == xv1_ + y0 = evaluate.(xv1) + x1(t,δω) = sin((ω0+δω)*t) + x2(t,δω) = (ω0+δω)*cos((ω0+δω)*t) + @test norm(y0[:,1]-x1.(tv,0.0), Inf) < 1E-11 + @test norm(y0[:,2]-x2.(tv,0.0), Inf) < 1E-11 + for i in 1:5 + δω=0.001*rand() + x0_disp = x0+[0.0,δω,δω] + xv = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0_disp, tv, _order, _abstol, maxsteps=2000)) + xv_ = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0_disp, ta, _order, _abstol, maxsteps=2000)) + @test xv == xv_ + y1 = evaluate.(xv1, δω) + @test norm(y1[:,1]-x1.(tv,δω), Inf) < 1E-11 + @test norm(y1[:,2]-x2.(tv,δω), Inf) < 1E-11 + @test norm(y1-xv, Inf) < 1E-11 + end end end - function harmosc!(dx, x, p, t) #the harmonic oscillator ODE - dx[1] = x[2] - dx[2] = -x[1] - nothing - end + let + function harmosc!(dx, x, p, t) #the harmonic oscillator ODE + dx[1] = x[2] + dx[2] = -x[1] + nothing + end - @testset "Test TaylorN jet transport (t0,tmax): harmonic oscillator" begin - p = set_variables("ξ", numvars=2, order=5) - x0 = [-1.0,0.45] - x0TN = x0 + p - tvTN, xvTN = taylorinteg(harmosc!, x0TN, 0.0, 10pi, _order, _abstol, maxsteps=1) - @test length(tvTN) == 2 - @test size(xvTN) == (length(tvTN), length(x0TN)) - tvTN, xvTN = taylorinteg(harmosc!, x0TN, 0.0, 10pi, _order, _abstol) - tv , xv = taylorinteg(harmosc!, x0 , 0.0, 10pi, _order, _abstol) - x_analyticsol(t,x0,p0) = p0*sin(t)+x0*cos(t) - p_analyticsol(t,x0,p0) = p0*cos(t)-x0*sin(t) - x_δsol = x_analyticsol(tvTN[end], x0TN[1], x0TN[2])-xvTN[end,1] - x_δcoeffs = map(y->y[1], map(x->x.coeffs, x_δsol.coeffs)) - p_δsol = p_analyticsol(tvTN[end], x0TN[1], x0TN[2])-xvTN[end,2] - p_δcoeffs = map(y->y[1], map(x->x.coeffs, p_δsol.coeffs)) + @testset "Test TaylorN jet transport (t0,tmax): harmonic oscillator" begin + p = set_variables("ξ", numvars=2, order=5) + x0 = [-1.0,0.45] + x0TN = x0 + p + tvTN, xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + harmosc!, x0TN, 0.0, 10pi, _order, _abstol, maxsteps=1)) + @test length(tvTN) == 2 + @test size(xvTN) == (length(tvTN), length(x0TN)) + tvTN, xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0TN, 0.0, 10pi, _order, _abstol)) + tv , xv = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, x0 , 0.0, 10pi, _order, _abstol)) + x_analyticsol(t,x0,p0) = p0*sin(t)+x0*cos(t) + p_analyticsol(t,x0,p0) = p0*cos(t)-x0*sin(t) + x_δsol = x_analyticsol(tvTN[end], x0TN[1], x0TN[2])-xvTN[end,1] + x_δcoeffs = map(y->y[1], map(x->x.coeffs, x_δsol.coeffs)) + p_δsol = p_analyticsol(tvTN[end], x0TN[1], x0TN[2])-xvTN[end,2] + p_δcoeffs = map(y->y[1], map(x->x.coeffs, p_δsol.coeffs)) - @test (length(tvTN), length(x0)) == size(xvTN) - @test isapprox(x_δcoeffs, zeros(6), atol=1e-10, rtol=0) - @test isapprox(x0, map( x->evaluate(x), xvTN[1,:])) - @test isapprox(xv[1,:], map( x->evaluate(x), xvTN[1,:])) # nominal solution must coincide with jet evaluated at ξ=(0,0) at initial time - @test isapprox(xv[end,:], map( x->evaluate(x), xvTN[end,:])) #nominal solution must coincide with jet evaluated at ξ=(0,0) at final time + @test (length(tvTN), length(x0)) == size(xvTN) + @test isapprox(x_δcoeffs, zeros(6), atol=1e-10, rtol=0) + @test isapprox(x0, map( x->evaluate(x), xvTN[1,:])) + @test isapprox(xv[1,:], map( x->evaluate(x), xvTN[1,:])) # nominal solution must coincide with jet evaluated at ξ=(0,0) at initial time + @test isapprox(xv[end,:], map( x->evaluate(x), xvTN[end,:])) #nominal solution must coincide with jet evaluated at ξ=(0,0) at final time - xvTN_0 = map( x->evaluate(x), xvTN ) # the jet evaluated at the nominal solution - @test isapprox(xv[end,:], xvTN_0[end,:]) # nominal solution must coincide with jet evaluated at ξ=(0,0) at final time - @test isapprox(xvTN_0[1,:], xvTN_0[end,:]) # end point must coincide with a full period + xvTN_0 = map( x->evaluate(x), xvTN ) # the jet evaluated at the nominal solution + @test isapprox(xv[end,:], xvTN_0[end,:]) # nominal solution must coincide with jet evaluated at ξ=(0,0) at final time + @test isapprox(xvTN_0[1,:], xvTN_0[end,:]) # end point must coincide with a full period + end end function pendulum!(dx, x, p, t) #the simple pendulum ODE @@ -370,14 +428,17 @@ using Test #the time range tr = t0:integstep:tmax; #xv is the solution vector representing the propagation of the initial condition q0 propagated until time T - xv = taylorinteg(pendulum!, q0, tr, _order, _abstol, maxsteps=100) + xv = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0, tr, _order, _abstol, maxsteps=100)) #xvTN is the solution vector representing the propagation of the initial condition q0 plus variations (ξ₁,ξ₂) propagated until time T #note that q0 is a Vector{Float64}, but q0TN is a Vector{TaylorN{Float64}} #but apart from that difference, we're calling `taylorinteg` essentially with the same parameters! #thus, jet transport is reduced to a beautiful application of Julia's multiple dispatch! - xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=1) + xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, q0TN, tr, _order, _abstol, maxsteps=1)) @test size(xvTN) == (5,2) - xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=100) + xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0TN, tr, _order, _abstol, maxsteps=100)) xvTN_0 = map( x->evaluate(x, [0.0, 0.0]), xvTN ) # the jet evaluated at the nominal solution @@ -385,7 +446,8 @@ using Test @test isapprox(xv, xvTN_0) #nominal solution must coincide with jet evaluated at ξ=(0,0) #testing another jet transport method: - tv, xv = taylorinteg(pendulum!, q0, t0, tmax, _order, _abstol, maxsteps=100) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0, t0, tmax, _order, _abstol, maxsteps=100)) @test isapprox(xv[1,:], xvTN_0[1,:]) #nominal solution must coincide with jet evaluated at ξ=(0,0) at initial time @test isapprox(xv[end,:], xvTN_0[end,:]) #nominal solution must coincide with jet evaluated at ξ=(0,0) at final time @@ -393,7 +455,8 @@ using Test disp = 0.0001 #works even with 0.001, but we're giving it some margin #compare the jet solution evaluated at various variation vectors ξ, wrt to full solution, at each time evaluation point - xv_disp = taylorinteg(pendulum!, q0+[disp,0.0], tr, _order, _abstol, maxsteps=1) + xv_disp = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, q0+[disp,0.0], tr, _order, _abstol, maxsteps=1)) @test size(xv_disp) == (5,2) for i in 1:10 # generate a random angle @@ -401,7 +464,8 @@ using Test # generate a random point on a circumference of radius disp ξ = disp*[cos(ϕ), sin(ϕ)] #propagate in time full solution with initial condition q0+ξ - xv_disp = taylorinteg(pendulum!, q0+ξ, tr, _order, _abstol, maxsteps=100) + xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0+ξ, tr, _order, _abstol, maxsteps=100)) #evaluate jet at q0+ξ xvTN_disp = map( x->evaluate(x, ξ), xvTN ) #the propagated exact solution at q0+ξ should be approximately equal to the propagated jet solution evaluated at the same point: diff --git a/test/lyapunov.jl b/test/lyapunov.jl index e08798fb7..8c8c534f0 100644 --- a/test/lyapunov.jl +++ b/test/lyapunov.jl @@ -3,9 +3,13 @@ using TaylorIntegration using LinearAlgebra: norm, tr, dot, istriu, diag, I using Test +using Logging +import Logging: Warn @testset "Testing `lyapunov.jl`" begin + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + local _order = 28 local _abstol = 1.0E-20 @@ -114,13 +118,15 @@ using Test xi = set_variables("δ", order=1, numvars=length(q0)) #Test lyap_taylorinteg with autodiff-computed Jacobian, maxsteps=5 - tv, xv, λv = lyap_taylorinteg(lorenz!, q0, t0, tmax, _order, _abstol, nothing; maxsteps=5) + tv, xv, λv = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, t0, tmax, _order, _abstol, nothing; maxsteps=5)) @test size(tv) == (6,) @test size(xv) == (6,3) @test size(λv) == (6,3) #Test lyap_taylorinteg with user-defined Jacobian, maxsteps=5 - tv2, xv2, λv2 = lyap_taylorinteg(lorenz!, q0, t0, tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=5) + tv2, xv2, λv2 = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, t0, tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=5)) @test size(tv2) == (6,) @test size(xv2) == (6,3) @test size(λv2) == (6,3) @@ -129,7 +135,8 @@ using Test @test λv == λv2 #Test lyap_taylorinteg with autodiff-computed Jacobian, maxsteps=2000 - tv, xv, λv = lyap_taylorinteg(lorenz!, q0, t0, tmax, _order, _abstol; maxsteps=2000) + tv, xv, λv = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, t0, tmax, _order, _abstol; maxsteps=2000)) @test xv[1,:] == q0 @test tv[1] == t0 @test size(xv) == size(λv) @@ -144,12 +151,14 @@ using Test @test t_ == tv @test x_ == xv # test backward integration - tvb, xvb, λvb = lyap_taylorinteg(lorenz!, q0, t0, -tmax, _order, _abstol; maxsteps=2000) + tvb, xvb, λvb = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, t0, -tmax, _order, _abstol; maxsteps=2000)) @test isapprox(sum(λvb[1,:]), lorenztr) == false @test isapprox(sum(λvb[end,:]), lorenztr) #Test lyap_taylorinteg with user-defined Jacobian, maxsteps=2000 - tv2, xv2, λv2 = lyap_taylorinteg(lorenz!, q0, t0, tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000) + tv2, xv2, λv2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, t0, tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) @test xv2[1,:] == q0 @test tv2[1] == t0 @test size(xv2) == size(λv2) @@ -166,7 +175,8 @@ using Test @test t_ == tv2 @test x_ == xv2 # test backward integration - tv2b, xv2b, λv2b = lyap_taylorinteg(lorenz!, q0, t0, -tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000) + tv2b, xv2b, λv2b = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, t0, -tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) @test isapprox(sum(λv2b[1,:]), lorenztr) == false @test isapprox(sum(λv2b[end,:]), lorenztr) end @@ -185,28 +195,32 @@ using Test @test lorenztr == -(1+σ+β) trange = t0:1.0:tmax - xw, λw = lyap_taylorinteg(lorenz!, q0, trange, _order, _abstol; maxsteps=5) + xw, λw = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, trange, _order, _abstol; maxsteps=5)) @test xw[1,:] == q0 @test size(xw) == (length(trange), 3) @test λw[1,:] == zeros(3) @test size(λw) == size(xw) @test all(isnan.(xw[2:end,:])) @test all(isnan.(λw[2:end,:])) - xw2, λw2 = lyap_taylorinteg(lorenz!, q0, vec(trange), _order, _abstol; maxsteps=5) + xw2, λw2 = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, vec(trange), _order, _abstol; maxsteps=5)) @test xw[1, :] == xw2[1, :] @test all(isnan.(xw2[2:end,:])) @test λw2[1, :] == zeros(3) @test all(isnan.(λw2[2:end,:])) @test size(xw2) == (length(trange), 3) @test size(λw2) == size(xw2) - xw_, λw_ = lyap_taylorinteg(lorenz!, q0, trange, _order, _abstol, nothing, lorenz_jac!; maxsteps=5) + xw_, λw_ = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, trange, _order, _abstol, nothing, lorenz_jac!; maxsteps=5)) @test xw_[1,:] == q0 @test all(isnan.(xw_[2:end,:])) @test size(xw_) == (length(trange), 3) @test size(λw_) == size(xw_) @test λw_[1, :] == zeros(3) @test all(isnan.(λw_[2:end,:])) - xw2_, λw2_ = lyap_taylorinteg(lorenz!, q0, vec(trange), _order, _abstol, nothing, lorenz_jac!; maxsteps=5) + xw2_, λw2_ = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + lorenz!, q0, vec(trange), _order, _abstol, nothing, lorenz_jac!; maxsteps=5)) @test xw_[1, :] == xw2_[1, :] @test all(isnan.(xw2_[2:end,:])) @test λw2_[1, :] == zeros(3) @@ -214,13 +228,15 @@ using Test @test size(xw2_) == (length(trange), 3) @test size(λw2_) == size(xw2_) - xw, λw = lyap_taylorinteg(lorenz!, q0, trange, _order, _abstol; maxsteps=2000) + xw, λw = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, trange, _order, _abstol; maxsteps=2000)) @test xw[1,:] == q0 @test size(xw) == (length(trange), length(q0)) @test size(λw) == (length(trange), length(q0)) @test isapprox(sum(λw[1,:]), lorenztr) == false @test isapprox(sum(λw[end,:]), lorenztr) - tz, xz, λz = lyap_taylorinteg(lorenz!, q0, trange[1], trange[end], _order, _abstol; maxsteps=2000) + tz, xz, λz = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, trange[1], trange[end], _order, _abstol; maxsteps=2000)) @test λw[end,:] == λz[end,:] @test xw[end,:] == xz[end,:] @test tz[end] == trange[end] @@ -228,22 +244,26 @@ using Test @test isapprox(λw[end,1], 1.47167, rtol=mytol, atol=mytol) @test isapprox(λw[end,2], -0.00831, rtol=mytol, atol=mytol) @test isapprox(λw[end,3], -22.46336, rtol=mytol, atol=mytol) - xw_, λw_ = lyap_taylorinteg(lorenz!, q0, trange, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000) + xw_, λw_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, trange, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) @test xw == xw_ @test λw == λw_ - tz_, xz_, λz_ = lyap_taylorinteg(lorenz!, q0, trange[1], trange[end], _order, _abstol, nothing, lorenz_jac!; maxsteps=2000) + tz_, xz_, λz_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, trange[1], trange[end], _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) @test λw_[end,:] == λz_[end,:] @test xw_[end,:] == xz_[end,:] @test tz_[end] == trange[end] - xw2, λw2 = lyap_taylorinteg(lorenz!, q0, vec(trange), _order, _abstol; maxsteps=2000) + xw2, λw2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, vec(trange), _order, _abstol; maxsteps=2000)) @test xw2 == xw @test λw2 == λw @test xw2[1,:] == q0 @test size(xw) == (length(trange), length(q0)) @test isapprox(sum(λw2[1,:]), lorenztr) == false @test isapprox(sum(λw2[end,:]), lorenztr) - tz2, xz2, λz2 = lyap_taylorinteg(lorenz!, q0, trange[1], trange[end], _order, _abstol; maxsteps=2000) + tz2, xz2, λz2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, trange[1], trange[end], _order, _abstol; maxsteps=2000)) @test λw2[end,:] == λz2[end,:] @test xw2[end,:] == xz2[end,:] @test tz2[end] == trange[end] @@ -251,10 +271,12 @@ using Test @test isapprox(λw2[end,1], 1.47167, rtol=mytol, atol=mytol) @test isapprox(λw2[end,2], -0.00831, rtol=mytol, atol=mytol) @test isapprox(λw2[end,3], -22.46336, rtol=mytol, atol=mytol) - xw2_, λw2_ = lyap_taylorinteg(lorenz!, q0, vec(trange), _order, _abstol, nothing, lorenz_jac!; maxsteps=2000) + xw2_, λw2_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, vec(trange), _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) @test xw2 == xw2_ @test λw2 == λw2_ - tz2_, xz2_, λz2_ = lyap_taylorinteg(lorenz!, q0, trange[1], trange[end], _order, _abstol, nothing, lorenz_jac!; maxsteps=2000) + tz2_, xz2_, λz2_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz!, q0, trange[1], trange[end], _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) @test λw2_[end,:] == λz2_[end,:] @test xw2_[end,:] == xz2_[end,:] @test tz2_[end] == trange[end] diff --git a/test/many_ode.jl b/test/many_ode.jl index e56a9eb07..feb664603 100644 --- a/test/many_ode.jl +++ b/test/many_ode.jl @@ -3,6 +3,8 @@ using TaylorIntegration using Test using LinearAlgebra: Transpose, norm +using Logging +import Logging: Warn @testset "Testing `many_ode.jl`" begin @@ -10,6 +12,8 @@ using LinearAlgebra: Transpose, norm local _abstol = 1.0E-20 local tT = Taylor1(_order) + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + @testset "Tests: dot{x}=x.^2, x(0) = [3, 1]" begin function eqs_mov!(Dx, x, p, t) for i in eachindex(x) @@ -30,14 +34,16 @@ using LinearAlgebra: Transpose, norm δt = (_abstol/q0T[1].coeffs[end-1])^inv(_order-1) @test TaylorIntegration.stepsize(q0T, _abstol) == δt - tv, xv = taylorinteg(eqs_mov!, q0, 0.0, 0.5, _order, _abstol, nothing) + tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov!, q0, 0.0, 0.5, _order, _abstol, nothing)) @test length(tv) == 501 @test isa(xv, SubArray) @test xv[1,:] == q0 @test tv[end] < 1/3 trange = 0.0:1/8:1.0 - xv = taylorinteg(eqs_mov!, q0, trange, _order, _abstol) + xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov!, q0, trange, _order, _abstol)) @test size(xv) == (9,2) @test q0 == [3.0, 1.0] @test typeof(xv) == Transpose{Float64, Array{Float64,2}} @@ -48,7 +54,8 @@ using LinearAlgebra: Transpose, norm @test abs(xv[2,1] - 4.8) ≤ eps(4.8) tarray = vec(trange) - xv2 = taylorinteg(eqs_mov!, q0, tarray, _order, _abstol, nothing) + xv2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov!, q0, tarray, _order, _abstol, nothing)) @test xv[1:3,:] == xv2[1:3,:] @test xv2[1:3,:] ≈ xv[1:3,:] atol=eps() rtol=0.0 @test size(xv2) == (9,2) @@ -61,7 +68,8 @@ using LinearAlgebra: Transpose, norm @test abs(xv2[2,1] - 4.8) ≤ eps(4.8) # Output includes Taylor polynomial solution - tv, xv, polynV = taylorinteg(eqs_mov!, q0, 0, 0.5, _order, _abstol, Val(true), nothing, maxsteps=2) + tv, xv, polynV = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov!, q0, 0, 0.5, _order, _abstol, Val(true), nothing, maxsteps=2)) @test size(polynV) == (3, 2) @test xv[1,:] == q0 @test polynV[1,:] == Taylor1.(q0, _order) @@ -80,7 +88,8 @@ using LinearAlgebra: Transpose, norm q0 = [3.0, 3.0] tmax = 0.3 - tv, xv = taylorinteg(eqs_mov!, q0, 0, tmax, _order, _abstol, nothing) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + eqs_mov!, q0, 0, tmax, _order, _abstol, nothing)) @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 @@ -95,7 +104,8 @@ using LinearAlgebra: Transpose, norm @test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 5e-14 tmax = 0.33 - tv, xv = taylorinteg(eqs_mov!, [3, 3], 0.0, tmax, _order, _abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + eqs_mov!, [3, 3], 0.0, tmax, _order, _abstol)) @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 @@ -129,7 +139,8 @@ using LinearAlgebra: Transpose, norm order = 25 x0 = [t0, 0.0] #initial conditions such that x(t)=sin(t) - tv, xv = taylorinteg(f!, x0, t0, tmax, order, abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + f!, x0, t0, tmax, order, abstol)) @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 @@ -139,7 +150,8 @@ using LinearAlgebra: Transpose, norm @test abs(sin(tmax)-xv[end,2]) < 1e-14 # Backward integration - tb, xb = taylorinteg(f!, [tmax, sin(tmax)], tmax, t0, order, abstol) + tb, xb = (@test_logs min_level=Logging.Warn taylorinteg( + f!, [tmax, sin(tmax)], tmax, t0, order, abstol)) @test length(tb) < 501 @test length(xb[:,1]) < 501 @test length(xb[:,2]) < 501 @@ -152,13 +164,15 @@ using LinearAlgebra: Transpose, norm tmax = 15*(2pi) Δt = (tmax-t0)/1024 tspan = t0:Δt:tmax - xv = taylorinteg(f!, x0, tspan, order, abstol, nothing) + xv = (@test_logs min_level=Logging.Warn taylorinteg( + f!, x0, tspan, order, abstol, nothing)) @test xv[1,1:end] == x0 @test tmax == xv[end,1] @test abs(sin(tmax)-xv[end,2]) < 1e-14 # Backward integration - xback = taylorinteg(f!, xv[end, :], reverse(tspan), order, abstol, nothing) + xback = (@test_logs min_level=Logging.Warn taylorinteg( + f!, xv[end, :], reverse(tspan), order, abstol, nothing)) @test xback[1,:] == xv[end, :] @test abs(xback[end,1]-x0[1]) < 5.0e-14 @test abs(xback[end,2]-x0[2]) < 5.0e-14 @@ -181,7 +195,8 @@ using LinearAlgebra: Transpose, norm abstol = 1e-20 order = 10 x0 = [10.0, 0.0] #initial conditions such that x(t)=sin(t) - tv, xv = taylorinteg(fallball!, x0, t0, tmax, order, abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + fallball!, x0, t0, tmax, order, abstol)) @test length(tv) < 501 @test length(xv[:,1]) < 501 @test exactsol.(tv, x0[1], x0[2]) ≈ xv[:,1] diff --git a/test/one_ode.jl b/test/one_ode.jl index c6692d878..042b4c014 100644 --- a/test/one_ode.jl +++ b/test/one_ode.jl @@ -3,6 +3,8 @@ using TaylorIntegration using Test using LinearAlgebra: norm +using Logging +import Logging: Warn @testset "Testing `one_ode.jl`" begin @@ -10,6 +12,8 @@ using LinearAlgebra: norm local _abstol = 1.0E-20 local tT = Taylor1(_order) + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + @testset "Tests: dot{x}=x^2, x(0) = 1" begin eqs_mov(x, p, t) = x^2 t0 = 0.0 @@ -21,46 +25,53 @@ using LinearAlgebra: norm δt = _abstol^inv(_order-1) @test TaylorIntegration.stepsize(x0T, _abstol) == δt - tv, xv = taylorinteg(eqs_mov, 1, 0.0, 1.0, _order, _abstol) + tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov, 1, 0.0, 1.0, _order, _abstol)) @test length(tv) == 501 @test length(xv) == 501 @test xv[1] == x0 @test tv[end] < 1.0 - tv, xv = taylorinteg(eqs_mov, x0, 0.0, 1.0, _order, _abstol, nothing) + tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov, x0, 0.0, 1.0, _order, _abstol, nothing)) @test length(tv) == 501 @test length(xv) == 501 @test xv[1] == x0 @test tv[end] < 1.0 trange = 0.0:1/8:1.0 - xv = taylorinteg(eqs_mov, 1, trange, _order, _abstol) + xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov, 1, trange, _order, _abstol)) @test length(xv) == length(trange) @test typeof(xv) == Array{typeof(x0),1} @test xv[1] == x0 @test isnan(xv[end]) @test abs(xv[5] - 2.0) ≤ eps(2.0) - tvr, xvr = taylorinteg(eqs_mov, x0, trange[1], trange[end-1], _order, _abstol) + tvr, xvr = (@test_logs min_level=Logging.Warn taylorinteg( + eqs_mov, x0, trange[1], trange[end-1], _order, _abstol)) @test tvr[1] == trange[1] @test tvr[end] == trange[end-1] @test xvr[1] == xv[1] @test xvr[end] == xv[end-1] trange = 0.0:1/8:1.0 - xv = taylorinteg(eqs_mov, x0, trange, _order, _abstol, nothing) + xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov, x0, trange, _order, _abstol, nothing)) @test length(xv) == length(trange) @test typeof(xv) == Array{typeof(x0),1} @test xv[1] == x0 @test isnan(xv[end]) @test abs(xv[5] - 2.0) ≤ eps(2.0) - tvr, xvr = taylorinteg(eqs_mov, x0, trange[1], trange[end-1], _order, _abstol, nothing) + tvr, xvr = (@test_logs min_level=Logging.Warn taylorinteg( + eqs_mov, x0, trange[1], trange[end-1], _order, _abstol, nothing)) @test tvr[1] == trange[1] @test tvr[end] == trange[end-1] @test xvr[1] == xv[1] @test xvr[end] == xv[end-1] tarray = vec(trange) - xv2 = taylorinteg(eqs_mov, x0, tarray, _order, _abstol) + xv2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov, x0, tarray, _order, _abstol)) @test xv[1:end-1] == xv2[1:end-1] @test xv2[1:end-1] ≈ xv[1:end-1] atol=eps() rtol=0.0 @test length(xv2) == length(tarray) @@ -70,7 +81,8 @@ using LinearAlgebra: norm @test abs(xv2[5] - 2.0) ≤ eps(2.0) # Output includes Taylor polynomial solution - tv, xv, polynV = taylorinteg(eqs_mov, x0, 0, 0.5, _order, _abstol, Val(true), maxsteps=2) + tv, xv, polynV = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov, x0, 0, 0.5, _order, _abstol, Val(true), maxsteps=2)) @test length(polynV) == 3 @test xv[1] == x0 @test polynV[1] == Taylor1(x0, _order) @@ -91,7 +103,8 @@ using LinearAlgebra: norm δt = (_abstol/x0T.coeffs[end-1])^inv(_order-1) @test TaylorIntegration.stepsize(x0T, _abstol) == δt - tv, xv = taylorinteg(eqs_mov, x0, 0, tmax, _order, _abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + eqs_mov, x0, 0, tmax, _order, _abstol)) @test length(tv) < 501 @test length(xv) < 501 @test length(tv) == 14 @@ -102,7 +115,8 @@ using LinearAlgebra: norm @test abs(xv[end]-exactsol(tv[end], xv[1])) < 5e-14 tmax = 0.33 - tv, xv = taylorinteg(eqs_mov, x0, t0, tmax, _order, _abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + eqs_mov, x0, t0, tmax, _order, _abstol)) @test length(tv) < 501 @test length(xv) < 501 @test length(tv) == 28 @@ -121,7 +135,8 @@ using LinearAlgebra: norm order = 25 x0 = 0.0 #initial conditions such that x(t)=sin(t) - tv, xv = taylorinteg(f, x0, t0, tmax, order, abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0, t0, tmax, order, abstol)) @test length(tv) < 501 @test length(xv) < 501 @test xv[1] == x0 @@ -129,7 +144,8 @@ using LinearAlgebra: norm @test abs(sin(tmax)-xv[end]) < 1e-14 # Backward integration - tb, xb = taylorinteg(f, sin(tmax), tmax, t0, order, abstol) + tb, xb = (@test_logs min_level=Logging.Warn taylorinteg( + f, sin(tmax), tmax, t0, order, abstol)) @test length(tb) < 501 @test length(xb) < 501 @test xb[1] == sin(tmax) @@ -140,12 +156,14 @@ using LinearAlgebra: norm tmax = 15*(2pi) Δt = (tmax-t0)/1024 tspan = t0:Δt:tmax - xv = taylorinteg(f, x0, tspan, order, abstol) + xv = (@test_logs min_level=Logging.Warn taylorinteg( + f, x0, tspan, order, abstol)) @test xv[1] == x0 @test abs(sin(tmax)-xv[end]) < 1e-14 # Backward integration - xback = taylorinteg(f, xv[end], reverse(tspan), order, abstol) + xback = (@test_logs min_level=Logging.Warn taylorinteg( + f, xv[end], reverse(tspan), order, abstol)) @test xback[1] == xv[end] @test abs(sin(t0)-xback[end]) < 5e-14 @test norm(xv[:]-xback[end:-1:1], Inf) < 5.0e-14 diff --git a/test/rootfinding.jl b/test/rootfinding.jl index 937bd5850..cd712fb60 100644 --- a/test/rootfinding.jl +++ b/test/rootfinding.jl @@ -3,9 +3,14 @@ using TaylorIntegration using Test using LinearAlgebra: norm +using Logging +import Logging: Warn @testset "Testing `rootfinding.jl`" begin + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + + err_newton_raphson() = "Newton-Raphson did not converge for prescribed tolerance and maximum allowed iterations.\n" local _order = 28 local _abstol = 1.0E-20 @@ -27,8 +32,8 @@ using LinearAlgebra: norm x01 = x0 + [ξ, ξ] #warm-up lap and preliminary tests - tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, t0, Tend, - _order, _abstol, maxsteps=1) + tv, xv, tvS, xvS, gvS = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, g, x0, t0, Tend, _order, _abstol, maxsteps=1)) @test size(tv) == (2,) @test tv[1] == t0 @test size(xv) == (2,2) @@ -37,8 +42,8 @@ using LinearAlgebra: norm @test size(xvS) == (0,2) @test size(gvS) == (0,) - tvN, xvN, tvSN, xvSN, gvSN = taylorinteg(pendulum!, g, x0N, t0, 3Tend, - _order, _abstol, maxsteps=1) + tvN, xvN, tvSN, xvSN, gvSN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, g, x0N, t0, 3Tend, _order, _abstol, maxsteps=1)) @test eltype(tvN) == Float64 @test eltype(xvN) == TaylorN{Float64} @test eltype(tvSN) == TaylorN{Float64} @@ -52,8 +57,8 @@ using LinearAlgebra: norm @test size(xvSN) == (0,2) @test size(gvSN) == (0,) - tv1, xv1, tvS1, xvS1, gvS1 = taylorinteg(pendulum!, g, x01, t0, 3Tend, - _order, _abstol, maxsteps=1) + tv1, xv1, tvS1, xvS1, gvS1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, g, x01, t0, 3Tend, _order, _abstol, maxsteps=1)) @test eltype(tv1) == Float64 @test eltype(xv1) == Taylor1{Float64} @test eltype(tvS1) == Taylor1{Float64} @@ -68,8 +73,8 @@ using LinearAlgebra: norm @test size(gvS1) == (0,) #testing 0-th order root-finding - tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, t0, 3Tend, - _order, _abstol, maxsteps=1000) + tv, xv, tvS, xvS, gvS = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000)) @test tv[1] == t0 @test xv[1,:] == x0 @test size(tvS) == (5,) @@ -80,8 +85,8 @@ using LinearAlgebra: norm tvr = [t0, Tend/2, Tend, 3Tend/2, 2Tend, 5Tend/2, 3Tend] @test_throws AssertionError taylorinteg(pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000, eventorder=_order+1) - xvr, tvSr, xvSr, gvSr = taylorinteg(pendulum!, g, x0, view(tvr, :), - _order, _abstol, maxsteps=1000) + xvr, tvSr, xvSr, gvSr = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000)) @test xvr[1,:] == x0 @test size(tvSr) == (5,) @test size(tvSr) == size(tvr[2:end-1]) @@ -92,8 +97,8 @@ using LinearAlgebra: norm @test norm(tvS-tvSr, Inf) < 5E-15 #testing 0-th order root-finding + TaylorN jet transport - tvN, xvN, tvSN, xvSN, gvSN = taylorinteg(pendulum!, g, x0N, t0, 3Tend, - _order, _abstol, maxsteps=1000) + tvN, xvN, tvSN, xvSN, gvSN = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0N, t0, 3Tend, _order, _abstol, maxsteps=1000)) @test size(tvSN) == size(tvS) @test size(xvSN) == size(xvS) @test size(gvSN) == size(gvS) @@ -102,8 +107,8 @@ using LinearAlgebra: norm @test norm( xvSN()-xvS, Inf ) < 1E-14 #testing 0-th root-finding + Taylor1 jet transport - tv1, xv1, tvS1, xvS1, gvS1 = taylorinteg(pendulum!, g, x01, t0, 3Tend, - _order, _abstol, maxsteps=1000) + tv1, xv1, tvS1, xvS1, gvS1 = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x01, t0, 3Tend, _order, _abstol, maxsteps=1000)) @test size(tvS1) == size(tvS) @test size(xvS1) == size(xvS) @test size(gvS1) == size(gvS) @@ -115,8 +120,8 @@ using LinearAlgebra: norm @test_throws AssertionError taylorinteg(pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=_order+1, newtoniter=2) - tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, t0, 3Tend, - _order, _abstol, maxsteps=1000, eventorder=2, newtoniter=2) + tv, xv, tvS, xvS, gvS = (@test_logs (Warn, err_newton_raphson()) match_mode=:any taylorinteg( + pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=2, newtoniter=2)) @test tv[1] < tv[end] @test tv[1] == t0 @test xv[1,:] == x0 @@ -125,8 +130,8 @@ using LinearAlgebra: norm @test norm(gvS[:]) < 1E-15 # testing backward integrations - tvb, xvb, tvSb, xvSb, gvSb = taylorinteg(pendulum!, g, xv[end,:], 3Tend, t0, - _order, _abstol, maxsteps=1000, eventorder=2, newtoniter=2) + tvb, xvb, tvSb, xvSb, gvSb = (@test_logs (Warn, err_newton_raphson()) match_mode=:any taylorinteg( + pendulum!, g, xv[end,:], 3Tend, t0, _order, _abstol, maxsteps=1000, eventorder=2, newtoniter=2)) @test tvb[1] > tvb[end] @test tvSb[1] > tvSb[end] @test norm(gvSb[:]) < 1E-14 @@ -134,8 +139,8 @@ using LinearAlgebra: norm @test norm( xvSb[2:end,:] .- xvS[end:-1:1,:], Inf ) < 1E-13 #testing higher order root-finding + TaylorN jet transport - tvN, xvN, tvSN, xvSN, gvSN = taylorinteg(pendulum!, g, x0N, t0, 3Tend, - _order, _abstol, maxsteps=1000, eventorder=2) + tvN, xvN, tvSN, xvSN, gvSN = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0N, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=2)) @test size(tvSN) == size(tvS) @test size(xvSN) == size(xvS) @test size(gvSN) == size(gvS) @@ -144,8 +149,8 @@ using LinearAlgebra: norm @test norm(xvSN()-xvS, Inf) < 1E-14 #testing higher root-finding + Taylor1 jet transport - tv1, xv1, tvS1, xvS1, gvS1 = taylorinteg(pendulum!, g, x01, t0, 3Tend, - _order, _abstol, maxsteps=1000, eventorder=2) + tv1, xv1, tvS1, xvS1, gvS1 = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x01, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=2)) @test size(tvS1) == size(tvS) @test size(xvS1) == size(xvS) @test size(gvS1) == size(gvS) @@ -156,10 +161,10 @@ using LinearAlgebra: norm # Tests if trange is properly sorted Δt = (3Tend-t0)/1000 tspan = t0:Δt:(3Tend-0.125) - xv1r, tvS1r, xvS1r, gvS1r = taylorinteg(pendulum!, g, x01, tspan, - _order, _abstol, maxsteps=1000, eventorder=2) - xv1rb, tvS1rb, xvS1rb, gvS1rb = taylorinteg(pendulum!, g, xv1r[end,:], reverse(tspan), - _order, _abstol, maxsteps=1000, eventorder=2) + xv1r, tvS1r, xvS1r, gvS1r = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x01, tspan, _order, _abstol, maxsteps=1000, eventorder=2)) + xv1rb, tvS1rb, xvS1rb, gvS1rb = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, xv1r[end,:], reverse(tspan), _order, _abstol, maxsteps=1000, eventorder=2)) @test size(xv1r) == size(xv1rb) @test size(tvS1r) == size(tvS1rb) @test size(xvS1r) == size(xvS1rb) diff --git a/test/taylorize.jl b/test/taylorize.jl index 3364283f1..242d0ec6e 100644 --- a/test/taylorize.jl +++ b/test/taylorize.jl @@ -7,37 +7,47 @@ using InteractiveUtils: methodswith using Elliptic using Base.Threads import Pkg +using Logging +import Logging: Warn @testset "Testing `taylorize.jl`" begin # Constants for the integrations + local TI = TaylorIntegration local _order = 20 local _abstol = 1.0e-20 local t0 = 0.0 local tf = 1000.0 + unable_to_parse(f) = """Unable to use the parsed method of `jetcoeffs!` for `$f`, + despite of having `parse_eqs=true`, due to some internal error. + Using `parse_eqs = false`.""" + + max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + # Scalar integration @testset "Scalar case: xdot(x, p, t) = b-x^2" begin + # Use a global constant b1 = 3.0 @taylorize xdot1(x, p, t) = b1-x^2 @test (@isdefined xdot1) x0 = 1.0 - tv1, xv1 = taylorinteg( xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000, - parse_eqs=false) - tv1p, xv1p = taylorinteg( xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000) + tv1, xv1 = taylorinteg( xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000, parse_eqs=false) + tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000)) @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @test iszero( norm(xv1-xv1p, Inf) ) - # Now using `local` constants + # Use `local` constants @taylorize xdot2(x, p, t) = (local b2 = 3; b2-x^2) @test (@isdefined xdot2) - tv2, xv2 = taylorinteg( xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000, - parse_eqs=false) - tv2p, xv2p = taylorinteg( xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000) + tv2, xv2 = taylorinteg( xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000, parse_eqs=false) + tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000)) @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -49,7 +59,8 @@ import Pkg tv3, xv3 = taylorinteg( xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000, parse_eqs=false) - tv3p, xv3p = taylorinteg( xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000) + tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000)) @test length(tv3) == length(tv3p) @test iszero( norm(tv3-tv3p, Inf) ) @@ -64,7 +75,8 @@ import Pkg xv4 = taylorinteg( xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000, parse_eqs=false) - xv4p = taylorinteg( xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000) + xv4p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000)) @test iszero( norm(xv4-xv4p, Inf) ) # Compare to exact solution @@ -72,37 +84,56 @@ import Pkg ((sqrt(b)+x0)+(sqrt(b)-x0)*exp(-2sqrt(b)*t)) @test norm(xv1p[end] - exact_sol(tv1p[end], b1, x0), Inf) < 1.0e-15 - # Check that the parsed `jetcoeffs` produces the correct series in `x` + # Check that the parsed `jetcoeffs` produces the correct series in `x` and no error + # TODO: Use metaprogramming here tT = t0 + Taylor1(_order) xT = x0 + zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), xdot1, tT, xT, nothing) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(xdot1), tT, xT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(xdot1), tT, xT, nothing, tmpTaylor, arrTaylor) @test xT ≈ exact_sol(tT, b1, x0) xT = x0 + zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), xdot2, tT, xT, nothing) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(xdot2), tT, xT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(xdot2), tT, xT, nothing, tmpTaylor, arrTaylor) @test xT ≈ exact_sol(tT, 3.0, x0) xT = x0 + zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), xdot2, tT, xT, b1) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(xdot2), tT, xT, b1)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(xdot2), tT, xT, b1, tmpTaylor, arrTaylor) @test xT ≈ exact_sol(tT, b1, x0) - # The macro returns a (parsed) jetcoeffs! function which yields an error and - # therefore it runs with the default jetcoeffs! method. A warning is issued. + # The macro returns a (parsed) jetcoeffs! function which yields a `MethodError` + # (TaylorSeries.identity! requires *two* Taylor series and an integer) and + # therefore it runs with the default jetcoeffs! method (i.e., parse_eqs=false). + # A warning is issued. @taylorize function xdot2_err(x, p, t) - b2 = 3 + b2 = 3 # This line is parsed as `identity!` and yields a MethodError; + # correct this by replacing it to `b2 = 3 + zero(t)` return b2-x^2 end @test (@isdefined xdot2_err) - @test !isempty(methodswith(Val{xdot2_err}, TaylorIntegration.jetcoeffs!)) - tv2e, xv2e = taylorinteg( xdot2_err, x0, t0, tf, _order, _abstol, maxsteps=1000) + @test !isempty(methodswith(Val{xdot2_err}, TI.jetcoeffs!)) + tv2e, xv2e = (@test_logs (Warn, unable_to_parse(xdot2_err)) taylorinteg( + xdot2_err, x0, t0, tf, _order, _abstol, maxsteps=1000)) @test length(tv2) == length(tv2e) @test iszero( norm(tv2-tv2e, Inf) ) @test iszero( norm(xv2-xv2e, Inf) ) - @test_throws UndefVarError TaylorIntegration.__jetcoeffs!(Val(true), xdot2_err, tT, qT, similar(qT), similar(qT), [1.0]) + xT = x0 + zero(tT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(xdot2_err), tT, xT, nothing)) + @test_throws MethodError TI.jetcoeffs!(Val(xdot2_err), tT, xT, nothing, + tmpTaylor, arrTaylor) # Output includes Taylor polynomial solution - tv3t, xv3t, polynV3t = taylorinteg(xdot3, x0, t0, tf, _order, _abstol, Val(true), 0.0, maxsteps=2) + tv3t, xv3t, polynV3t = (@test_logs (Warn, max_iters_reached()) taylorinteg( + xdot3, x0, t0, tf, _order, _abstol, Val(true), 0.0, maxsteps=2)) @test length(polynV3t) == 3 @test xv3t[1] == x0 @test polynV3t[1] == Taylor1(x0, _order) @@ -112,27 +143,30 @@ import Pkg @testset "Scalar case: xdot(x, p, t) = -10" begin - xdot1(x, p, t) = -10 + zero(t) # `zero(t)` is needed; cf #20 - @taylorize xdot1_parsed(x, p, t) = -10# `zero(t)` can be avoided here + # Return an expression + xdot1(x, p, t) = -10 + zero(t) # `zero(t)` is needed; cf #20 + @taylorize xdot1_parsed(x, p, t) = -10 # `zero(t)` can be avoided here ! @test (@isdefined xdot1_parsed) tv1, xv1 = taylorinteg( xdot1, 10, 1, 20.0, _order, _abstol) - tv1p, xv1p = taylorinteg( xdot1_parsed, 10.0, 1.0, 20.0, _order, _abstol) + tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot1_parsed, 10.0, 1.0, 20.0, _order, _abstol)) @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @test iszero( norm(xv1-xv1p, Inf) ) - # Now using `local` + # Use `local` expression @taylorize function xdot2(x, p, t) - local ggrav = -10 + zero(t) # zero(t) is needed for the parse_eqs=false case + local ggrav = -10 + zero(t) # zero(t) is needed for the `parse_eqs=false` case tmp = ggrav # needed to avoid an error when parsing end @test (@isdefined xdot2) tv2, xv2 = taylorinteg( xdot2, 10.0, 1.0, 20.0, _order, _abstol, parse_eqs=false) - tv2p, xv2p = taylorinteg( xdot2, 10.0, 1, 20.0, _order, _abstol) + tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot2, 10.0, 1, 20.0, _order, _abstol)) @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -145,7 +179,8 @@ import Pkg tv3, xv3 = taylorinteg( xdot3, 10.0, 1.0, 20.0, _order, _abstol, -10, parse_eqs=false) - tv3p, xv3p = taylorinteg( xdot3, 10, 1.0, 20.0, _order, _abstol, -10) + tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot3, 10, 1.0, 20.0, _order, _abstol, -10)) @test length(tv3) == length(tv3p) @test iszero( norm(tv3-tv3p, Inf) ) @@ -162,18 +197,27 @@ import Pkg exact_sol(t, g, x0) = x0 + g*(t-1.0) @test norm(xv1p[end] - exact_sol(20, -10, 10), Inf) < 10eps(exact_sol(20, -10, 10)) - # Check that the parsed `jetcoeffs` produces the correct series in `x` + # Check that the parsed `jetcoeffs` produces the correct series in `x` and no errors tT = 1.0 + Taylor1(_order) xT = 10.0 + zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), xdot1_parsed, tT, xT, nothing) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(xdot1_parsed), tT, xT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(xdot1_parsed), tT, xT, nothing, tmpTaylor, arrTaylor) @test xT ≈ exact_sol(tT, -10, 10) xT = 10.0 + zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), xdot2, tT, xT, nothing) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(xdot2), tT, xT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(xdot2), tT, xT, nothing, tmpTaylor, arrTaylor) @test xT ≈ exact_sol(tT, -10, 10) xT = 10.0 + zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), xdot3, tT, xT, -10) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(xdot3), tT, xT, -10)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(xdot3), tT, xT, -10, tmpTaylor, arrTaylor) @test xT ≈ exact_sol(tT, -10, 10) end @@ -190,7 +234,8 @@ import Pkg q0 = [pi-0.001, 0.0] tv2, xv2 = taylorinteg(pendulum!, q0, t0, tf, _order, _abstol, parse_eqs=false, maxsteps=5000) - tv2p, xv2p = taylorinteg(pendulum!, q0, t0, tf, _order, _abstol, maxsteps=5000) + tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0, t0, tf, _order, _abstol, maxsteps=5000)) @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -204,11 +249,14 @@ import Pkg @test sol1.t == sol2.t == sol3.t == tv2p @test sol1.u[end] == sol2.u[end] == sol3.u[end] == xv2p[end,1:2] - # Check that the parsed `jetcoeffs` produces the correct series in `x`; - # we check that it does not throws an error - tT = 0.0 + Taylor1(_order) + # Check that the parsed `jetcoeffs!` produces the correct series in `x` and no errors + tT = t0 + Taylor1(_order) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), pendulum!, tT, qT, similar(qT), similar(qT), nothing) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(pendulum!), tT, qT, dqT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(pendulum!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor) end @@ -221,7 +269,8 @@ import Pkg cx0 = complex(1.0, 0.0) tv1, xv1 = taylorinteg(eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500, parse_eqs=false) - tv1p, xv1p = taylorinteg(eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500) + tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( + eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500)) @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -233,7 +282,8 @@ import Pkg @test (@isdefined eqscmplx2) tv2, xv2 = taylorinteg(eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500, parse_eqs=false) - tv2p, xv2p = taylorinteg(eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500) + tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( + eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500)) @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -244,7 +294,8 @@ import Pkg @test (@isdefined eqscmplx3) tv3, xv3 = taylorinteg(eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500, parse_eqs=false) - tv3p, xv3p = taylorinteg(eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500) + tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( + eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500)) @test length(tv3) == length(tv3p) @test iszero( norm(tv3-tv3p, Inf) ) @@ -278,14 +329,21 @@ import Pkg zz0 = [z0, z0] ts = 0.0:pi:2pi zsol = taylorinteg(eqs3!, zz0, ts, _order, _abstol, parse_eqs=false, maxsteps=10) - zsolp = taylorinteg(eqs3!, zz0, ts, _order, _abstol, maxsteps=10) + zsolp = (@test_logs min_level=Logging.Warn taylorinteg( + eqs3!, zz0, ts, _order, _abstol, maxsteps=10)) @test length(tv3) == length(tv3p) @test iszero( norm(tv3-tv3p, Inf) ) @test iszero( norm(xv3-xv3p, Inf) ) tT = t0 + Taylor1(_order) zT = zz0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), eqs3!, tT, zT, similar(zT), similar(zT), nothing) + dzT = similar(zT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(eqs3!), tT, zT, dzT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!( + Val(eqs3!), tT, zT, dzT, nothing, tmpTaylor, arrTaylor) + @test typeof(tmpTaylor) == Vector{Taylor1{ComplexF64}} + @test typeof(arrTaylor) == Vector{Vector{Taylor1{ComplexF64}}} @test zT[1] == exact_sol(tT, -1, z0) @test zT[2] == exact_sol(tT, z0, z0) end @@ -306,13 +364,15 @@ import Pkg @test (@isdefined integ_cos2) tv11, xv11 = taylorinteg(integ_cos1, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false) - tv12, xv12 = taylorinteg(integ_cos1, 0.0, 0.0, pi, _order, _abstol) + tv12, xv12 = (@test_logs min_level=Logging.Warn taylorinteg( + integ_cos1, 0.0, 0.0, pi, _order, _abstol)) @test length(tv11) == length(tv12) @test iszero( norm(tv11-tv12, Inf) ) @test iszero( norm(xv11-xv12, Inf) ) tv21, xv21 = taylorinteg(integ_cos2, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false) - tv22, xv22 = taylorinteg(integ_cos2, 0.0, 0.0, pi, _order, _abstol) + tv22, xv22 = (@test_logs min_level=Logging.Warn taylorinteg( + integ_cos2, 0.0, 0.0, pi, _order, _abstol)) @test length(tv21) == length(tv22) @test iszero( norm(tv21-tv22, Inf) ) @test iszero( norm(xv21-xv22, Inf) ) @@ -335,7 +395,8 @@ import Pkg @test (@isdefined integ_vec) x0 = [0.0, 1.0] tv11, xv11 = taylorinteg(integ_vec, x0, 0.0, pi, _order, _abstol, parse_eqs=false) - tv12, xv12 = taylorinteg(integ_vec, x0, 0.0, pi, _order, _abstol) + tv12, xv12 = (@test_logs min_level=Logging.Warn taylorinteg( + integ_vec, x0, 0.0, pi, _order, _abstol)) @test length(tv11) == length(tv12) @test iszero( norm(tv11-tv12, Inf) ) @test iszero( norm(xv11-xv12, Inf) ) @@ -359,10 +420,11 @@ import Pkg @test (@isdefined harm_osc!) q0 = [1.0, 0.0] p = [2.0] + local tf = 300.0 tv1, xv1 = taylorinteg(harm_osc!, q0, t0, tf, _order, _abstol, p, maxsteps=1000, parse_eqs=false) - tv1p, xv1p = taylorinteg(harm_osc!, q0, t0, tf, _order, _abstol, - p, maxsteps=1000) + tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( + harm_osc!, q0, t0, tf, _order, _abstol, p, maxsteps=1000)) @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -372,15 +434,20 @@ import Pkg @test norm(xv1p[end, 1] - cos(p[1]*tv1p[end]), Inf) < 1.0e-12 @test norm(xv1p[end, 2] + p[1]*sin(p[1]*tv1p[end]), Inf) < 3.0e-12 - # Check that the parsed `jetcoeffs` produces the correct series in `x` + # Check that the parsed `jetcoeffs` produces the correct series in `x` and no errors tT = t0 + Taylor1(_order) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), harm_osc!, tT, qT, similar(qT), similar(qT), [1.0]) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(harm_osc!), tT, qT, dqT, [1.0])) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(harm_osc!), tT, qT, dqT, [1.0], tmpTaylor, arrTaylor) @test qT[1] ≈ cos(tT) @test qT[2] ≈ -sin(tT) - # The macro returns a (parsed) jetcoeffs! function which yields an error and - # therefore it runs with the default jetcoeffs! methdo. A warining is issued. + # The macro returns a (parsed) jetcoeffs! function which yields a `MethodError` + # (TaylorSeries.pow! requires *two* Taylor series as first arguments) and + # therefore it runs with the default jetcoeffs! methdo. A warning is issued. + # To solve it, define as a local variable `ω2 = ω^2`. @taylorize function harm_osc_error!(dx, x, p, t) local ω = p[1] dx[1] = x[2] @@ -388,16 +455,24 @@ import Pkg return nothing end - @test !isempty(methodswith(Val{harm_osc_error!}, TaylorIntegration.jetcoeffs!)) - tv2e, xv2e = taylorinteg(harm_osc_error!, q0, t0, tf, _order, _abstol, - p, maxsteps=1000, parse_eqs=true) + @test !isempty(methodswith(Val{harm_osc_error!}, TI.jetcoeffs!)) + tv2e, xv2e = ( + @test_logs (Warn, unable_to_parse(harm_osc_error!)) taylorinteg( + harm_osc_error!, q0, t0, tf, _order, _abstol, p, maxsteps=1000, parse_eqs=true)) @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv2e, Inf) ) @test iszero( norm(xv1-xv2e, Inf) ) - @test_throws MethodError TaylorIntegration.__jetcoeffs!(Val(true), harm_osc_error!, tT, qT, similar(qT), similar(qT), [1.0]) - end + tT = t0 + Taylor1(_order) + qT = q0 .+ zero(tT) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(harm_osc_error!), tT, qT, dqT, [1.0])) + @test_throws MethodError TI.jetcoeffs!(Val(harm_osc_error!), tT, qT, dqT, [1.0], + tmpTaylor, arrTaylor) + end + local tf = 100.0 @testset "Multiple pendula" begin NN = 3 nnrange = 1:3 @@ -413,8 +488,8 @@ import Pkg q0 = [pi-0.001, 0.0, pi-0.001, 0.0, pi-0.001, 0.0] tv1, xv1 = taylorinteg(multpendula1!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000, parse_eqs=false) - tv1p, xv1p = taylorinteg(multpendula1!, q0, t0, tf, _order, _abstol, - [NN, nnrange], maxsteps=1000) + tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( + multpendula1!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000)) @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -433,15 +508,16 @@ import Pkg @test (@isdefined multpendula2!) tv2, xv2 = taylorinteg(multpendula2!, q0, t0, tf, _order, _abstol, maxsteps=1000, parse_eqs=false) - tv2p, xv2p = taylorinteg(multpendula2!, q0, t0, tf, _order, _abstol, - maxsteps=1000) + tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( + multpendula2!, q0, t0, tf, _order, _abstol, maxsteps=1000)) @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @test iszero( norm(xv2-xv2p, Inf) ) @taylorize function multpendula3!(dx, x, p, t) - nn, nnrng = p # `local` is not needed to reassign `p`; internally is treated as local) + nn, nnrng = p # `local` is not needed to reassign `p`; + # internally is treated as local) for i in nnrng dx[i] = x[nn+i] dx[i+nn] = -sin( x[i] ) @@ -452,8 +528,8 @@ import Pkg @test (@isdefined multpendula3!) tv3, xv3 = taylorinteg(multpendula3!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000, parse_eqs=false) - tv3p, xv3p = taylorinteg(multpendula3!, q0, t0, tf, _order, _abstol, - [NN, nnrange], maxsteps=1000) + tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( + multpendula3!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000)) # Comparing integrations @test length(tv1) == length(tv2) == length(tv3) @@ -462,17 +538,23 @@ import Pkg @test iszero( norm(xv1-xv2, Inf) ) @test iszero( norm(xv1-xv3, Inf) ) - # Check that the parsed `jetcoeffs` produces the correct series in `x`; - # we check that it does not throws an error + # Check that the parsed `jetcoeffs` produces the correct series in `x` and no errors tT = 0.0 + Taylor1(_order) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), multpendula1!, tT, qT, similar(qT), similar(qT), (NN, nnrange)) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(multpendula1!), tT, qT, dqT, (NN, nnrange))) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(multpendula1!), tT, qT, dqT, (NN, nnrange), tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), multpendula2!, tT, qT, similar(qT), similar(qT), nothing) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(multpendula2!), tT, qT, dqT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(multpendula2!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), multpendula3!, tT, qT, similar(qT), similar(qT), [NN, nnrange]) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(multpendula3!), tT, qT, dqT, [NN, nnrange])) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(multpendula3!), tT, qT, dqT, [NN, nnrange], tmpTaylor, arrTaylor) end @@ -515,13 +597,13 @@ import Pkg q0 = [0.2, 0.0, 0.0, 3.0] tv1, xv1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000, parse_eqs=false) - tv1p, xv1p = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0, - maxsteps=500000) + tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0, + maxsteps=500000)) tv6p, xv6p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars, maxsteps=500000, parse_eqs=false) - tv7p, xv7p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars, - maxsteps=500000) + tv7p, xv7p = (@test_logs min_level=Logging.Warn taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars, + maxsteps=500000)) @test length(tv1) == length(tv1p) == length(tv6p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -564,13 +646,13 @@ import Pkg @test (@isdefined kepler4!) tv3, xv3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, maxsteps=500000, parse_eqs=false) - tv3p, xv3p = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, - maxsteps=500000) + tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, + maxsteps=500000)) tv4, xv4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, maxsteps=500000, parse_eqs=false) - tv4p, xv4p = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, - maxsteps=500000) + tv4p, xv4p = (@test_logs min_level=Logging.Warn taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, + maxsteps=500000)) @test length(tv3) == length(tv3p) == length(tv4) @test iszero( norm(tv3-tv3p, Inf) ) @@ -588,16 +670,28 @@ import Pkg # we check that it does not throws an error tT = 0.0 + Taylor1(_order) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler1!, tT, qT, similar(qT), similar(qT), -1.0) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler1!), + tT, qT, dqT, -1.0)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler1!), tT, qT, dqT, -1.0, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler2!, tT, qT, similar(qT), similar(qT), pars) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler2!), + tT, qT, dqT, pars)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler2!), tT, qT, dqT, pars, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler3!, tT, qT, similar(qT), similar(qT), nothing) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler3!), + tT, qT, dqT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler3!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler4!, tT, qT, similar(qT), similar(qT), nothing) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler4!), + tT, qT, dqT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler4!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor) end @@ -638,13 +732,13 @@ import Pkg q0 = [0.2, 0.0, 0.0, 3.0] tv1, xv1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, maxsteps=500000, parse_eqs=false) - tv1p, xv1p = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, - maxsteps=500000) + tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler1!, q0, t0, tf, _order, _abstol, maxsteps=500000)) tv2, xv2 = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000, parse_eqs=false) - tv2p, xv2p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, -1.0, - maxsteps=500000) + tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler2!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000)) @test length(tv1) == length(tv1p) == length(tv2) @test iszero( norm(tv1-tv1p, Inf) ) @@ -690,13 +784,13 @@ import Pkg @test (@isdefined kepler4!) tv3, xv3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000, parse_eqs=false) - tv3p, xv3p = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, -1.0, - maxsteps=500000) + tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler3!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000)) tv4, xv4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, maxsteps=500000, parse_eqs=false) - tv4p, xv4p = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, - maxsteps=500000) + tv4p, xv4p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler4!, q0, t0, tf, _order, _abstol, maxsteps=500000)) @test length(tv3) == length(tv3p) == length(tv4) @test iszero( norm(tv3-tv3p, Inf) ) @@ -714,16 +808,25 @@ import Pkg # we check that it does not throws an error tT = 0.0 + Taylor1(_order) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler1!, tT, qT, similar(qT), similar(qT), nothing) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler1!), + tT, qT, dqT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler1!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler2!, tT, qT, similar(qT), similar(qT), -1.0) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler2!), + tT, qT, dqT, -1.0)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler2!), tT, qT, dqT, -1.0, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler3!, tT, qT, similar(qT), similar(qT), -1.0) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler3!), + tT, qT, dqT, -1.0)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler3!), tT, qT, dqT, -1.0, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), kepler4!, tT, qT, similar(qT), similar(qT), nothing) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!(Val(kepler4!), + tT, qT, dqT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(kepler4!), tT, qT, dqT, nothing, tmpTaylor, arrTaylor) end @@ -762,20 +865,21 @@ import Pkg xi = set_variables("δ", order=1, numvars=length(q0)) tv1, lv1, xv1 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, - params, maxsteps=2000, parse_eqs=false); + params, maxsteps=2000, parse_eqs=false) - tv1p, lv1p, xv1p = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, - params, maxsteps=2000); + tv1p, lv1p, xv1p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz1!, q0, t0, tf, _order, _abstol, params, maxsteps=2000)) @test tv1 == tv1p @test lv1 == lv1p @test xv1 == xv1p tv2, lv2, xv2 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, - params, lorenz1_jac!, maxsteps=2000, parse_eqs=false); + params, lorenz1_jac!, maxsteps=2000, parse_eqs=false) - tv2p, lv2p, xv2p = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, - params, lorenz1_jac!, maxsteps=2000, parse_eqs=true); + tv2p, lv2p, xv2p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz1!, q0, t0, tf, _order, _abstol, params, lorenz1_jac!, maxsteps=2000, + parse_eqs=true)) @test tv2 == tv2p @test lv2 == lv2p @@ -819,20 +923,21 @@ import Pkg end tv3, lv3, xv3 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, - maxsteps=2000, parse_eqs=false); + maxsteps=2000, parse_eqs=false) - tv3p, lv3p, xv3p = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, - maxsteps=2000); + tv3p, lv3p, xv3p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz2!, q0, t0, tf, _order, _abstol, maxsteps=2000)) @test tv3 == tv3p @test lv3 == lv3p @test xv3 == xv3p tv4, lv4, xv4 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, - lorenz2_jac!, maxsteps=2000, parse_eqs=false); + lorenz2_jac!, maxsteps=2000, parse_eqs=false) - tv4p, lv4p, xv4p = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, - lorenz2_jac!, maxsteps=2000, parse_eqs=true); + tv4p, lv4p, xv4p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz2!, q0, t0, tf, _order, _abstol, lorenz2_jac!, maxsteps=2000, + parse_eqs=true)) @test tv4 == tv4p @test lv4 == lv4p @@ -850,29 +955,39 @@ import Pkg # Using ranges lv5, xv5 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, - maxsteps=2000, parse_eqs=false); + maxsteps=2000, parse_eqs=false) - lv5p, xv5p = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, - maxsteps=2000, parse_eqs=true); + lv5p, xv5p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz2!, q0, t0:0.125:tf, _order, _abstol, maxsteps=2000, parse_eqs=true)) @test lv5 == lv5p @test xv5 == xv5p lv6, xv6 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, - lorenz2_jac!, maxsteps=2000, parse_eqs=false); + lorenz2_jac!, maxsteps=2000, parse_eqs=false) - lv6p, xv6p = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, - lorenz2_jac!, maxsteps=2000, parse_eqs=true); + lv6p, xv6p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + lorenz2!, q0, t0:0.125:tf, _order, _abstol, lorenz2_jac!, maxsteps=2000, + parse_eqs=true)) @test lv6 == lv6p @test xv6 == xv6p + # Check that no errors are thrown tT = 0.0 + Taylor1(_order) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), lorenz1!, tT, qT, similar(qT), similar(qT), params) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(lorenz1!), tT, qT, dqT, params)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(lorenz1!), + tT, qT, dqT, params, tmpTaylor, arrTaylor) qT = q0 .+ zero(tT) - TaylorIntegration.__jetcoeffs!(Val(true), lorenz2!, tT, qT, similar(qT), similar(qT), nothing) + dqT = similar(qT) + tmpTaylor, arrTaylor = (@test_logs min_level=Logging.Warn TI._allocate_jetcoeffs!( + Val(lorenz2!), tT, qT, dqT, nothing)) + @test_logs min_level=Logging.Warn TI.jetcoeffs!(Val(lorenz2!), + tT, qT, dqT, nothing, tmpTaylor, arrTaylor) end @@ -882,39 +997,39 @@ import Pkg dx[1] = x[2] dx[2] = -sin( x[1] ) end) - @test_throws ArgumentError TaylorIntegration._make_parsed_jetcoeffs(ex) + @test_throws ArgumentError TI._make_parsed_jetcoeffs(ex) # `&&` is not yet implemented ex = :(function f_p!(x, p, t) true && x end) - @test_throws ArgumentError TaylorIntegration._make_parsed_jetcoeffs(ex) + @test_throws ArgumentError TI._make_parsed_jetcoeffs(ex) # a is not an Expr; String ex = :(function f_p!(x, p, t) "a" end) - @test_throws ArgumentError TaylorIntegration._make_parsed_jetcoeffs(ex) + @test_throws ArgumentError TI._make_parsed_jetcoeffs(ex) # KeyError: key :fname not found ex = :(begin x=1 x+x end) - @test_throws KeyError TaylorIntegration._make_parsed_jetcoeffs(ex) + @test_throws KeyError TI._make_parsed_jetcoeffs(ex) # BoundsError; no return variable defined ex = :(function f_p!(x, p, t) local cos(t) end) - @test_throws BoundsError TaylorIntegration._make_parsed_jetcoeffs(ex) + @test_throws BoundsError TI._make_parsed_jetcoeffs(ex) # ArgumentError; .= is not yet implemented ex = :(function err_bbroadcasting!(Dz, z, p, t) Dz .= z nothing end) - @test_throws ArgumentError TaylorIntegration._make_parsed_jetcoeffs(ex) + @test_throws ArgumentError TI._make_parsed_jetcoeffs(ex) # The macro works fine, but the `jetcoeffs!` method does not work # (so the integration would silently run using `parse_eqs=false`) @@ -927,7 +1042,8 @@ import Pkg end tT = t0 + Taylor1(_order) qT = [1.0, 0.0] .+ zero(tT) - @test_throws MethodError TaylorIntegration.__jetcoeffs!(Val(true), harm_osc!, tT, qT, similar(qT), similar(qT), [2.0]) + @test_throws MethodError TI.__jetcoeffs!( + Val(true), harm_osc!, tT, qT, similar(qT), similar(qT), [2.0]) @taylorize function kepler1!(dq, q, p, t) μ = p @@ -942,7 +1058,8 @@ import Pkg end tT = t0 + Taylor1(_order) qT = [0.2, 0.0, 0.0, 3.0] .+ zero(tT) - @test_throws MethodError TaylorIntegration.__jetcoeffs!(Val(true), kepler1!, tT, qT, similar(qT), similar(qT), -1.0) + @test_throws MethodError TI.__jetcoeffs!( + Val(true), kepler1!, tT, qT, similar(qT), similar(qT), -1.0) end @@ -964,22 +1081,28 @@ import Pkg #the time range tr = t0:integstep:T; #note that as called below, taylorinteg uses the parsed jetcoeffs! method by default - xvp = taylorinteg(pendulum!, q0, tr, _order, _abstol, maxsteps=100) + xvp = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0, tr, _order, _abstol, maxsteps=100)) # "warmup" for jet transport integration - xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=1, parse_eqs=false) - xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=1) + xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, q0TN, tr, _order, _abstol, maxsteps=1, parse_eqs=false)) + xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, q0TN, tr, _order, _abstol, maxsteps=1)) @test size(xvTN) == (5,2) #jet transport integration with parsed jetcoeffs! - xvTNp = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=100) + xvTNp = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0TN, tr, _order, _abstol, maxsteps=100)) #jet transport integration with non-parsed jetcoeffs! - xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=100, parse_eqs=false) + xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0TN, tr, _order, _abstol, maxsteps=100, parse_eqs=false)) @test xvTN == xvTNp @test norm(xvTNp[:,:]() - xvp, Inf) < 1e-15 dq = 0.0001rand(2) q1 = q0 + dq - yv = taylorinteg(pendulum!, q1, tr, _order, _abstol, maxsteps=100) + yv = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q1, tr, _order, _abstol, maxsteps=100)) yv_jt = xvTNp[:,:](dq) @test norm(yv-yv_jt, Inf) < 1e-11 @@ -987,9 +1110,12 @@ import Pkg t = Taylor1([0.0, 1.0], 10) x0T1 = q0+[0t,t] q1 = q0+[0.0,dq] - tv, xv = taylorinteg(pendulum!, q1, t0, 2T, _order, _abstol) - tvT1, xvT1 = taylorinteg(pendulum!, x0T1, t0, 2T, _order, _abstol, parse_eqs=false) - tvT1p, xvT1p = taylorinteg(pendulum!, x0T1, t0, 2T, _order, _abstol) + tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q1, t0, 2T, _order, _abstol)) + tvT1, xvT1 = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, x0T1, t0, 2T, _order, _abstol, parse_eqs=false)) + tvT1p, xvT1p = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, x0T1, t0, 2T, _order, _abstol)) @test tvT1 == tvT1p @test xvT1 == xvT1p xv_jt = xvT1p[:,:](dq) @@ -1024,13 +1150,14 @@ import Pkg x01 = x0 + [ξ, ξ] #warm-up lap and preliminary tests - taylorinteg(pendulum!, g, x0, t0, Tend, _order, _abstol, maxsteps=1) - @test_throws AssertionError taylorinteg(pendulum!, g, x0, t0, Tend, - _order, _abstol, maxsteps=1, eventorder=_order+1) + @test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, g, x0, t0, Tend, _order, _abstol, maxsteps=1) + @test_throws AssertionError taylorinteg( + pendulum!, g, x0, t0, Tend, _order, _abstol, maxsteps=1, eventorder=_order+1) #testing 0-th order root-finding - tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, t0, 3Tend, - _order, _abstol, maxsteps=1000) + tv, xv, tvS, xvS, gvS = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000)) @test tv[1] == t0 @test xv[1,:] == x0 @test size(tvS) == (2,) @@ -1039,11 +1166,12 @@ import Pkg #testing 0-th order root-finding with time ranges/vectors tvr = [t0, Tend/2, Tend, 3Tend/2, 2Tend, 5Tend/2, 3Tend] - taylorinteg(pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1) - @test_throws AssertionError taylorinteg(pendulum!, g, x0, view(tvr, :), - _order, _abstol, maxsteps=1, eventorder=_order+1) - xvr, tvSr, xvSr, gvSr = taylorinteg(pendulum!, g, x0, view(tvr, :), - _order, _abstol, maxsteps=1000) + @test_logs (Warn, max_iters_reached()) taylorinteg( + pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1) + @test_throws AssertionError taylorinteg( + pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1, eventorder=_order+1) + xvr, tvSr, xvSr, gvSr = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000)) @test xvr[1,:] == x0 @test size(tvSr) == (2,) @test norm(tvSr-tvr[3:2:end-1], Inf) < 1e-13 @@ -1071,45 +1199,48 @@ import Pkg end end nothing - end - ) - newex = TaylorIntegration._make_parsed_jetcoeffs(ex) + end) + + newex1, newex2 = TI._make_parsed_jetcoeffs(ex) # Ignore declarations in the function - @test newex.args[1] == :( + @test newex1.args[1] == :( TaylorIntegration.jetcoeffs!(::Val{f!}, t::Taylor1{_T}, q::AbstractArray{Taylor1{_S}, _N}, - dq::AbstractArray{Taylor1{_S}, _N}, p) where + dq::AbstractArray{Taylor1{_S}, _N}, p, __tmpTaylor, __arrTaylor) where {_T <: Real, _S <: Number, _N}) # Include not recognized functions as they appear - @test newex.args[2].args[2] == :(aa = my_simple_function(q, p, t)) - @test newex.args[2].args[5].args[2].args[2] == - :(aa = my_simple_function(q, p, t)) + @test newex1.args[2].args[2] == :(__tmpTaylor[1] = my_simple_function(q, p, t)) + @test newex2.args[2].args[2] == :(aa = my_simple_function(q, p, t)) + @test newex1.args[2].args[5].args[2].args[2] == + :(__tmpTaylor[1] = my_simple_function(q, p, t)) + + # Return line + @test newex1.args[2].args[end] == :(return nothing) + @test newex2.args[2].args[end] == :(return ([aa], [Vector{Taylor1{_S}}(undef, 0)])) # Issue 96: deal with `elseif`s, `continue` and `break` - @test newex.args[2].args[5].args[2].args[3] == - :(for i = 1:length(q) - if i == 1 - TaylorSeries.mul!(dq[i], 2, q[i], ord) - else - if i == 2 - TaylorSeries.identity!(dq[i], q[i], ord) - else - if i == 3 - TaylorSeries.identity!(dq[i], aa, ord) - continue - else - dq[i] = my_complicate_function(q) - break - end - end - end - end - ) + @test newex1.args[2].args[5].args[2].args[3] == :( + for i = 1:length(q) + if i == 1 + TaylorSeries.mul!(dq[i], 2, q[i], ord) + else + if i == 2 + TaylorSeries.identity!(dq[i], q[i], ord) + else + if i == 3 + TaylorSeries.identity!(dq[i], __tmpTaylor[1], ord) + continue + else + dq[i] = my_complicate_function(q) + break + end + end + end + end) end - @testset "Test @taylorize with @threads" begin @taylorize function f1!(dq, q, params, t) for i in 1:10 @@ -1156,10 +1287,27 @@ import Pkg dx03 = similar(x01) t3 = deepcopy(t1) - TaylorIntegration.jetcoeffs!(f1!, t1, x01, dx01, xaux1, nothing) - TaylorIntegration.jetcoeffs!(Val(f1_parsed!), t1p, x01p, dx01p, nothing) - TaylorIntegration.jetcoeffs!(Val(f2!), t2, x02, dx02, nothing) - TaylorIntegration.jetcoeffs!(Val(f3!), t3, x03, dx03, nothing) + # No error is thrown + parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!( + true, f1!, t1, x01, dx01, nothing)) + @test parse_eqs + @test typeof(tmpTaylor) == Vector{Taylor1{Float64}} + @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}} + parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!( + true, f1_parsed!, t1p, x01p, dx01p, nothing)) + @test parse_eqs + @test typeof(tmpTaylor) == Vector{Taylor1{Float64}} + @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}} + parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!( + true, f2!, t2, x02, dx02, nothing)) + @test parse_eqs + @test typeof(tmpTaylor) == Vector{Taylor1{Float64}} + @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}} + parse_eqs, tmpTaylor, arrTaylor = (@test_logs min_level=Warn TI._determine_parsing!( + true, f3!, t3, x03, dx03, nothing)) + @test parse_eqs + @test typeof(tmpTaylor) == Vector{Taylor1{Float64}} + @test typeof(arrTaylor) == Vector{Vector{Taylor1{Float64}}} @test x01 == x01p @test x01p == x02 @@ -1233,23 +1381,23 @@ import Pkg @show Threads.nthreads() - TaylorIntegration.jetcoeffs!(Val(harmosc1dchain!), t, x, dx, μ) - @time TaylorIntegration.jetcoeffs!(Val(harmosc1dchain!), t, x, dx, μ) - - TaylorIntegration.jetcoeffs!(Val(harmosc1dchain_threads!), t_, x_, dx_, μ) - @time TaylorIntegration.jetcoeffs!(Val(harmosc1dchain_threads!), t_, x_, dx_, μ) + parse_eqs, _, _ = (@test_logs min_level=Warn TI._determine_parsing!( + true, harmosc1dchain!, t, x, dx, μ)) + @test parse_eqs + parse_eqs_, _, _ = (@test_logs min_level=Warn TI._determine_parsing!( + true, harmosc1dchain!, t, x_, dx_, μ)) + @test parse_eqs @test x == x_ @test dx == dx_ - tv, xv = taylorinteg(harmosc1dchain!, x0, t0, 10000.0, _order, _abstol, μ) - tv_, xv_ = taylorinteg(harmosc1dchain_threads!, x0, t0, 10000.0, _order, _abstol, μ) + tv, xv = (@test_logs min_level=Warn taylorinteg(harmosc1dchain!, x0, t0, 10000.0, _order, _abstol, μ)) + tv_, xv_ = (@test_logs min_level=Warn taylorinteg(harmosc1dchain_threads!, x0, t0, 10000.0, _order, _abstol, μ)) @test tv == tv_ @test xv == xv_ end - # Issue 106: allow calls to macro from Julia packages @testset "Test @taylorize use in modules/packages" begin # TestPkg is a local (unregistered) Julia package which tests the use @@ -1260,9 +1408,9 @@ import Pkg # are equivalent to (nex_1, nex_2, nex_3) generated here Pkg.develop( Pkg.PackageSpec( path=joinpath(@__DIR__, "TestPkg") ) ) using TestPkg - nex1_ = TaylorIntegration._make_parsed_jetcoeffs(TestPkg.ex1) - nex2_ = TaylorIntegration._make_parsed_jetcoeffs(TestPkg.ex2) - nex3_ = TaylorIntegration._make_parsed_jetcoeffs(TestPkg.ex3) + nex1_, nall1_ = TI._make_parsed_jetcoeffs(TestPkg.ex1) + nex2_, nall2_ = TI._make_parsed_jetcoeffs(TestPkg.ex2) + nex3_, nall3_ = TI._make_parsed_jetcoeffs(TestPkg.ex3) @test length(TestPkg.nex1.args[2].args) == length(nex1_.args[2].args) @test length(TestPkg.nex2.args[2].args) == length(nex2_.args[2].args) @test length(TestPkg.nex3.args[2].args) == length(nex3_.args[2].args)