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

Author: PerezHz

.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

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"
 

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.

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

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

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,25 +148,25 @@ 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
     # bkkeep: book-keeping structure having info of the variables
     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,27 +1067,21 @@ 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
         throw(ArgumentError(
         "Wrong number of arguments in the definition of the function $fn"))
     end
 
-    return rec_preamb, rec_fnbody
+    return rec_fnbody
 end
 
 

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 = :(