How to use Base.rationalize on whole symbolic expression?

New to Julia
,
1

fyi, issue opened

I obtain a symbolic expression from call to one of Julia symbolic functions (integrate) that has a mix of floating points and integers.

I’d like to use Base.rationalize to rationalize these floating points numbers within the expression to becomes numer//denom instead of number.decimal format.

I can’t do Base.rationalize(expression) as this gives an error.

Is there a way to map Base.rationalize only on the real numbers within the expression? Here is a full example.

using Symbolics
@syms x
expr=(0.014378665475956662 - 0.0645731363698366im)*((72 + 132x + 2262(x^4) + 6743(x^5) - 754(x^2) - 234(x^6) - 1641(x^3) - 9180(x^7) - 5400(x^8))^-1) + (-0.0995313498923154 + 0.07922575111068406im)*x*((72 + 132x + 2262(x^4) + 6743(x^5) - 754(x^2) - 234(x^6) - 1641(x^3) - 9180(x^7) - 5400(x^8))^-1)

I could do

julia> Base.rationalize(0.014378665475956662 - 0.0645731363698366im)
26488241//1842190504 - 36390080//563548281*im

But I am doing this in program, and would like to just do

julia> Base.rationalize(expr)

ERROR: MethodError: no method matching rationalize(::SymbolicUtils.Add{Number, Int64, Dict{Any, Number}, Nothing})
Closest candidates are:
  rationalize(::Type{T}, ::AbstractFloat; tol) where T<:Integer at rational.jl:217
  rationalize(::Type{T}, ::AbstractFloat, ::Real) where T<:Integer at rational.jl:157
  rationalize(::AbstractFloat; kvs...) at rational.jl:218
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[16]:1

Any suggestion how to do this? In Mathematica, this is possible to do on the whole expression. Here is an example in Mathematica

expr = (0.014378665475956662 - 0.0645731363698366*I)/(72 + 132*x - 754*x^2 - 1641*x^3 + 2262*x^4 + 6743*x^5 - 234*x^6 - 9180*x^7 - 5400*x^8) - 
  ((0.0995313498923154 - 0.07922575111068406*I)*x)/(72 + 132*x - 754*x^2 - 1641*x^3 + 2262*x^4 + 6743*x^5 - 234*x^6 - 9180*x^7 - 5400*x^8)

SetAccuracy[expr, Infinity]

(1036092039673655/72057594037927936 - (4652984846293445*I)/72057594037927936)/(72 + 132*x - 754*x^2 - 1641*x^3 + 2262*x^4 + 6743*x^5 - 234*x^6 - 9180*x^7 - 5400*x^8) - 
  ((7171989604587425/72057594037927936 - (2854408505441795*I)/36028797018963968)*x)/(72 + 132*x - 754*x^2 - 1641*x^3 + 2262*x^4 + 6743*x^5 - 234*x^6 - 9180*x^7 - 5400*x^8)

Can this be done in Julia? I am using 1.8

1 Like
2

Open an issue.

3

Chris’ answer implies that it should work.

But… I am not in any way an expert on this but Symbolics always fascinates me, so I tried to find a way to do this and I thought it should work somehow recursively. This is the solution I created. It’s not complete and generic enough to work on all expressions, but it can easily be expanded (so I hope).

And maybe someone is motivated to show us a more Symbolics-Way of doing it.

using Symbolics

sr = @rule ( +(~~x) ) => ~~x
mr = @rule ( *(~~x) ) => ~~x
pr = @rule ( ^(~x,~y) ) => [~x,~y]

function rationalize!(a)
	for index in eachindex(a)
		if typeof(a[index]) <: AbstractFloat || typeof(a[index]) <: ComplexF64
			a[index] = Base.rationalize(a[index])
		end
	end
	a
end

function decompose_compose(expr)
	if typeof(expr) <: SymbolicUtils.Sym || typeof(expr) <: Number
		return expr
	end
	if ! isnothing(sr(expr))
		subexpr = sr(expr)
		compose=sum
	elseif ! isnothing(mr(expr))
		subexpr = mr(expr)
		compose=prod
	elseif ! isnothing(pr(expr))
		subexpr = pr(expr)
		compose=^
	else
		display(expr)
		error("no rule defined")
	end
	subexpr = rationalize!(subexpr)
	subexpr = decompose_compose.(subexpr)
	(compose == ^) ? compose(subexpr...) : compose(subexpr)
end

@syms x

expr=(0.014378665475956662 - 0.0645731363698366im)*((72 + 132x + 2262(x^4) + 6743(x^5) - 754(x^2) - 234(x^6) - 1641(x^3) - 9180(x^7) - 5400(x^8))^-1) + (-0.0995313498923154 + 0.07922575111068406im)*x*((72 + 132x + 2262(x^4) + 6743(x^5) - 754(x^2) - 234(x^6) - 1641(x^3) - 9180(x^7) - 5400(x^8))^-1)

