- Join over
**1.5M+ people** - Join over
**100K+ communities** - Free
**without limits** - Create
**your own community**

- Jan 29 2019 01:29chakravala labeled #21
- Jan 29 2019 01:27chakravala closed #10
- Jan 29 2019 01:27chakravala closed #13
- Jan 29 2019 01:25chakravala commented #13
- Jan 29 2019 01:22chakravala labeled #23
- Jan 28 2019 17:39flipgthb starred chakravala/Reduce.jl
- Jan 27 2019 15:59chakravala commented #23
- Jan 27 2019 13:50atthom commented #23
- Jan 26 2019 16:30chakravala commented #23
- Jan 26 2019 01:39atthom opened #23
- Jan 22 2019 11:27chakravala labeled #22
- Jan 22 2019 09:35chakravala commented #22
- Jan 22 2019 09:34chakravala closed #22
- Jan 22 2019 09:34chakravala commented #22
- Jan 22 2019 02:12cscherrer commented #22
- Jan 22 2019 00:02chakravala commented #22
- Jan 21 2019 21:08cscherrer opened #22
- Jan 21 2019 19:53cscherrer commented #21
- Jan 21 2019 19:29chakravala commented #21
- Jan 21 2019 19:23chakravala commented #21

Hi! I'm trying to approximate probability distributions with weighted polynomials weighted with exp(- x'x). Basically, the idea is that polynomials are "conjugate" (multiply two polynomials, and you get another polynomial), and easy to integrate -- two things that could make Bayesian inference very easy. So I'm now trying to write a few things to compare how well it does in accuracy & speed vs MCMC.

The integration is where Reduce comes in. It lets me generate functions for marginalizing out parameters or CDFs without having to write out all the different expressions for the integral of `a*exp(-x^2)`

, `b*x*exp(-x^2)`

, `c*x^2*exp(-x^2)`

, etc....

Two changes I would like to move towards eventually are splines in place of polynomials, and fitting log densities instead of densities. There's way too much math I don't understand before I can get there, so even though both of those things (especially splines) would be a massive improvement, it'll be months before I get there. Just trying a proof of concept now.

I'm not far at all yet, so I don't have much to say at the moment, but I think manipulating Julia expressions is an awesome idea that could be really useful in the future for writing powerful code

Noticed in your issue that you used

`rcall(:(int( x^3 * exp(-x^2), x)))`

. It's not well documented yet, but you can also get the same result by doing this: `int(:(x^3 * exp(-x^2)), :x)`

.
Cool, I like that syntax more.

Definitely looking forward to seeing where it goes. Also depends on how powerful it is as a CAS.

For example, I'd eventually like to write a macro that lets people specify a statistical model. I could see Reduce being helpful, potentially simplifying user-provided expressions, or even marginalizing out nuisance parameters entirely.

