How to use Base.rationalize on whole symbolic expression?

Nasser · September 4, 2022, 6:58am

fyi, issue opened

Could `Base.rationalize` be made to work on symbolic expression?

opened 09:21PM - 04 Sep 22 UTC

reference https://discourse.julialang.org/t/how-to-use-base-rationalize-on-whole…-symbolic-expression/86748 Using 1.8 `Base.rationalize(expr)` gives an error on symbolic expression ``` 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) Base.rationalize(expr) ``` gives ``` 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 ``` Could `Base.rationalize` be made to work on whole symbolic expression?

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

ChrisRackauckas · September 4, 2022, 6:46pm

Open an issue.

oheil · September 4, 2022, 6:56pm

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)

Nasser · September 4, 2022, 9:23pm

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

oheil · September 5, 2022, 7:30am

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)

Nasser · September 5, 2022, 7:49am

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

Thanks for the update also.

oheil · September 5, 2022, 7:57am

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

j_verzani · September 5, 2022, 6:32pm

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)

Topic		Replies	Views
Simplify Julia Expr with Symbolics.jl General Usage question , symbolics	2	269	February 25, 2024
Symbolics.jl expressions that don't reduce to single type Modelling & Simulations question , symbolics	5	367	September 3, 2023
Expression as single fraction in Symbolics.jl General Usage symbolics	4	1216	September 28, 2021
Imaginary unit in Julia Symbolics General Usage question	1	294	January 17, 2023
Simplifying Symbolics.jl expression General Usage symbolics	2	118	September 3, 2024

How to use Base.rationalize on whole symbolic expression?

Related topics