decompose_compose(expr)
1 Like
4

Thanks. This seems to work OK. It is ok if it fails on some expressions. Will add try/catch to catch the error in those cases.

Here is one case where it can’t handle

julia> expr=0.8103075638223697(((3//5) + x)^-1) + 0.0012552102783267068(((1//4) + x^2 - x)^-1) + 1.3505126063706157log((3//5) + x) + 0.01769427781144674(x^2)*(((1//4) + x^2 - x)^-1) - 0.804168197111734(((9//25) + (6//5)*x + x^2)^-1) - 0.06487901864197139log(x - (1//2)) - 0.011357559462376793x*(((1//4) + x^2 - x)^-1) - 1.3402803285195568x*(((9//25) + (6//5)*x + x^2)^-1)

julia> decompose_compose(expr)

gives

log(x - (1//2))
ERROR: no rule defined
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] decompose_compose(expr::SymbolicUtils.Term{Number, Nothing})
    @ Main .\REPL[6]:16
  [3] _broadcast_getindex_evalf
    @ .\broadcast.jl:670 [inlined]
  [4] _broadcast_getindex
    @ .\broadcast.jl:643 [inlined]
  [5] getindex
    @ .\broadcast.jl:597 [inlined]
  [6] copyto_nonleaf!(dest::Vector{Rational{Int64}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(decompose_compose), Tuple{Base.Broadcast.Extruded{SubArray{Any, 1, Vector{Any}, Tuple{UnitRange{Int64}}, true}, Tuple{Bool}, Tuple{Int64}}}}, iter::Base.OneTo{Int64}, state::Int64, count::Int64)
    @ Base.Broadcast .\broadcast.jl:1055
  [7] copy
    @ .\broadcast.jl:907 [inlined]
  [8] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(decompose_compose), Tuple{SubArray{Any, 1, Vector{Any}, Tuple{UnitRange{Int64}}, true}}})
    @ Base.Broadcast .\broadcast.jl:860
  [9] decompose_compose(expr::SymbolicUtils.Mul{Number, Float64, Dict{Any, Number}, Nothing})
    @ Main .\REPL[6]:19
 [10] _broadcast_getindex_evalf
    @ .\broadcast.jl:670 [inlined]
 [11] _broadcast_getindex
    @ .\broadcast.jl:643 [inlined]
 [12] getindex
    @ .\broadcast.jl:597 [inlined]
 [13] copyto_nonleaf!(dest::Vector{SymbolicUtils.Mul{Number, Rational{Int64}, Dict{Any, Number}, Nothing}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(decompose_compose), Tuple{Base.Broadcast.Extruded{SubArray{SymbolicUtils.Mul{Number, Float64, Dict{Any, Number}, Nothing}, 1, Vector{SymbolicUtils.Mul{Number, Float64, Dict{Any, Number}, Nothing}}, Tuple{UnitRange{Int64}}, true}, Tuple{Bool}, Tuple{Int64}}}}, iter::Base.OneTo{Int64}, state::Int64, count::Int64)
    @ Base.Broadcast .\broadcast.jl:1055
 [14] copy
    @ .\broadcast.jl:907 [inlined]
 [15] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(decompose_compose), Tuple{SubArray{SymbolicUtils.Mul{Number, Float64, Dict{Any, Number}, Nothing}, 1, Vector{SymbolicUtils.Mul{Number, Float64, Dict{Any, Number}, Nothing}}, Tuple{UnitRange{Int64}}, true}}})
    @ Base.Broadcast .\broadcast.jl:860
 [16] decompose_compose(expr::SymbolicUtils.Add{Number, Int64, Dict{Any, Number}, Nothing})
    @ Main .\REPL[6]:19
 [17] top-level scope
    @ REPL[9]:1

julia>

Is it possible to make your function not fail when there is no real numbers in the expression? Here is an example