That'd be more accurate than doing it numerically, but it'd be complicated to implement and would have to be pretty efficient to actually save time (don't want things to take an eternity to compile).

"there is still a lot of functionality that can be built on top of it. For now, I just wanted to make sure that the basic functionality is robust, but more features will come."

What sort of things do you have in mind?

Yes, the compile-time is currently slower than I would prefer, but once the parsers are compiled the performance is quicker. This is because I had to add in features to the parser to make it more robust for handling more complicated nested expressions, functions, begin/end blocks (so the code for it is large). I did this at the expense of initial compilation time. It's possible to improve the performance of it, but it requires a lot of concentration since it grew to hundreds of lines of code. I have some ideas for improving the performance but it will not happen over night.

You can take a look at the Reduce documentation for more information on its capabilities: http://www.reduce-algebra.com/manual/manual.html

Recently I implemented the parsing of matrices and there are some REDUCE packages that manipulate matrices, so I would like to build some more of those features into the Julia package. also, I would like to try out some of the differential equations solving capability in the future. I don't have a specific road plan laid out, but as the need for more features arises and as time passes, it is becoming more clear how it can be integrated into Julia.

The REDUCE software is actually one of the fastest computer algebra systems on the planet, unfortunately the text communcation via `Pipe`

and parsing into Julia aspect of it slows it down significantly; but it still might be faster overall than some of the other symbolic packages available for Julia. One notable difference is that it doesn't define a new "symbolic type" within Julia, instead the package works on Julia expression objects.

Julia does not have any fully native computer algebra packages, so this can help work around that. I wouldn't recommend using the package for doing things that you could otherwise do natively with Julia; but it should be a good tool for those situations where an external symbolic computation is essential for code generation.

```
julia> rcall(:(int( x^2 * exp(-x^2/2 ) , x, -Inf, Inf)))
:(sqrt(π) * sqrt(2))
julia> int(:( x^2 * exp(-x^2/2 ) ), :x, -Inf, Inf)
:((sqrt(π) * e ^ (inf ^ 2 / 2) * sqrt(2) * erf(inf / sqrt(2)) - 2inf) / e ^ (inf ^ 2 / 2))
```

Seems to be a bug. Besides the second not being what I'm looking for, "inf" isn't defined, so it also throws errors.

Although, most people may actually prefer errors over Infs, so they know asap something didn't go quite right.

Also, while Reduce itself doesn't need SpecialFunctions.jl, a lot of people may want it. Not sure for how long the helpful error will stick around, so may be worth at least mentioning it somewhere in the docs.

Was hoping going back to rcall would fix this:

```
julia> ex2 = :(exp(-a*exp(-x)-b*x)*x^2)
:(exp(-a * exp(-x) - b * x) * x ^ 2)
julia> int(ex2,:x,-Inf,Inf)
:(((((((((((e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=inf, int(x ^ 2 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b ^ 2 + 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=inf, int(x / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b + 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=inf, int(1 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a) - e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=-inf, int(x ^ 2 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b ^ 2) - 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=-inf, int(x / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b) - 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=-inf, int(1 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a) + e ^ ((2 * e ^ inf * b * inf + a) / e ^ inf) * b ^ 2 * inf ^ 2) - 2 * e ^ ((2 * e ^ inf * b * inf + a) / e ^ inf) * b * inf) + 2 * e ^ ((2 * e ^ inf * b * inf + a) / e ^ inf)) - e ^ (e ^ inf * a) * b ^ 2 * inf ^ 2) - 2 * e ^ (e ^ inf * a) * b * inf) - 2 * e ^ (e ^ inf * a)) / (e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * b ^ 3))
```

But:

```
julia> rcall(:(int(exp(-a*exp(-x)-b*x)*x^2),x,-Inf,Inf))
ERROR: Reduce: Nested :tuple block structure not supported
Stacktrace:
[1] show_expr(::Base.AbstractIOBuffer{Array{UInt8,1}}, ::Expr) at /home/chris/.julia/v0.6/Reduce/src/parser.jl:505
[2] unparse(::Expr) at /home/chris/.julia/v0.6/Reduce/src/parser.jl:547
[3] Type at /home/chris/.julia/v0.6/Reduce/src/rexpr.jl:38 [inlined]
[4] #rcall#14(::Array{Symbol,1}, ::Array{Symbol,1}, ::Function, ::Expr) at /home/chris/.julia/v0.6/Reduce/src/rexpr.jl:292
[5] rcall(::Expr) at /home/chris/.julia/v0.6/Reduce/src/rexpr.jl:292
julia> rcall(:(int($ex2,x,-Inf,Inf) ))
:(int(exp(-a * exp(-x) - b * x) * x ^ 2, x, -inf, inf))
```

The answer to that involves the di- and tri-gamma functions, and requires `a > 0 && b > 0`

.

I'd have to learn more about Reduce to know what is all possible. But I also definitely like the idea of it being able to work on matrices.

"Julia does not have any fully native computer algebra packages, so this can help work around that. I wouldn't recommend using the package for doing things that you could otherwise do natively with Julia; but it should be a good tool for those situations where an external symbolic computation is essential for code generation."

That's probably true for generated functions in general. Your code would probably compile faster if instead you wrote out explicit versions for all of them.

Maybe if you're expecting some numbers to be especially popular, it would be best to explicitly write those versions (possibly using Reduce or some other CAS to derive the expression for you , and then copy + paste -- I know Chris Rackauckas has often used Mathematica for that)

so the bulk of the work doesn't have to be redone upon every compilation, and then make a generated version as a fallback for rarer cases.

```
julia> int(:( x^2 * exp(-x^2/2 ) ), :x, "-infinity", "infinity")
:(sqrt(π) * sqrt(2))
julia> rcall(:(int( x^2 * exp(-x^2/2 ) , x, -Inf, Inf)))
:(sqrt(π) * sqrt(2))
```

Based on reading the REDUCE manual, it would seem that the di-gamma and tri-gamma functions are not supported, you'd have to communicate with the upstream REDUCE developers about that.

As for the `Nested :tuple block structure not supported`

error, it seems you made a mistake with your parenthesis, it should have been entered like this:

```
julia> rcall(:(int(exp(-a*exp(-x)-b*x)*x^2,x,-Inf,Inf)))
:(int(exp(-a * exp(-x) - b * x) * x ^ 2, x, -inf, inf))
julia> int(ex2,:x,"-infinity","infinity")
:(int(exp(-a * exp(-x) - b * x) * x ^ 2, x, -inf, inf))
```

which also appears to be an unsupported integral (from the upstream REDUCE end) at first glance

when an integral cannot be evaluated by REDUCE, it just returns the expression back

However, the REDUCE documentation does describe how to set the positivity of variables:

http://www.reduce-algebra.com/manual/manualse68.html#x83-15900010.2

On the di- and tri-gamma functions, I see your Julia docs already mention support for polygamma. `Polygamma(0,x) == digamma(x)`

and `polygamma(1,x) == trigramma(x)`

.

So the answer to the likely unsupported integral is: `a^-b*Gamma(b)*( (log(a)-polygamma(0,b))^2 + polygamma(1,b) )`

.

Not sure how many special functions like the bessel, hypergeometrics, etc, it does support. Not sure how likely Reduce developers are to be able to increase support for things like that.

Likely unsupported, because it requires a > 0 && b > 0. Perhaps it could return an answer, given those conditions. From the Wolfram sandbox:

https://sandbox.open.wolframcloud.com/app/objects/384d34c1-ac46-46a5-ad28-b91ed2876d5e#sidebar=compute`Assuming[a > 0 && b > 0, Integrate[Exp[-(a Exp[-x]) - b x] x^2, {x, -Infinity, Infinity}]]`

Is where I got the answer. If I leave off the constraint, it'll return the same answer within a conditional expression.

Not sure how well Reduce handles conditionals; I'll get around to looking through the docs more thoroughly as I try to use it more. For now, http://www.reduce-algebra.com/manual/manualse16.html#x25-380005.3 appears to only cover user-provided conditionals.

Also, I'm glad I can pass reduce strings as a fall back. Trying to follow the example on the conditionals:

```
julia> rcall( :(abs(a*b*c)), "let sign(a) => 1, sign(b) => 1" )
ERROR: Reduce:
***** nil not defined as switch
abs(c)*a*b
***** nil not defined as switch
Stacktrace:
[1] ReduceCheck(::String) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/Reduce.jl:62
[2] read(::Reduce.PSL) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/Reduce.jl:84
[3] readsp(::Reduce.PSL) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/Reduce.jl:88
[4] #rcall#13(::Array{Symbol,1}, ::Array{Symbol,1}, ::Function, ::Reduce.RExpr) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/rexpr.jl:245
[5] (::Reduce.#kw##rcall)(::Array{Any,1}, ::Reduce.#rcall, ::Reduce.RExpr) at ./<missing>:0
[6] #rcall#14(::Array{Symbol,1}, ::Array{Symbol,1}, ::Function, ::Expr) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/rexpr.jl:278
[7] (::Reduce.#kw##rcall)(::Array{Any,1}, ::Reduce.#rcall, ::Expr) at ./<missing>:0
[8] rcall(::Expr, ::String, ::Vararg{String,N} where N) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/rexpr.jl:282
```

So, I see Reduce itself returned the right answer, but the parser failed.

You are using

`rcall`

incorrectly there. Extra arguments to `rcall`

have to be REDUCE switch names passed as Julia symbols, you can't just send it any random string. That is why the REDUCE error tells you that you did not enter a defined switch. For `int`

and `df`

it's different, since those can accept extra args. You have to make 2 separate calls to `rcall`

if you want to first evaluate a REDUCE string and then a Julia expression object like this:```
julia> rcall("let sign(a) => 1, sign(b) => 1")
""
julia> rcall( :(abs(a*b*c)))
:(abs(c) * a * b)
```

```
julia> solve(((:m - :x) ^ 2 - :s) / :s ^ 2, :s)
1-element Array{Any,1}:
:(s = (m ^ 2 - 2 * m * x) + x ^ 2)
```