## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
• Jan 29 2019 01:29
chakravala labeled #21
• Jan 29 2019 01:27
chakravala closed #10
• Jan 29 2019 01:27
chakravala closed #13
• Jan 29 2019 01:25
chakravala commented #13
• Jan 29 2019 01:22
chakravala labeled #23
• Jan 28 2019 17:39
flipgthb starred chakravala/Reduce.jl
• Jan 27 2019 15:59
chakravala commented #23
• Jan 27 2019 13:50
atthom commented #23
• Jan 26 2019 16:30
chakravala commented #23
• Jan 26 2019 01:39
atthom opened #23
• Jan 22 2019 11:27
chakravala labeled #22
• Jan 22 2019 09:35
chakravala commented #22
• Jan 22 2019 09:34
chakravala closed #22
• Jan 22 2019 09:34
chakravala commented #22
• Jan 22 2019 02:12
cscherrer commented #22
• Jan 22 2019 00:02
chakravala commented #22
• Jan 21 2019 21:08
cscherrer opened #22
• Jan 21 2019 19:53
cscherrer commented #21
• Jan 21 2019 19:29
chakravala commented #21
• Jan 21 2019 19:23
chakravala commented #21
Dream Scatter
@chakravala
Hi, let me know how you are using Reduce.jl or if a feature you need is missing
chriselrod
@chriselrod

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
Dream Scatter
@chakravala
Cool, thanks for sharing. Yea, this package is also a proof of concept at the moment, 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.
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).
chriselrod
@chriselrod

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?

Dream Scatter
@chakravala

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.

chriselrod
@chriselrod
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.

chriselrod
@chriselrod

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. chriselrod @chriselrod "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. Dream Scatter @chakravala Thanks for your feedback. I will look into those errors, it should be possible to handle the infinities. Dream Scatter @chakravala To correctly calculate the infinities, use int(:( x^2 * exp(-x^2/2 ) ), :x, "-infinity", "infinity") because the parser currently does not recognize the Inf Julia type at the moment, but REDUCE strings work. 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)) Dream Scatter @chakravala 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. Dream Scatter @chakravala 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 Dream Scatter @chakravala However, the REDUCE documentation does describe how to set the positivity of variables: http://www.reduce-algebra.com/manual/manualse68.html#x83-15900010.2 Dream Scatter @chakravala The newest master branch of Reduce.jl parser now supports the Julia Inf and -Inf values of Float64: julia> int(:( x^2 * exp(-x^2/2 ) ), :x, -Inf, Inf) :(sqrt(π) * sqrt(2)) chriselrod @chriselrod 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. Dream Scatter @chakravala 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) chriselrod @chriselrod Okay, cool, thanks Dream Scatter @chakravala The only extra arguments that rcall takes are REDUCE switch names that are defined, for example: julia> rcall(:(x^2+2x+1),:factor) :((x + 1) ^ 2) julia> factor(:(x^2+2x+1)) :((x + 1) ^ 2) But you can also call REDUCE switches with the syntax in multiple ways chriselrod @chriselrod Unfortunately, still get the same result for that integral. Also, perhaps it would be worth adding the "solve" switch? But I can still use it through rcall: julia> rcall(:(solve(((m - x) ^ 2 - s) / s ^ 2, s))) :($(Expr(:cell1d, :(s = (m ^ 2 - 2 * m * x) + x ^ 2))))
Dream Scatter
@chakravala
Thanks for sending me all these examples, the solve and cell syntax is not handled yet, but I can look into that. At the moment I don't have time available to implement more features, but this weekend I might.
Dream Scatter
@chakravala
Note that solve is not a switch. Please take a look at the REDUCE manual to understand what a switch is. Also, the REDUCE developers are worth contacting on SourceForge, I can't speak for what they will or will not implement.
chriselrod
@chriselrod
No worries -- thanks for addressing my concerns. I'll look at the manual more before bugging you further.