[BREAKING] feat: update to SymbolicUtils@4

Project.toml

@@ -1,7 +1,7 @@
 name = "Symbolics"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 authors = ["Shashi Gowda <gowda@mit.edu>"]
-version = "6.56.0"
+version = "7.0.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -11,25 +11,24 @@ CommonWorldInvalidations = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
-LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+Moshi = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d"
+MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
-SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLPublic = "431bcebd-1456-4ced-9d72-93c2757fff0b"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -42,20 +41,27 @@ TermInterface = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
 
 [weakdeps]
 D3Trees = "e3df1716-f71e-5df9-9e2d-98e193103c45"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
+PlutoRunner = "dc6b355a-2368-4481-ae6d-ae0351418d79"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
 
 [extensions]
 SymbolicsD3TreesExt = "D3Trees"
+SymbolicsDistributionsExt = "Distributions"
 SymbolicsForwardDiffExt = "ForwardDiff"
 SymbolicsGroebnerExt = "Groebner"
+SymbolicsLatexifyExt = ["Latexify", "LaTeXStrings"]
 SymbolicsLuxExt = "Lux"
 SymbolicsNemoExt = "Nemo"
+SymbolicsPlutoRunnerExt = "PlutoRunner"
 SymbolicsPreallocationToolsExt = ["PreallocationTools", "ForwardDiff"]
 SymbolicsSymPyExt = "SymPy"
 SymbolicsSymPyPythonCallExt = "SymPyPythonCall"
@@ -84,19 +90,21 @@ Latexify = "0.16"
 LogExpFunctions = "0.3"
 Lux = "1"
 MacroTools = "0.5"
+Moshi = "0.3.7"
+MultivariatePolynomials = "0.5.12"
+MutableArithmetics = "1.6.5"
 NaNMath = "1"
 Nemo = "0.46, 0.47, 0.48, 0.49, 0.52"
-OffsetArrays = "1.15.0"
 PkgBenchmark = "0.2"
 PreallocationTools = "0.4"
 PrecompileTools = "1"
+Preferences = "1.5.0"
 Primes = "0.5"
 RecipesBase = "1.1"
 Reexport = "1"
 ReferenceTests = "0.10"
 RuntimeGeneratedFunctions = "0.5.9"
 SafeTestsets = "0.1.0"
-SciMLBase = "2"
 SciMLPublic = "1.0.0"
 Setfield = "1"
 SnoopCompile = "3.1"
@@ -106,7 +114,7 @@ SymPy = "2.2"
 SymPyPythonCall = "0.5"
 SymbolicIndexingInterface = "0.3.14"
 SymbolicLimits = "0.2.2"
-SymbolicUtils = "3.24"
+SymbolicUtils = "4.2"
 TermInterface = "2"
 julia = "1.10"
 
@@ -117,7 +125,9 @@ D3Trees = "e3df1716-f71e-5df9-9e2d-98e193103c45"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LambertW = "984bce1d-4616-540c-a9ee-88d1112d94c9"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
@@ -126,8 +136,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SnoopCompile = "aa65fe97-06da-5843-b5b1-d5d13cad87d2"
-SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "Groebner", "BenchmarkTools", "ReferenceTests", "SymPyPythonCall", "Random", "LambertW", "Lux", "ComponentArrays", "Nemo", "DynamicQuantities", "D3Trees", "SnoopCompile"]
+test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "Groebner", "BenchmarkTools", "ReferenceTests", "Random", "LambertW", "Lux", "ComponentArrays", "Nemo", "DynamicQuantities", "D3Trees", "SnoopCompile", "Latexify", "LaTeXStrings"]

docs/Project.toml

@@ -19,5 +19,5 @@ Nemo = "0.46, 0.47, 0.48, 0.49, 0.52"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
 StaticArrays = "1.5"
-SymbolicUtils = "2.0.2, 3"
-Symbolics = "6"
+SymbolicUtils = "4.2"
+Symbolics = "7"

ext/SymbolicsDistributionsExt.jl

