How to generate sums effectively via creating an Add object

[PerformanceCoder]

Hi!

In Python, if I have an array of terms, a sum element can be created in this way:

res = sp.Add(*list_of_terms)

Usually it works much faster than

for i in range(n):
    res += list_of_terms[i]

Is there any possibility to create SymPy sums in Julia the same way?

[jverzani]

So maybe through PyCall, I wouldn't know. This wouldn't be necessary in Julia but might be helpful here, as each + call dispatches to a Python call, so I'd guess performance would be like the sum example. Does SymPy itself provide such a function? If so, just calling that would be the better solution.

[PerformanceCoder]

Thank you for your reply. Created a minimal working example using SymPy via PyCall. It demonstrates the difference in performance between the two ways of sum generation.

using BenchmarkTools

using SymPy
const SP = SymPy

using PyCall
sympy = pyimport("sympy")

function generate_sum_sp_trivial(n::Int)
    SP.@syms x
    sum = 0.0

    for i in 1:n
        sum += x^i*sin(i*x)
    end
    return sum
end

function generate_sum_sp_array(n::Int)
    SP.@syms x
    terms = map(t -> x^t*sin(t*x), 1:n)

    return +(terms...)
end

function generate_sum_sp_pycall(n::Int)
    x = sympy.Symbol("x")
    terms = map(t -> x^t*sympy.sin(t*x), 1:n)
    return sympy.Add(terms...)
end

for i in [1, 10, 100, 1000]
    println("----\n$i\n----\n")
    println("Trivial generation")
    display(@benchmark generate_sum_sp_trivial($i))

    println("\n\nMap and sum")
    display(@benchmark generate_sum_sp_array($i))

    println("\n\nSymPy.Add")
    display(@benchmark generate_sum_sp_pycall($i))
end

In this test, for i = 1000 I can achieve a 40x performance improvement with the use of sympy.Add. Could you please add a wrapper for this function in SymPyCore?

[jverzani]

Sure, what interface. Here is some test of a possible one that seems promising:

 using SymPy
import SymPy: SymbolicObject

Base.:+(x::SymbolicObject, y, z, zs...) = sympy.Add(x, y, z, zs...) # <-- this is the addition
function f1(x)
    t = first(x)
    for i in 2:length(x)
        t += x[i]
    end
    t
end

@syms x[1:50]

julia> @time sum(x); @time sum(x);
  0.064616 seconds (41.41 k allocations: 2.823 MiB, 67.47% compilation time)
  0.000073 seconds (50 allocations: 800 bytes)

julia> @time +(x...); @time +(x...);
  0.090991 seconds (39.29 k allocations: 2.621 MiB, 98.48% compilation time)
  0.000078 seconds (166 allocations: 5.609 KiB)

julia> @time sympy.Add(x...); @time sympy.Add(x...);
  0.000209 seconds (114 allocations: 4.406 KiB)
  0.000080 seconds (114 allocations: 4.406 KiB)

julia> @time f1(x); @time f1(x);
  0.014017 seconds (3.66 k allocations: 249.180 KiB, 99.08% compilation time)
  0.000075 seconds (50 allocations: 800 bytes)

julia> @syms x[1:5000];

julia> @time sum(x); @time sum(x);
  9.424470 seconds (5.00 k allocations: 78.125 KiB)
  9.642760 seconds (5.00 k allocations: 78.125 KiB)

julia> @time +(x...); @time +(x...);
  0.059711 seconds (26.08 k allocations: 1.509 MiB, 32.12% compilation time)
  0.001912 seconds (15.02 k allocations: 743.609 KiB)

julia> @time sympy.Add(x...); @time sympy.Add(x...);
  0.002470 seconds (10.02 k allocations: 587.141 KiB)
  0.001666 seconds (10.02 k allocations: 587.141 KiB)

julia> @time f1(x); @time f1(x);
 69.795652 seconds (5.00 k allocations: 78.125 KiB)
77.302930 seconds (5.00 k allocations: 78.125 KiB)

[jverzani]

I checked in a solution using +(x...). Let me know if that isn't helpful.

[PerformanceCoder]

Thank you, this solution works great!