# SymPy/Symbolics and HCubature

**URL:** <https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513>\
**Category:** New to Julia\
**Tags:** symbolic\
**Created:** [March 4, 2021, 11:24pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513 "2021-03-04T23:24:32Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 4, 2021, 11:24pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/1 "2021-03-04T23:24:32Z")

</div>

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:

```julia
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.,

```julia
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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 4, 2021, 11:43pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/2 "2021-03-04T23:43:06Z")

</div>

> [@Aleksandr\_Mikheev](#):
>
> **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?

See [`SymPy.lambdify`](https://juliahub.com/docs/SymPy/KzewI/1.0.41/Tutorial/basic_operations/#lambdify-1).

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 5, 2021, 12:11am UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/3 "2021-03-05T00:11:07Z")

</div>

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

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

```

instead of

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

```

I also tried to first make an array of variables,

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 5, 2021, 1:30am UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/4 "2021-03-05T01:30:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 5, 2021, 2:04am UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/5 "2021-03-05T02:04:03Z")

</div>

> [@Aleksandr\_Mikheev](#):
>
> However, `g = lambdify(f(x,y))` returns a function that takes a tuple (x,y)(x,y) as an argument, not a vector, i.e., if I understand correctly, it gives `function g(x,y)`

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 5, 2021, 2:25am UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/6 "2021-03-05T02:25:57Z")

</div>

Got to a computer to write the example:

```julia
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})`

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 5, 2021, 2:49pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/7 "2021-03-05T14:49:34Z")

</div>

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](https://docs.julialang.org/en/v1/manual/variables-and-scoping/#Let-Blocks), let statements allocate a new binding for a variable that does not affect the binding outside the let block. E.g.,

```julia
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)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 5, 2021, 3:04pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/8 "2021-03-05T15:04:16Z")

</div>

> [@Aleksandr\_Mikheev](#):
>
> But what role does it play in your code?

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.

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 5, 2021, 8:09pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/9 "2021-03-05T20:09:47Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 5, 2021, 8:14pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/10 "2021-03-05T20:14:38Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 5, 2021, 8:15pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/11 "2021-03-05T20:15:15Z")

</div>

Thanks a lot again!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 5, 2021, 9:36pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/12 "2021-03-05T21:36:50Z")

</div>

> [@Aleksandr\_Mikheev](#):
>
> 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(?)

“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](https://en.wikipedia.org/wiki/Embarrassingly_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.

---

<div class="post-metadata">

**Author:** ![bad\_at\_math](https://avatars.discourse-cdn.com/v4/letter/b/51bf81/32.png) [@bad\_at\_math](https://discourse.julialang.org/u/bad_at_math)\
**Post date:** [August 30, 2021, 11:03pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/13 "2021-08-30T23:03:50Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 30, 2021, 11:11pm UTC](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513/14 "2021-08-30T23:11:43Z")

</div>

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.