fix: use term() guards to match master behavior for commutative ops

Adds SymbolicUtils.term() usage to prevent SUv4 from canonicalizing
x*x to x^2 during conversion. This makes the backport behaviorally
equivalent to upstream/master.

- Import term from SymbolicUtils
- Use term(op, ...) for * and + in parse_tree_to_eqs
- Use term(*, ...) in multiply_powers to build explicit products
- Restore test to check x1*x1 stays as explicit multiplication

Author: MilesCranmerBot

ext/DynamicExpressionsSymbolicUtilsExt.jl

@@ -8,7 +8,7 @@ using DynamicExpressions.OperatorEnumModule: AbstractOperatorEnum
 using DynamicExpressions.UtilsModule: deprecate_varmap
 
 using SymbolicUtils
-using SymbolicUtils: BasicSymbolic, SymReal, iscall, issym, isconst, unwrap_const
+using SymbolicUtils: BasicSymbolic, SymReal, iscall, issym, isconst, unwrap_const, term
 
 import DynamicExpressions.ExtensionInterfaceModule: node_to_symbolic, symbolic_to_node
 import DynamicExpressions.ValueInterfaceModule: is_valid
@@ -69,6 +69,15 @@ function parse_tree_to_eqs(
     # Convert children to symbolic form
     sym_children = map(x -> parse_tree_to_eqs(x, operators, index_functions), children)
 
+    # SymbolicUtils v4 may canonicalize some commutative operations at construction time
+    # (e.g. `x*x` -> `x^2`, or reordering `a + b`).
+    #
+    # For stable round-trips (and to avoid introducing `^` when it's not in the operator set),
+    # construct commutative ops as explicit terms.
+    if op === (*) || op === (+)
+        return subs_bad(term(op, sym_children...))
+    end
+
     return subs_bad(op(sym_children...))
 end
 
@@ -320,11 +329,22 @@ function multiply_powers(
             if n_val == 1
                 return l, true
             elseif n_val == -1
-                return 1.0 / l, true
+                return term(/, 1.0, l), true
             elseif n_val > 1
-                return reduce(*, [l for i in 1:n_val]), true
+                # IMPORTANT: use `term(*, ...)` to prevent SymbolicUtils from immediately
+                # canonicalizing `l*l` back into `l^2`.
+                out = l
+                for _ in 2:n_val
+                    out = term(*, out, l)
+                end
+                return out, true
             elseif n_val < -1
-                return reduce(/, vcat([1], [l for i in 1:abs(n_val)])), true
+                # Build 1/(l*l*...) using explicit multiplication terms.
+                denom = l
+                for _ in 2:abs(n_val)
+                    denom = term(*, denom, l)
+                end
+                return term(/, 1.0, denom), true
             else
                 return 1.0, true
             end
@@ -340,6 +360,11 @@ function multiply_powers(
         r, complete2 = multiply_powers(args[2])
         @return_on_false complete2 eqn
         @return_on_false is_valid(r) eqn
+        # SymbolicUtils v4 normalizes `x*x` into `x^2` via the `*` method; preserve
+        # explicit multiplication terms so we don't introduce `^` during conversion.
+        if op == *
+            return term(op, l, r), true
+        end
         return op(l, r), true
     else
         # return tree_mapreduce(multiply_powers, op, args)
@@ -351,7 +376,8 @@ function multiply_powers(
         end
         cumulator = out[1][1]
         for i in 2:size(out, 1)
-            cumulator = op(cumulator, out[i][1])
+            cumulator =
+                (op == *) ? term(op, cumulator, out[i][1]) : op(cumulator, out[i][1])
             @return_on_false is_valid(cumulator) eqn
         end
         return cumulator, true

test/test_simplification.jl

@@ -41,36 +41,37 @@ tree = convert(Node, eqn2, operators)
 # Make sure the other node is x1:
 @test (!tree.l.constant ? tree.l : tree.r).feature == 1
 
-# SymbolicUtils v4 automatically simplifies x1*x1 to x1^2
-# For round-trip to work, we need ^ in the operator set
-operators_with_pow = OperatorEnum(; binary_operators=(+, -, /, *, ^))
+# Finally, let's try converting a product, and ensure
+# that SymbolicUtils does not convert it to a power:
 tree = Node("x1") * Node("x1")
-eqn = convert(BasicSymbolic, tree, operators_with_pow)
-# The symbolic repr will be x1^2 in SymbolicUtils v4
-@test occursin("x1", repr(eqn))
-# Test converting back (x^2 comes back as x^2 since ^ is in operators):
-tree_copy = convert(Node, eqn, operators_with_pow)
-# The structure is preserved as a power in v4
-@test occursin("x1", repr(tree_copy))
-
-# Let's test a more complex function with supported operators
-# (Custom operators are not supported in SymbolicUtils v4+)
-operators = OperatorEnum(; binary_operators=(+, *, -, /), unary_operators=(cos, exp, sin))
+eqn = convert(BasicSymbolic, tree, operators)
+@test repr(eqn) ≈ "x1*x1"
+# Test converting back:
+tree_copy = convert(Node, eqn, operators)
+@test repr(tree_copy) ≈ "(x1*x1)"
+
+# Let's test a more complex function. In SymbolicUtils v4+, custom operators need
+# `index_functions=true` to round-trip.
 
 x1, x2, x3 = Node("x1"), Node("x2"), Node("x3")
+pow_abs2(x, y) = abs(x)^y
+
+operators = OperatorEnum(;
+    binary_operators=(+, *, -, /, pow_abs2), unary_operators=(custom_cos, exp, sin)
+)
 @extend_operators operators
 tree = (
-    ((x2 + x2) * ((-0.5982493 / (x1 * x2)) / -0.54734415)) + (
+    ((x2 + x2) * ((-0.5982493 / pow_abs2(x1, x2)) / -0.54734415)) + (
         sin(
-            cos(
+            custom_cos(
                 sin(1.2926733 - 1.6606787) /
                 sin(((0.14577048 * x1) + ((0.111149654 + x1) - -0.8298334)) - -1.2071426),
-            ) * (cos(x3 - 2.3201916) + ((x1 - (x1 * x2)) / x2)),
-        ) / (0.14854191 - ((cos(x2) * -1.6047639) - 0.023943262))
+            ) * (custom_cos(x3 - 2.3201916) + ((x1 - (x1 * x2)) / x2)),
+        ) / (0.14854191 - ((custom_cos(x2) * -1.6047639) - 0.023943262))
     )
 )
 # Convert to symbolic form
-eqn = convert(BasicSymbolic, tree, operators)
+eqn = convert(BasicSymbolic, tree, operators; index_functions=true)
 
 tree_copy = convert(Node, eqn, operators)
 tree_copy2 = convert(Node, simplify(eqn), operators)