Precision–agnostic function with powered mathematical constants generated from symbolic SymPy expression

General Usage
1

Hello Julia users,

I use Mathematica to generate analytic expressions for some quite complicated integrals I calculate.

In Mathematica I generate tables of indexes and symbolic expression e.g.

{{-1, -1, -1, -1, -1, -1, 0}, {(256*Pi^4*(Log[b] - Log[a + b] - Log[b + d]))/(a^2*c^2*d^2)}}
{{-1, -1, -1, -1, -1, 0, -1}, {(256*Pi^4*(Log[b] - Log[a + b] + Log[a + b + c] + Log[a + b + d]))/(a^2*c^2*d^2)}}
{{-1, -1, -1, -1, 0, -1, -1}, {(256*Pi^4*(Log[b] - 2*Log[b + c] - Log[b + d] + Log[b + c + d]))/(a^2*c^2*d^2)}}
{{-1, -1, -1, -1, 0, -1, 0}, {(-256*(b^2 - c*d)*Pi^4)/(a^2*b*c^2*(b + c)*d^2*(b + d))}}
{{-1, -1, -1, -1, 0, 0, -1}, {(256*Pi^4)/(a^2*b*c^2*d^2)}}

and save it to the text file (in production the formulas are often much longer, I chose some of the shortest ones for this MWE).

Then, using SymPy.jl and some manual parsing I create function to compute values of integral with given indexes ({i,j,k,l,m,n,p} that are integers) and parameters (a,b,c,d which are floats).

    f = open("I4O2s.txt")
    lines = readlines(f)
    using SymPy
    @vars a b c d
    const sympy_parsing_mathematica = SymPy.PyCall.pyimport("sympy.parsing.mathematica")
    # create new array with pairs as arguments
    A = []
    for i in 1:length(lines)
          push!(A, (tuple(parse.(Int, split(split(strip(lines[i], ['{', '}']), "}, {")[1], ","))...), SymPy.lambdify(sympy_parsing_mathematica."mathematica"(split(strip(lines[i], ['{', '}']), "}, {")[2]), (a,b,c,d))))
     end   
    # create dict from array 
    i4dict = Dict(A)
    A = 0
    # use dict to create function
    i4(i,j,k,l,m,n,p,a,b,c,d) = i4dict[(i,j,k,l,m,n,p)](a,b,c,d)

It may not be elegant, but it works (all sugestions of improvement are more than welcome!).

Offtopic question: How to view Julia expression from SymPy.lambdify?

One, a bit off-topic question would be how to look–up the Julia expression for the givent formula, since after SymPy.lambdify i only get:

x = SymPy.lambdify(sympy_parsing_mathematica."mathematica"(split(strip(lines[6], ['{', '}']), "}, {")[2]), (a,b,c,d))

#79 (generic function with 1 method)

when creating dictionary i get expressions such as ##79#81(Box(##418)) instead of formula

 i4dict = Dict(A)
Dict{NTuple{7,Int64},getfield(SymPy, Symbol("##79#81"))} with 36 entries:
  (-1, -1, 0, 0, -1, -1, -1)   => ##79#81(Box(##418))
  (-1, -1, -1, 0, -1, -1, 0)   => ##79#81(Box(##410))
  etc...

and trying to view it in the dictionary gets me:

i4dict(-1,-1,-1,-1,-1,-1,-1)

ERROR: MethodError: objects of type Dict{NTuple{7,Int64},getfield(SymPy, Symbol("##79#81"))} are not callable
Stacktrace: [1] top-level scope at none:0

My problem is following:
When I try to compute the expression e.g.
i4(-1, -1, -1, -1, 0, 0, -1, BigFloat(1.0), BigFloat(1.0), BigFloat(1.0), BigFloat(1.0))
i get different result than when I do it with higher precision in Mathematica (and different from the article with some high-precision values for these expressions, they agree with the ones calculated in Mathematica).

I managed to narrow it down to the following: the problem is with the precision of pi^4. After some tests I realised that the pi is converted to Float64 and then powered. Of course this is expected behaviour but it leads to precision loss. However, due to numerical instabilities in some of the more complicated expressions I would need extended precision: DoubleDouble, Quad or BigFloats (or ArbFloats, maybe I will even think about some port of quad double and/or double quad to Julia in the future but BigFloats are probably good enough for now).

So my question is: is there a nice, Julian way to generate function that computes with arbitrary precision floating, from symbolic expression, with automatically adjusted precision of the constant (pi in this case, but I believe the question is more broad and applies to many mathematical constans in Base.MathConstants).

In my case I can collect pi^4 using mathematica and manually multiply the function by BigFloat(pi)^4 or Float128(pi)^4 (or taking the inherited type from a,b,c,d arguments) but I am asking for general case, when in principle pi or/and e could be in various powers and therefore such manual manipulation would be infeasible. I believe that there is a better way to do it and therefore asking this question.

I would also like to add that it is amazing that such things are possible in Julia and thank all autors for this amazing language :smiley:

Why don't you use SymPy/Reduce from the start to generate expressions?