@@ -0,0 +1,37 @@
+module SymbolicsDistributionsExt
+
+using Symbolics
+using Distributions
+using Symbolics: Num, Arr, VartypeT, unwrap, BasicSymbolic, shape
+
+
+for f in [pdf, logpdf, cdf, logcdf, quantile]
+    @eval function (::$(typeof(f)))(dist::Distributions.Distribution, x::Num)
+        $f(dist, unwrap(x))
+    end
+    @eval function (::$(typeof(f)))(dist::Distributions.Distribution, x::Arr)
+        $f(dist, unwrap(x))
+    end
+    @eval function (::$(typeof(f)))(dist::BasicSymbolic{VartypeT}, x::Num)
+        $f(dist, unwrap(x))
+    end
+    @eval function (::$(typeof(f)))(dist::BasicSymbolic{VartypeT}, x::Arr)
+        $f(dist, unwrap(x))
+    end
+    @eval function (::$(typeof(f)))(dist::BasicSymbolic{VartypeT}, x)
+        $f(dist, unwrap(x))
+    end
+end
+
+for f in [Distributions.Uniform, Distributions.Normal]
+    for (T1, T2) in Iterators.product(Iterators.repeated([Real, BasicSymbolic{VartypeT}, Num], 2)...)
+        if T1 != Num && T2 != Num
+            continue
+        end
+        @eval function (::Type{$f})(a::$T1, b::$T2)
+            $f(unwrap(a), unwrap(b))
+        end
+    end
+end
+
+end

ext/SymbolicsGroebnerExt.jl

@@ -5,19 +5,20 @@ const Nemo = Groebner.Nemo
 using Symbolics
 using Symbolics: Num, symtype, BasicSymbolic
 import Symbolics.PrecompileTools
+import Symbolics.Bijections
+import Symbolics: SymbolicUtils
+import Symbolics: DP, MP
 
 function Symbolics.groebner_basis(polynomials::Vector{Num}; ordering=InputOrdering(), kwargs...)
-    polynoms, pvar2sym, sym2term = Symbolics.symbol_to_poly(polynomials)
-    sym2term_for_groebner = Dict{Any,Any}(v1 => k for (k, (v1, v2)) in sym2term)
-    all_sym_vars = Groebner.ordering_variables(ordering)
-    missed = setdiff(all_sym_vars, Set(collect(keys(sym2term_for_groebner))))
-    for var in missed
-        sym2term_for_groebner[var] = var
-    end
-    ordering = Groebner.ordering_transform(ordering, sym2term_for_groebner )
-    basis = Groebner.groebner(polynoms; ordering=ordering, kwargs...)
-    PolyType = symtype(first(polynomials))
-    Symbolics.poly_to_symbol(basis, pvar2sym, sym2term, PolyType)
+    polynoms, poly_to_bs = Symbolics.symbol_to_poly(polynomials)
+    basis = groebner_basis_poly(polynoms, poly_to_bs; ordering, kwargs...)
+    Symbolics.poly_to_symbol(basis, poly_to_bs)
+end
+
+function groebner_basis_poly(polynoms::Vector{<:DP.Polynomial}, poly_to_bs::Bijections.Bijection; ordering=InputOrdering(), kwargs...)
+    bs_to_poly = Bijections.active_inv(poly_to_bs)
+    ordering = Groebner.ordering_transform(ordering, bs_to_poly)
+    return Groebner.groebner(polynoms; ordering=ordering, kwargs...)
 end
 
 """
@@ -41,8 +42,8 @@ julia> @variables x y;
 julia> is_groebner_basis([x^2 - y^2, x*y^2 + x, y^3 + y])
 ```
 """
-function Symbolics.is_groebner_basis(polynomials::Vector{<:Union{Num, BasicSymbolic{<:Number}}}; kwargs...)
-    polynoms, _, _ = Symbolics.symbol_to_poly(polynomials)
+function Symbolics.is_groebner_basis(polynomials::Vector{<:Union{Num, BasicSymbolic}}; kwargs...)
+    polynoms, _ = Symbolics.symbol_to_poly(polynomials)
     Groebner.isgroebner(polynoms; kwargs...)
 end
 
@@ -66,7 +67,13 @@ function nemo_crude_evaluate(poly::Nemo.MPolyRingElem, varmap)
 end
 
 function nemo_crude_evaluate(poly::Nemo.FracElem, varmap)
-    nemo_crude_evaluate(numerator(poly), varmap) // nemo_crude_evaluate(denominator(poly), varmap)
+    num = nemo_crude_evaluate(numerator(poly), varmap)
+    den = nemo_crude_evaluate(denominator(poly), varmap)
+    if num isa Num || den isa Num
+        num / den
+    else
+        num // den
+    end
 end
 
 function nemo_crude_evaluate(poly::Nemo.ZZRingElem, varmap)
@@ -86,27 +93,27 @@ function gen_separating_var(vars)
 end
 
 # Given a GB in k[params][vars] produces a GB in k(params)[vars]