expr=(1//2)*log(1 + x) - (1//4)*log(1 + x^2) - (1//2)*atan(-x)
decompose_compose(expr)

gives

log(1 + x^2)
ERROR: no rule defined
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] decompose_compose(expr::SymbolicUtils.Term{Number, Nothing})
    @ Main .\REPL[8]:16
  [3] _broadcast_getindex_evalf
    @ .\broadcast.jl:670 [inlined]
  [4] _broadcast_getindex
    @ .\broadcast.jl:643 [inlined]
  [5] getindex
    @ .\broadcast.jl:597 [inlined]
  [6] copyto_nonleaf!(dest::Vector{Rational{Int64}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(decompose_compose), Tuple{Base.Broadcast.Extruded{SubArray{Any, 1, Vector{Any}, Tuple{UnitRange{Int64}}, true}, Tuple{Bool}, Tuple{Int64}}}}, iter::Base.OneTo{Int64}, state::Int64, count::Int64)
    @ Base.Broadcast .\broadcast.jl:1055
  [7] copy
    @ .\broadcast.jl:907 [inlined]
  [8] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(decompose_compose), Tuple{SubArray{Any, 1, Vector{Any}, Tuple{UnitRange{Int64}}, true}}})
    @ Base.Broadcast .\broadcast.jl:860
  [9] decompose_compose(expr::SymbolicUtils.Mul{Number, Rational{Int64}, Dict{Any, Number}, Nothing})
    @ Main .\REPL[8]:19
 [10] _broadcast_getindex_evalf
    @ .\broadcast.jl:670 [inlined]
 [11] _broadcast_getindex
    @ .\broadcast.jl:643 [inlined]
 [12] getindex
    @ .\broadcast.jl:597 [inlined]
 [13] copy
    @ .\broadcast.jl:899 [inlined]
 [14] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(decompose_compose), Tuple{SubArray{SymbolicUtils.Mul{Number, Rational{Int64}, Dict{Any, Number}, Nothing}, 1, Vector{SymbolicUtils.Mul{Number, Rational{Int64}, Dict{Any, Number}, Nothing}}, Tuple{UnitRange{Int64}}, true}}})
    @ Base.Broadcast .\broadcast.jl:860
 [15] decompose_compose(expr::SymbolicUtils.Add{Number, Int64, Dict{Any, Number}, Nothing})
    @ Main .\REPL[8]:19
 [16] top-level scope
    @ REPL[9]:1

Thanks

5

Just add two more special rules for log and atan.

But this shows clearly that this solution isn’t very scaleable. Perhaps this is what Chris had in mind, that opening an issue would help to keep track, what people need and make it possible in future versions.

So I repeat: please open an issue here Issues · JuliaSymbolics/SymbolicUtils.jl · GitHub , it’s easy and you can refer to this discussion.

Adding two rules:

using Symbolics

sr = @rule ( +(~~x) ) => ~~x
mr = @rule ( *(~~x) ) => ~~x
pr = @rule ( ^(~x,~y) ) => [~x,~y]
logr = @rule ( log(~x) ) => [~x]
atanr = @rule ( atan(~x) ) => [~x]

function rationalize!(a)
	for index in eachindex(a)
		if typeof(a[index]) <: AbstractFloat || typeof(a[index]) <: ComplexF64
			a[index] = Base.rationalize(a[index])
		end
	end
	a
end

function decompose_compose(expr)
	if typeof(expr) <: SymbolicUtils.Sym || typeof(expr) <: Number
		return expr
	end
	if ! isnothing(sr(expr))
		subexpr = sr(expr)
		compose=sum
	elseif ! isnothing(mr(expr))
		subexpr = mr(expr)
		compose=prod
	elseif ! isnothing(pr(expr))
		subexpr = pr(expr)
		compose=^
	elseif ! isnothing(logr(expr))
		subexpr = logr(expr)
		compose=log
	elseif ! isnothing(atanr(expr))
		subexpr = atanr(expr)
		compose=log
	else
		display(expr)
		error("no rule defined")
	end
	subexpr = rationalize!(subexpr)
	subexpr = decompose_compose.(subexpr)
	(compose == ^) || (compose == log) || (compose == atan) ? compose(subexpr...) : compose(subexpr)
end

@syms x

expr=(1//2)*log(1 + x) - (1//4)*log(1 + x^2) - (1//2)*atan(-x)

decompose_compose(expr)
6

But I already did, few hrs ago. You can see the link at the top of my question above.

Thanks for the update also.

1 Like
7

Ah, great. Sorry, I missed it and didn’t checked.

8

Maybe something like this, which works with the tree structure returned by SymbolicUtils, will work for you.

_rat(x::Num) = _rat(Symbolics.value(x))

function _rat(x::SymbolicUtils.Add)
    +(_rat.(arguments(x))...)
end

function _rat(x::SymbolicUtils.Mul)
    *(_rat.(arguments(x))...)
end

function _rat(x::SymbolicUtils.Div)
    /(_rat.(arguments(x))...)
end

function _rat(x::SymbolicUtils.Pow)
    ^(_rat.(arguments(x))...)
end

function _rat(x::SymbolicUtils.Term)
    op = operation(x)
    op(_rat.(arguments(x))...)
end

function _rat(x::SymbolicUtils.Sym)
    x
end

_rat(x::Integer) = x
_rat(x::Rational) = x
_rat(x) = rationalize(x)
3 Likes