Make sure coeff_table is not changed by differentiate or integrate (#319

)

* Make a copy so coeff_table is not changed by differentiate or integrate

* Add at_isonethread to avoid some overhead

due to copying the coeff_table in differentiate and integrate

* Add tests and document at_isonethread

* Adapt ci.yml to use more threads

* Correct indentation in ci.yml

* Fixes in ci.yml

* Bump minor version

.github/workflows/ci.yml

@@ -14,8 +14,11 @@ on:
 
 jobs:
   test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.julia-threads }} thread(s) - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     if: "!contains(github.event.head_commit.message, 'skip ci')"
+    env:
+      JULIA_NUM_THREADS: ${{ matrix.julia-threads }}
     strategy:
       fail-fast: false
       matrix:
@@ -26,15 +29,22 @@ jobs:
         # exclude:
         #   - os: macOS-latest
         #     julia-arch: x86
+        julia-threads:
+          - '1'
+        include:
+          - os: ubuntu-latest
+            julia-version: '1'
+            julia-arch: x64
+            julia-threads: '2'
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - name: "Set up Julia"
         uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.julia-version }}
           arch: ${{ matrix.julia-arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v3
         env:
           cache-name: cache-artifacts
         with:

Project.toml

@@ -1,9 +1,10 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.13.2"
+version = "0.14.0"
 
 [deps]
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"

docs/src/api.md

@@ -96,6 +96,7 @@ _dict_unary_ops
 _dict_binary_calls
 _dict_unary_calls
 _dict_binary_ops
+@isonethread
 ```
 
 ## Index

src/auxiliary.jl

@@ -336,3 +336,25 @@ nonlinear_polynomial(a::AbstractSeries) = a - constant_term(a) - linear_polynomi
 nonlinear_polynomial(a::Vector{T}) where {T<:Number} = nonlinear_polynomial.(a)
 
 nonlinear_polynomial(a::Number) = zero(a)
+
+
+"""
+    @isonethread (expr)
+
+Internal macro used to check the number of threads in use, to prevent a data race
+that modifies `coeff_table` when using `differentiate` or `integrate`; see
+https://github.com/JuliaDiff/TaylorSeries.jl/issues/318.
+
+This macro is inspired by the macro `@threaded`; see https://github.com/trixi-framework/Trixi.jl/blob/main/src/auxiliary/auxiliary.jl;
+and https://github.com/trixi-framework/Trixi.jl/pull/426/files.
+"""
+macro isonethread(expr)
+    return esc(quote
+        if Threads.nthreads() == 1
+            $(expr)
+        else
+            copy($(expr))
+        end
+    end)
+end
+

src/calculus.jl

@@ -135,11 +135,11 @@ function differentiate(a::HomogeneousPolynomial, r::Int)
     @inbounds num_coeffs = size_table[a.order+1]
 
     @inbounds for i = 1:num_coeffs
-        iind = coeff_table[a.order+1][i]
+        iind = @isonethread coeff_table[a.order+1][i]
         n = iind[r]
         n == 0 && continue
         iind[r] -= 1
-        kdic = in_base(get_order(),iind)
+        kdic = in_base(get_order(), iind)
         pos = posTb[kdic]
         coeffs[pos] = n * a[i]
         iind[r] += 1
@@ -364,7 +364,7 @@ function integrate(a::HomogeneousPolynomial, r::Int)
     coeffs = zeros(T, size_table[a.order+2])
 
     @inbounds for i = 1:num_coeffs
-        iind = coeff_table[a.order+1][i]
+        iind = @isonethread coeff_table[a.order+1][i]
         n = iind[r]
         n == order_max && continue
         iind[r] += 1

test/manyvariables.jl

@@ -255,6 +255,7 @@ end
     vr = rand(2)
     @test hp(vr) == evaluate(hp, vr)
 
+    ctab = copy(TS.coeff_table)
     @test integrate(yH,1) == integrate(xH, :x₂)
     p = (xT-yT)^6
     @test integrate(differentiate(p, 1), 1, yT^6) == p
@@ -270,6 +271,7 @@ end
     @test integrate(xT^17, 1, yT) == yT
     @test integrate(xT^17, 1, 2.0) == TaylorN(2.0, get_order())
     @test integrate(xT^17, :x₁, 2.0) == TaylorN(2.0, get_order())
+    @test ctab == TS.@isonethread(TS.coeff_table)
     @test_throws AssertionError integrate(xT, 1, xT)
     @test_throws AssertionError integrate(xT, :x₁, xT)
     @test_throws AssertionError differentiate(xT, (1,))
@@ -773,3 +775,67 @@ end
     @test integrate(y, 2) == x * y
 
 end
+
+@testset "Consistency of coeff_table (differentiate)" begin
+    order = 20
+    x, y, z, w = set_variables(Int128, "x y z w", numvars=4, order=2order)
+    ctab = deepcopy(TS.coeff_table);
+
+    function fun(degree::Int)
+        s = x + y + z + w + 1
+        return s^degree
+    end
+
+    function diffs1(f)
+        f1 = differentiate(f, 1)
+        f2 = differentiate(f, 2)
+        f3 = differentiate(f, 3)
+        f4 = differentiate(f, 4)
+        return f1 + f2 + f3 + f4
+    end
+
+    function diffs2(f)
+        vder = zeros(typeof(f), get_numvars())
+        Threads.@threads for i = 1:get_numvars()
+            vder[i] = differentiate(f, i)
+        end
+        return sum(vder)
+    end
+
+    f = fun(order);
+    df_exact = 4*20*fun(order-1);
+    fint_exact = 4*fun(order+1)/(order+1);
+
+    df1 = diffs1(f);
+    @test ctab == TS.coeff_table
+    @test df1 == df_exact
+
+    df2 = diffs2(f);
+    @test ctab == TS.coeff_table
+    @test df1 == df_exact
+
+    # function integ1(f)
+    #     f1 = integrate(f, 1)
+    #     f2 = integrate(f, 2)
+    #     f3 = integrate(f, 3)
+    #     f4 = integrate(f, 4)
+    #     return f1 + f2 + f3 + f4
+    # end
+    #
+    # function integ2(f)
+    #     vder = zeros(typeof(f), get_numvars())
+    #     Threads.@threads for i = 1:get_numvars()
+    #         vder[i] = integrate(f, i)
+    #     end
+    #     return sum(vder)
+    # end
+    #
+    # ii1 = integ1(f);
+    # @test ctab == TS.coeff_table
+    # @test ii1 == fint_exact
+    #
+    # ii2 = integ2(f);
+    # @test ctab == TS.coeff_table
+    # @test ii2 == fint_exact
+
+end