-function demote(gb, vars::Vector{Num}, params::Vector{Num})
-    isequal(gb, [1]) && return gb 
-
-    gb = Symbolics.wrap.(SymbolicUtils.toterm.(gb))
+function demote(gb, gb_as_poly, bs_to_poly, vars::Vector{Num}, params::Vector{Num})
+    length(gb) == 1 && SymbolicUtils._isone(gb[1]) && return gb
+    isequal(gb, [1]) && return gb
+    # gb = Symbolics.wrap.(SymbolicUtils.toterm.(gb))
     Symbolics.check_polynomial.(gb)
 
     all_vars = [vars..., params...]
     nemo_ring, nemo_all_vars = Nemo.polynomial_ring(Nemo.QQ, map(string, all_vars))
 
     sym_to_nemo = Dict(all_vars .=> nemo_all_vars)
     nemo_to_sym = Dict(v => k for (k, v) in sym_to_nemo)
-    nemo_gb = Symbolics.substitute(gb, sym_to_nemo)
-    nemo_gb = Symbolics.substitute(nemo_gb, sym_to_nemo)
+    pvar_to_nemo = [bs_to_poly[v] for v in all_vars] => nemo_all_vars
+    nemo_gb = [poly(pvar_to_nemo) for poly in gb_as_poly]
 
-    nemo_vars = filter(v -> string(v) in string.(vars), nemo_all_vars)
-    nemo_params = filter(v -> string(v) in string.(params), nemo_all_vars)
+    nemo_vars = view(nemo_all_vars, 1:length(vars))
+    nemo_params = view(nemo_all_vars, length(vars)+1:length(all_vars))
 
     ring_flat = parent(nemo_vars[1])
     ring_param, params_demoted = Nemo.polynomial_ring(Nemo.base_ring(ring_flat), map(string, nemo_params))
     ring_demoted, vars_demoted = Nemo.polynomial_ring(Nemo.fraction_field(ring_param), map(string, nemo_vars), internal_ordering=:lex)
-    varmap = Dict((nemo_vars .=> vars_demoted)..., (nemo_params .=> params_demoted)...)
+    varmap = Dict(vcat(nemo_vars .=> vars_demoted, nemo_params .=> params_demoted))
     gb_demoted = map(f -> ring_demoted(nemo_crude_evaluate(f, varmap)), nemo_gb)
     result = empty(gb_demoted)
     while true
@@ -124,7 +131,7 @@ function demote(gb, vars::Vector{Num}, params::Vector{Num})
     end
     @assert all(f -> isone(Nemo.leading_coefficient(f)), result)
 
-    sym_to_nemo = Dict(sym => nem for sym in all_vars for nem in [vars_demoted..., params_demoted...] if isequal(string(sym),string(nem)))
+    sym_to_nemo = Dict(all_vars .=> [vars_demoted; params_demoted])
     nemo_to_sym = Dict(v => k for (k, v) in sym_to_nemo)
 
     final_result = Num[]
@@ -150,9 +157,8 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
     # AAECC 9, 433–461 (1999). https://doi.org/10.1007/s002000050114
 
     rng = Groebner.Random.Xoshiro(42)
-
-    all_indeterminates = reduce(union, map(Symbolics.get_variables, eqs))
-    params = map(Symbolics.Num ∘ Symbolics.wrap, setdiff(all_indeterminates, vars))
+    all_indeterminates = collect(reduce(union!, map(Symbolics.get_variables, eqs)))
+    params = map(Symbolics.Num, setdiff(all_indeterminates, vars))
 
     # Use a new variable to separate the input polynomials (Reference above)
     new_var = gen_separating_var(vars)
@@ -179,23 +185,24 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
 
         push!(new_eqs, separating_form)
 
-        new_eqs = Symbolics.groebner_basis(new_eqs, ordering=Lex(vcat(vars, params)))
+        polynoms, poly_to_bs = Symbolics.symbol_to_poly(new_eqs)
+        basis = groebner_basis_poly(polynoms, poly_to_bs; ordering = Lex(vcat(vars, params)))
+        new_eqs = Symbolics.poly_to_symbol(basis, poly_to_bs)
 
         # handle "unsolvable" case
-        if isequal(1, new_eqs[1])
+        if SymbolicUtils._iszero(new_eqs[1])
             return []
         end
 
         for i in reverse(eachindex(new_eqs))
             all_present = Symbolics.get_variables(new_eqs[i])
             if length(intersect(all_present, vars)) < 1
                 deleteat!(new_eqs, i)