I would prefer to use open source tool for generation of formulas, but at the moment it is not really possible:
– Simplification of expressions is unfortunately much better in Mathematica,
– Mathematica has better support for some special functions that sometimes appear in these formulas,
– I know it better and my collaborators use it.

Hopefully latter I will be able to move to SymPy and FORM (https://github.com/HaraldHofstaetter/FORM.jl) but at the moment it is not possible and I would need the results from Mathematica for benchmarks anyway.

I am also looking forward to learn and try Reduce since there’s a great Julia package, but the same problems as above apply.

2

first, welcome to the Julia Language Forum!

from what i know, pi can be calculated with arbitrary precision in julia, to use as a number, you have to give pi a concrete type or use it in operations with concrete types (Float64, BigFloat). for example, BigFloat(pi) has all the correct digits. this post has a deeper treatment of the subject

2 Likes
3

Thank You very much for the link. I read it and it is quite didactic :slight_smile: Unfortunately I was not able to solve my problem using these informations.

Maybe I will try to clarify a bit what my problem is. Key point in my question is that the expressions are very long and generated with Computer Algebra System (in my case there are few hundreds but in the future there can be much more).

Therefore, I would like to automate code generation as far as possible (also for reproduction and error reasons) and with that amount of expressions it is nearly impossible to manually write concrete type in every expression. I am looking for a way to do it in some automated way, possibly using the existing workflow, but I’m open for other ways to do it.

4

What does Mathematica’s Pi get converted to via your parsing code? It is the symbolic SymPy.PI or Julia’s pi::Irrational? I ask because:

julia> using SymPy

julia> 256*PI^4 |> N |> typeof
BigFloat

julia> 256*pi^4 |> typeof
Float64

As the blog post above shows, you can convert pi to a BigFloat and get all those bits accurately (e.g. 256*big(pi)^4) but it sounds like you’d have to do some parsing to replace the pi’s with big(pi). It seems like however if you could just parse Mathematica’s Pi’s into SymPy’s PI, then it might just work out (since PI gets evaluated to a BigFloat automatically, instead of being converted based on the type of the other argument in * or + like Julia’s pi).

2 Likes
5

Thank You very much for answer.

I checked that mathematica Pi gets converted to something weird:

julia> sympy_parsing_mathematica."mathematica"("Pi")
Π

julia> typeof(sympy_parsing_mathematica."mathematica"("Pi"))
Sym

julia> sympy_parsing_mathematica."mathematica"("Pi") == SymPy.PI
false

julia> sympy_parsing_mathematica."mathematica"("Pi") == pi
false

julia> sympy_parsing_mathematica."mathematica"("Pi") |> N |> typeof
Sym

I suppose this is just a symbol Pi.

What I managed to do is I used SymPy abilities to add Pi as new symbol and transform it to SymPy.PI as You suggested. Thank You very much! Here is the code that follows Your advice:

julia> using SymPy
julia> @vars Pi
julia> sympy_parsing_mathematica."mathematica"("Pi^4").subs(Pi, SymPy.PI) |> N |> typeof
BigFloat

Unfortunately this way of code evaluation is too slow (i.e. 0.001350 vs 0.000016 s for simple 256*Pi^4 function calculated with SymPy’s N and native Julia generated by lambdify function).

However, I came to the conclusion that it is the “problem” lies in the lambdify function. When the function is lambdified SymPy.PI is converted to :pi and the code is effectively 256*pi^4 which turns Pi into Float64. At the moment I don’t see any workaround using SymPy and automatic expressions generation in Julia.

I will try to work around this by manually factoring out the pi^4 from all expressions and multiplying the expressions by Type(pi)^4 in function definition.

Thanks for help :slight_smile:

BTW I have found following problem with evaluation of the values:

Is this a bug in SymPy? 
julia> N(SymPy.PI, 50)
ERROR: UndefVarError: is_integer not defined
1 Like
6

Thanks, this issue brought about two bug fixes (https://github.com/JuliaPy/SymPy.jl/pull/285) that help:

  • N(PI, 50) now works
  • lambdify was a bit broken, thought your case didn’t hit it.

In general, you may benefit from the following pattern:

@vars Pi x  # Pi is what sympy_parsing_mathematica."mathematica"("Pi") returns, a symbol
ex = 256 * Pi * x # some expression

To lambdify this, you can call lambidfy. Digging in gives you the control you need though.

The expression lambdified is given by:

SymPy.walk_expression(ex)  # :(256 * pi * x)

The lambdified function then will be x -> 256 * pi * x and the float comes from the 256 * pi conversion. You need to treat 256 as a big int. Enter, ChangePrecision.jl:

using ChangePrecision
body = SymPy.walk_expression(ex)
cpbody = ChangePrecision.changeprecision(BigFloat, body)
vars = (x,)
fn = eval(Expr(:function, Expr(:call, gensym(), Symbol.(vars)...), cpbody))
fn(1) # use BigFloat now

I didn’t test timings, but I don’t see why it would be slow.

1 Like