SymPy/Symbolics and HCubature

New to Julia
1

Hi everyone! I was not sure whether to post it here or in the Numerics subforum, but given that my question is very basic I decided to go for this section.

Let’s say I want to integrate some function f(x,y) over [a_1,b_1] \times [a_2,b_2]. For definitness, suppose f(x,y) = x^2 + y and a_1 = -1, b_1 = 3, a_2 = 0, b_2 = 2. Using HCubature I could achieve this as follows:

using HCubature

function f(x)
   return x[1]^2 + x[2]
end

xmin = [-1,0];
xmax = [3,2];

val = hcubature(f,xmin,xmax);

Now, suppose the analytic expression for f(x,y) is not as simple and is given by some product of other functions and derivatives of thereof. In that case, it is much more convenient to first get the expression for f(x,y) using some in-built CAS such as SymPy (which I am a bit more familiar with) or Symbolics and then transform it into a Julia-type function. For simplicity, let’s say that the resulting expression is once again x^2 + y, i.e.,

using SymPy

x,y = symbols("x,y", real=true);

f(x,y) = x^2 + y;

How can I now convert it into f(x) from the first example? Thank you in advance!

TL;DR: Given symbolic a f(x,y) is there any way to convert it into the function f(x), with x being a 2-dimensional array?

2

See SymPy.lambdify.

3

Ah, yes, I was actually aware of SymPy.lambdify. However, g = lambdify(f(x,y)) returns a function that takes a tuple (x,y) as an argument, not a vector, i.e., if I understand correctly, it gives

function g(x,y)
   return x^2 + y
end

instead of

function g(x)
   return x[1]^2 + x[2]
end

I also tried to first make an array of variables,

using SymPy

x,y = symbols("x,y", real=true);
v    = [x,y];

f(v) = v[1]^2 + v[2];
g    = lambdify(f(v));

but that didn’t help either.

4

See build_function in Symbolics.jl. It builds the Julia function with an array input if you do build_function(symf,[x,y]), i.e. put the arguments in the style you want the input built.

2 Likes
5

Can’t you just wrap this in another function? For example: g = let gg = lambdify(f(x,y)); x -> gg(x...); end

1 Like
6

Got to a computer to write the example:

using Symbolics
@variables x y
_f = x^2 + y
f = eval(build_function(_f,[x,y]))
f([1.0,2.0])

or, to avoid world-age issues, use f = build_function(_f,[x,y],expression=Val{false})

2 Likes
7

Thanks a lot! Although I have to admit I don’t fully understand what’s going on in your code. From what I understand, lambdify(f(x,y)) returns a lambda (or anonymous) function, which for our example reads (x,y) -> x^2 + y. So, if you introduce gg = lambdify(f(x,y)), then gg(a,b) will return a^2 + b. So far so good. Now, we want to pass vectors as an argument of our function, so instead of (x,y) -> x^2 + y we want x -> x[1]^2 + x[2]. Using the splat operator ... we can define a new lambda function x -> g(x...), which when passing [a,b] to it returns g([a,b]...) = g(a,b) = a^2 + b. So I would imagine that g = v -> lambdify(f(x,y))(v...) should do the job, but it clearly doesn’t: ERROR: MethodError: no method matching ##258(::Int64, ::Int64). Therefore, this let block is necessary and serves some important purpose that I fail to understand. From what I could understand from Scope of Variables · The Julia Language, let statements allocate a new binding for a variable that does not affect the binding outside the let block. E.g.,

x = 4

let x = 6
   println("x = $x");     # prints x = 6
end

println("x = $x");        # prints x = 4

But what role does it play in your code? (Sorry for bothering you with such simple questions)

8

The let block is there so that you lambdify once and cache the result in a local variable gg that is only used within the x -> gg(x...) function and is not visible elsewhere.

3 Likes
9

Thanks! It does work nicely. While we are at it, I was also wondering if one could also convert a symbolic function f(x,y,z) into a Julia-type function that takes [x,y],z as an argument and returns an array of length length(z). For instance, suppose I want to compute g(z) = \int_{a_1}^{b_1}\mathrm{d}x\int_{a_2}^{b_2}\mathrm{d}y f(x,y,z) using HCubature. Since the hcubature function accepts vector-valued functions, it seems suggestive to construct a Julia-type function that f, for a given array z = [z_1,z_2,...,z_n], returns [f(x,y,z_1),f(x,y,z_2),...f(x,y,z_n)]. If that was the case, I could compute g(z) on a grid \lbrace z_i\rbrace in one line instead of doing a for-loop computing it at each z_i individually. Presumably, that would be not only more convenient but also faster (although I am not sure how “vectorized” HCubature is) and more suitable for parallelization if needed(?). Alternatively, z could be definied globally and the desired function then may only take [x,y] as an argument and still return [f(x,y,z_1),f(x,y,z_2),...f(x,y,z_n)] (I would guess that this approach is easier to implement but is less elegant since it involves global variables).

Hope it somewhat makes sense.

10

It matches the form you write. So if you have an array of equations, then it gives you a function that outputs an array of equations. If you have build_function([f1,f2,f3],[x,y],z), then it’s a function f([x,y],z) that outputs an array [f1,f2,f3]. If you make those arrays sparse matrices, etc. I think you can see the pattern from there.

1 Like
11

Thanks a lot again!

12

“Vectorization” is not really the right conceptual model here — scalar code in Julia is fast. It is also may be less parallelizable to vectorize your integrand because the alternative is to compute each integral independently, which is embarassingly parallel.

However, you might get some serial speedup from vector-valued integrands since it e.g. reduces the bookkeeping overhead of the quadrature, under the following conditions:

  1. You should use StaticArrays.jl for your vector-valued output, so that HCubature doesn’t incur heap allocations for working with these vectors. (This only really matters if the vectors are small.)
  2. The different integrands have qualitatively similar behaviors, so that the adaptively refined quadrature mesh is similar for all of the integrands.
  3. Ideally, you can share some computation between the different integrands.

In practice, I would normally only use vector-valued integrands for multiple scalar integrands when they share a lot of computations in common.

1 Like
13

I have a similar use case to @Aleksandr_Mikheev , and this thread has been useful.

Returning to this example: Can @ChrisRackauckas (or anyone else) elaborate on a few things that are unclear to me:

  • Why does calling build_function with expression=Val{true} introduce world-age issues, and is this something that will likely be fixed by later releases of Symbolics.jl?
  • Why does build_function return a function instead of an expression that needs to be eval’d before being executed?
14

Both are the same question. First, that doesn’t introduce world-age issues. Using eval gives world-age issues. But since most people don’t seem to understand that, build_function returns a function that uses a RuntimeGeneratedFunction so that this question doesn’t need to be asked and everything just naturally works for users.