From e32cd4d9295fc8988d0714fb6232a30d82c06dc4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Mon, 22 Jul 2024 22:49:42 +0200 Subject: [PATCH] Updates due to JuliaDiff/TaylorSeries.jl#361 (#194) * Updates due to JuliaDiff/TaylorSeries.jl#361 * Remove comments and small fix * Update CI (use TaylorSeries branch from fork); use `diffeq!` in non-parsed jetcoeffs! methods * Update CI * Use diffeq! in preamble * Retain only `local` declarations in preamble * Trigger CI * Update CI * Update Project.toml * Update CI * Rename diffeq! -> ode! * Update ci.yml * Trigger CI * Update README.md * Rename ode! -> solcoeff! * Bump patch version --- .github/workflows/ci.yml | 6 ++++ Project.toml | 4 +-- README.md | 2 +- src/TaylorIntegration.jl | 23 +++++++++++++ src/integrator.jl | 7 ++-- src/parse_eqs.jl | 72 ++++++++++++++++++++++++++-------------- test/taylorize.jl | 8 ++--- 7 files changed, 85 insertions(+), 37 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8e46312b1..66b2f5e03 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -5,6 +5,10 @@ on: - 'LICENSE.md' - 'README.md' pull_request: + branches: + - main + tags: '*' + jobs: pre_job: runs-on: ubuntu-latest @@ -54,6 +58,8 @@ jobs: ${{ runner.os }}-test-${{ env.cache-name }}- ${{ runner.os }}-test- ${{ runner.os }}- + ### Uncomment line below to test using a dependency branch + # - run: julia --project=. -e 'import Pkg; Pkg.add(url="https://github.com/PerezHz/TaylorSeries.jl", rev="jp/evaluate-tn"); Pkg.instantiate()' - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 diff --git a/Project.toml b/Project.toml index 443f77e57..dd0d5da43 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.15.2" +version = "0.15.3" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" @@ -38,7 +38,7 @@ RecursiveArrayTools = "2, 3" Reexport = "1" Requires = "1" StaticArrays = "0.12.5, 1" -TaylorSeries = "0.17" +TaylorSeries = "0.18" Test = "<0.0.1, 1" julia = "1.6" diff --git a/README.md b/README.md index becb5f832..eee64d9ff 100644 --- a/README.md +++ b/README.md @@ -44,4 +44,4 @@ Comments, suggestions, contributions and improvements are welcome and appreciate ## Acknowledgments -We acknowledge financial support from DGAPA-PAPIIT grants IG-100616 and IG-100819. +We acknowledge financial support from DGAPA-PAPIIT grants IG-100616, IG-100819 and IG-101122. diff --git a/src/TaylorIntegration.jl b/src/TaylorIntegration.jl index a6134a524..0686b00d7 100644 --- a/src/TaylorIntegration.jl +++ b/src/TaylorIntegration.jl @@ -28,4 +28,27 @@ function __init__() end end +@inline function solcoeff!(a::Taylor1{T}, b::Taylor1{T}, k::Int) where {T<:TaylorSeries.NumberNotSeries} + a[k] = b[k-1]/k + return nothing +end + +@inline function solcoeff!(res::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, + k::Int) where {T<:TaylorSeries.NumberNotSeries} + for l in eachindex(a[k-1]) + res[k][l] = a[k-1][l]/k + end + return nothing +end + +@inline function solcoeff!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + k::Int) where {T<:TaylorSeries.NumberNotSeries} + for l in eachindex(a[k-1]) + for m in eachindex(a[k-1][l]) + res[k][l][m] = a[k-1][l][m]/k + end + end + return nothing +end + end #module diff --git a/src/integrator.jl b/src/integrator.jl index d19b8204c..06fee0965 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -27,14 +27,13 @@ function jetcoeffs!(eqsdiff::Function, t::Taylor1{T}, x::Taylor1{U}, params) whe ordnext = ord+1 # # Set `taux`, `xaux`, auxiliary Taylor1 variables to order `ord` - # @inbounds taux = Taylor1( t.coeffs[1:ordnext] ) @inbounds xaux = Taylor1( x.coeffs[1:ordnext] ) # Equations of motion dx = eqsdiff(xaux, params, t) # Recursion relation - @inbounds x[ordnext] = dx[ord]/ordnext + @inbounds solcoeff!(x, dx, ordnext) end nothing end @@ -67,8 +66,6 @@ function jetcoeffs!(eqsdiff!::Function, t::Taylor1{T}, for ord in 0:order-1 ordnext = ord+1 - # # Set `taux`, auxiliary Taylor1 variable to order `ord` - # @inbounds taux = Taylor1( t.coeffs[1:ordnext] ) # Set `xaux`, auxiliary vector of Taylor1 to order `ord` for j in eachindex(x) @inbounds xaux[j] = Taylor1( x[j].coeffs[1:ordnext] ) @@ -79,7 +76,7 @@ function jetcoeffs!(eqsdiff!::Function, t::Taylor1{T}, # Recursion relations for j in eachindex(x) - @inbounds x[j].coeffs[ordnext+1] = dx[j].coeffs[ordnext]/ordnext + @inbounds solcoeff!(x[j], dx[j], ordnext) end end nothing diff --git a/src/parse_eqs.jl b/src/parse_eqs.jl index db082cb13..c4f85979a 100644 --- a/src/parse_eqs.jl +++ b/src/parse_eqs.jl @@ -35,11 +35,12 @@ mutable struct BookKeeping v_array4 :: Vector{Symbol} v_preamb :: Vector{Union{Symbol,Expr}} retvar :: Symbol + numaux :: Int function BookKeeping() return new(Dict{Symbol, Expr}(), Dict{Union{Symbol,Expr}, Number}(), Dict{Symbol, Expr}(), Symbol[], Symbol[], Symbol[], Symbol[], Symbol[], Symbol[], - Union{Symbol,Expr}[], :nothing) + Union{Symbol,Expr}[], :nothing, 0) end end @@ -123,7 +124,7 @@ const _HEAD_ALLOC_TAYLOR1_VECTOR = sanitize(:( # Constants for the initial declaration and initialization of arrays const _DECL_ARRAY = sanitize( Expr(:block, :(__var1 = Array{Taylor1{_S}}(undef, __var2)), - :( for i in CartesianIndices(__var1) __var1[i] = Taylor1( zero(constant_term(__x[1])), order ) end )) + :( for i in eachindex(__var1) __var1[i] = Taylor1( zero(constant_term(__x[1])), order ) end )) ); @@ -147,7 +148,7 @@ function _make_parsed_jetcoeffs(ex::Expr) new_jetcoeffs, new_allocjetcoeffs = _newhead(fn, fnargs) # Transform the graph representation of the body of the functions: - # defspreamble: inicializations used for the zeroth order (preamble) + # defspreamble: initializations 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 @@ -155,17 +156,17 @@ function _make_parsed_jetcoeffs(ex::Expr) defspreamble, defsprealloc, fnbody, bkkeep = _preamble_body(fnbody, fnargs) # 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) + # rec_fnbody: recursion loop for the body-function (recursion loop for all orders) + rec_fnbody = _recursionloop(fnargs, bkkeep) # 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)) ) - # Add rec_fnbody to forloopblock + forloopblock = Expr(:for, :(ord = 0:order-1), Expr(:block, :(ordnext = ord + 1)) ) + + # Add rec_fnbody to `forloopblock` push!(forloopblock.args[2].args, fnbody.args[1].args..., rec_fnbody) # Add preamble and recursion body to `new_jetcoeffs` - push!(new_jetcoeffs.args[2].args, defspreamble..., rec_preamb) + push!(new_jetcoeffs.args[2].args, defspreamble...) # Push preamble and forloopblock to `new_jetcoeffs` and return line push!(new_jetcoeffs.args[2].args, forloopblock, Meta.parse("return nothing")) @@ -337,11 +338,11 @@ of the original diferential equations function. """ function _preamble_body(fnbody, fnargs) - # Inicialize BookKeeping struct + # Initialize BookKeeping struct bkkeep = BookKeeping() # Rename vars to have the body in non-indexed form; bkkeep has different entries - # for bookkeeping variables/symbolds, including indexed ones + # for bookkeeping variables/symbols, including indexed ones fnbody, bkkeep.d_indx = _rename_indexedvars(fnbody) # Create `newfnbody` which corresponds to `fnbody`, cleaned (without irrelevant comments) @@ -362,15 +363,22 @@ function _preamble_body(fnbody, fnargs) preamble = subs(preamble, bkkeep.d_assign) prealloc = subs(prealloc, bkkeep.d_assign) - # Include the assignement of indexed auxiliary variables + # Include the assignment of indexed auxiliary variables defsprealloc = _defs_allocs!(prealloc, fnargs, bkkeep, false) preamble = subs(preamble, bkkeep.d_indx) - defspreamble = Expr[preamble.args...] + + # Retain only `local` declarations from `preamble` in `new_preamble` + new_preamble = Expr(:block,) + for (i, arg) in enumerate(preamble.args) + (arg.head == :local) && push!(new_preamble.args, arg) + end + defspreamble = Expr[new_preamble.args...] + # Bring back substitutions newfnbody = subs(newfnbody, bkkeep.d_indx) # 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] + bkkeep.retvar = length(fnargs) == 3 ? subs(bkkeep.v_newvars[end-bkkeep.numaux], bkkeep.d_indx) : fnargs[1] return defspreamble, defsprealloc, newfnbody, bkkeep end @@ -834,7 +842,7 @@ function _replacecalls!(bkkeep::BookKeeping, fnold::Expr, newvar::Symbol) end elseif ll == 3 - # Binary call; no auxiliary expressions needed + # Binary call newarg2 = fnold.args[3] # Replacements @@ -859,6 +867,28 @@ function _replacecalls!(bkkeep::BookKeeping, fnold::Expr, newvar::Symbol) # Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), # :_arg2 => :(constant_term($(newarg2))), :_k => :ord)) + # Auxiliary expression + if aux_fnexpr.head != :nothing + newaux = genname() + + 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, + :(TaylorSeries.zero!(_res)), + :(_res.coeffs[1] = $(aux_fnexpr.args[2])) ) + aux_fnexpr = subs(aux_fnexpr, + Dict(:_res => newaux, :_arg1 => :(constant_term($(newarg1))), :_aux => newaux)) + + fnexpr = subs(fnexpr, Dict(:_aux => newaux)) + if newvar ∈ bkkeep.v_arraydecl + push!(bkkeep.v_arraydecl, newaux) + else + bkkeep.numaux += 1 + push!(bkkeep.v_newvars, newaux) + end + end else # Recognized call, but not a unary or binary call; copy expression fnexpr = :($newvar = $fnold) @@ -1037,19 +1067,13 @@ 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_fnbody = sanitize( :( TaylorIntegration.solcoeff!($(fnargs[1]), $(bkkeep.retvar), ordnext) ) ) elseif ll == 4 bkkeep.retvar = fnargs[1] - rec_preamb = sanitize(:( - for __idx in eachindex($(fnargs[2])) - $(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] = - $(bkkeep.retvar)[__idx].coeffs[ordnext]/ordnext + TaylorIntegration.solcoeff!($(fnargs[2])[__idx], $(bkkeep.retvar)[__idx], ordnext) end)) else @@ -1057,7 +1081,7 @@ function _recursionloop(fnargs, bkkeep::BookKeeping) "Wrong number of arguments in the definition of the function $fn")) end - return rec_preamb, rec_fnbody + return rec_fnbody end diff --git a/test/taylorize.jl b/test/taylorize.jl index f438ad30e..2dfffbcf0 100644 --- a/test/taylorize.jl +++ b/test/taylorize.jl @@ -221,7 +221,7 @@ import Logging: Warn end - # Pendulum integrtf = 100.0 + # Pendulum integration @testset "Integration of the pendulum and DiffEqs interface" begin @taylorize function pendulum!(dx::Array{T,1}, x::Array{T,1}, p, t) where {T} dx[1] = x[2] @@ -1297,10 +1297,8 @@ import Logging: Warn # Include not recognized functions as they appear @test newex1.args[2].args[2] == :(aa = __ralloc.v0[1]) - @test newex1.args[2].args[3] == :(aa = my_simple_function(q, p, t)) @test newex2.args[2].args[2] == :(aa = my_simple_function(q, p, t)) - @test newex1.args[2].args[6].args[2].args[2] == - :(aa = my_simple_function(q, p, t)) + @test newex1.args[2].args[3].args[2].args[2] == :(aa = my_simple_function(q, p, t)) # Return line @test newex1.args[2].args[end] == :(return nothing) @@ -1331,7 +1329,7 @@ import Logging: Warn end end) - @test newex1.args[2].args[6].args[2].args[3] == Base.remove_linenums!(ex) + @test newex1.args[2].args[3].args[2].args[3] == Base.remove_linenums!(ex) # Throws no error ex = :(