+                deleteat!(basis, i)
             end
         end
-
-        new_eqs = demote(new_eqs, vars, params)
+        new_eqs = demote(new_eqs, basis, Bijections.active_inv(poly_to_bs), vars, params)
         new_eqs = map(Symbolics.unwrap, new_eqs)
-
         # condition for positive dimensionality, i.e. infinite solutions
         if length(new_eqs) < length(vars)
             warns && @warn("Infinite number of solutions")
@@ -209,13 +216,15 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
         # xn - fn(T, params) = 0
         generating = !(length(new_eqs) == length(vars))
         if length(new_eqs) == length(vars)
-            generating |= !(isequal(setdiff(Symbolics.get_variables(new_eqs[1]), params), [new_var]))
+            vars_in_1 = Symbolics.get_variables(new_eqs[1])
+            setdiff!(vars_in_1, params)
+            generating |= !(length(vars_in_1) == 1 && isequal(first(vars_in_1), new_var))
             for i in eachindex(new_eqs)[2:end]
                 present_vars = setdiff(Symbolics.get_variables(new_eqs[i]), new_var)
                 present_vars = setdiff(present_vars, params)
                 isempty(present_vars) && (generating = false; break;)
-                var_i = present_vars[1]
-                condition1 = isequal(present_vars, [var_i])
+                var_i = first(present_vars)
+                condition1 = length(present_vars) == 1
                 condition2 = Symbolics.degree(new_eqs[i], var_i) == 1
                 generating |= !(condition1 && condition2)
             end
@@ -234,7 +243,9 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
 
     # first, solve the first minimal polynomial
     @assert length(new_eqs) == length(vars)
-    @assert isequal(setdiff(Symbolics.get_variables(new_eqs[1]), params), [new_var])
+    vars_in_1 = Symbolics.get_variables(new_eqs[1])
+    setdiff!(vars_in_1, params)
+    @assert length(vars_in_1) == 1 && isequal(first(vars_in_1), new_var)
     minpoly_sols = Symbolics.symbolic_solve(Symbolics.wrap(new_eqs[1]), new_var, dropmultiplicity=dropmultiplicity)
     solutions = [Dict{Num, Any}(new_var => sol) for sol in minpoly_sols]
 
@@ -243,15 +254,15 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
     # second, iterate over eqs and sub each found solution
     # then add the roots of the remaining unknown variables 
     for (i, eq) in enumerate(new_eqs)
-        present_vars = setdiff(Symbolics.get_variables(eq), params)
-        present_vars = setdiff(present_vars, new_var)
+        present_vars = setdiff!(Symbolics.get_variables(eq), params)
+        present_vars = setdiff!(present_vars, new_var)
         @assert length(present_vars) == 1
-        var_tosolve = present_vars[1]
+        var_tosolve = first(present_vars)
         @assert Symbolics.degree(eq, var_tosolve) == 1
         @assert !isempty(solutions)
         for roots in solutions
-            subbded_eq = Symbolics.substitute(eq, Dict([new_var => roots[new_var]]); fold=false)
-            subbded_eq = Symbolics.substitute(subbded_eq, Dict([var_tosolve => 0]); fold=false)
+            subbded_eq = Symbolics.substitute(eq, Dict([new_var => roots[new_var]]); fold=Val(false))
+            subbded_eq = Symbolics.substitute(subbded_eq, Dict([var_tosolve => 0]); fold=Val(false))
             new_var_sols = [-subbded_eq]
             @assert length(new_var_sols) == 1
             root = new_var_sols[1]
@@ -270,7 +281,7 @@ end
 function transendence_basis(sys, vars)
     J = Symbolics.jacobian(sys, vars)
     x0 = Dict(v => rand(-10:10) for v in vars)
-    J_x0 = substitute(J, x0)
+    J_x0 = map(Symbolics.value, substitute(J, x0))
     rk, rref = Nemo.rref(Nemo.matrix(Nemo.QQ, J_x0))
     pivots = Int[]
     for i in 1:length(sys)
@@ -287,7 +298,7 @@ function Symbolics.solve_multivar(eqs::Vector, vars::Vector{Num}; dropmultiplici
     isempty(tr_basis) && return nothing
     vars_gen = setdiff(vars, tr_basis)
     sol = solve_zerodim(eqs, vars_gen; dropmultiplicity=dropmultiplicity, warns=warns)
-
+    sol === nothing && return nothing
     for roots in sol
         for x in tr_basis
             roots[x] = x