SciML
diff --git a/‎.github/workflows/Downgrade.yml‎
Lines changed: 4 additions & 0 deletions b/‎.github/workflows/Downgrade.yml‎
Lines changed: 4 additions & 0 deletions
diff --git a/‎Project.toml‎
Lines changed: 3 additions & 9 deletions b/‎Project.toml‎
Lines changed: 3 additions & 9 deletions
diff --git a/‎docs/Project.toml‎
Lines changed: 2 additions & 2 deletions b/‎docs/Project.toml‎
Lines changed: 2 additions & 2 deletions
diff --git a/‎src/SymbolicAnalysis.jl‎
Lines changed: 12 additions & 2 deletions b/‎src/SymbolicAnalysis.jl‎
Lines changed: 12 additions & 2 deletions
diff --git a/‎src/atoms.jl‎
Lines changed: 61 additions & 6 deletions b/‎src/atoms.jl‎
Lines changed: 61 additions & 6 deletions
diff --git a/‎src/canon.jl‎
Lines changed: 4 additions & 8 deletions b/‎src/canon.jl‎
Lines changed: 4 additions & 8 deletions
diff --git a/‎src/gdcp/gdcp_rules.jl‎
Lines changed: 12 additions & 4 deletions b/‎src/gdcp/gdcp_rules.jl‎
Lines changed: 12 additions & 4 deletions
diff --git a/‎src/gdcp/lorentz.jl‎
Lines changed: 36 additions & 18 deletions b/‎src/gdcp/lorentz.jl‎
Lines changed: 36 additions & 18 deletions
@@ -15,4 +15,8 @@ jobs:
     uses: "SciML/.github/.github/workflows/downgrade.yml@v1"
     with:
       julia-version: "lts"
+      # Downgrade the test project too: its deps (Optimization/Zygote/ForwardDiff)
+      # have floors that pull transitive deps (SciMLBase) far above the package's
+      # own minima, so pinning only the root project leaves an unsatisfiable env.
+      projects: ".,test"
     secrets: "inherit"
@@ -6,34 +6,28 @@ version = "0.3.7"
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
-IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 DSP = "0.7, 0.8"
 DataStructures = "0.18, 0.19"
-Dictionaries = "0.4"
 Distributions = "0.25"
 DomainSets = "0.7, 0.8"
-IfElse = "0.1"
 LinearAlgebra = "1.10"
 LogExpFunctions = "0.3, 1.0"
-Manifolds = "0.9, 0.10, 0.11"
+Manifolds = "0.10, 0.11"
 PDMats = "0.11"
 PrecompileTools = "1"
-RecursiveArrayTools = "3, 4"
 StatsBase = "0.34"
-SymbolicUtils = "3.1.2, 4.33"
-Symbolics = "6, 7.25"
+SymbolicUtils = "4.35"
+Symbolics = "7.28"
 julia = "1.10"
@@ -33,5 +33,5 @@ PDMats = "0.11"
 RecursiveArrayTools = "3, 4.3"
 StatsBase = "0.34"
 SymbolicAnalysis = "0.3.6"
-SymbolicUtils = "3, 4.33"
-Symbolics = "6, 7.25"
+SymbolicUtils = "4.35"
+Symbolics = "7.28"
@@ -9,11 +9,21 @@ using Distributions
 using DSP, DataStructures
 
 using Symbolics
-import Symbolics: Symbolic, issym, Term
-using SymbolicUtils: iscall
+import Symbolics: issym, Term
+using SymbolicUtils: iscall, BasicSymbolic
 using SymbolicUtils.Rewriters
 SymbolicUtils.inspect_metadata[] = true
 
+# Symbolics v7 / SymbolicUtils v4 removed the `Symbolic` abstract type: every
+# symbolic — scalar or array — is now a `BasicSymbolic{SymReal}`.
+const Symbolic = BasicSymbolic
+
+# The scalar-symbolic type Symbolics uses when dispatching `in(::symbolic, ::Domain)`.
+# Matching it exactly lets the `in(::_, ::CustomDomain)` disambiguators below stay
+# strictly more specific than Symbolics' `IntervalSets.Domain` method (which is
+# keyed on `BasicSymbolic{SymReal}`).
+const InDomainSymbolic = BasicSymbolic{SymbolicUtils.SymReal}
+
 struct VarDomain end
 
 include("rules.jl")
 
@@ -231,7 +231,10 @@ function quad_over_lin(x::AbstractVector{<:Real}, y::Real)
     return sum(x .^ 2) / y
 end
 
-Symbolics.@register_symbolic quad_over_lin(x::AbstractVector, y::Real) false
+# On Symbolics v7 both scalar- and array-valued registered symbolics share the
+# concrete type `BasicSymbolic{SymReal}`, so registering the vector form would
+# generate the same `quad_over_lin(::BasicSymbolic{SymReal}, ::Real)` method as
+# the scalar one below and collide. The scalar registration alone covers both.
 
 """
     quad_over_lin(x::Real, y::Real)
@@ -393,13 +396,57 @@ add_dcprule(imag, ℂ, AnySign, Affine, AnyMono)
 add_dcprule(inv, HalfLine{Real, :open}(), Positive, Convex, Decreasing)
 add_dcprule(log, HalfLine{Real, :open}(), AnySign, Concave, Increasing)
 
-@register_symbolic Base.log(A::Symbolics.Arr)
+# Matrix-valued atoms (`log`, `inv`, `sqrt` of a symbolic matrix). On Symbolics
+# v7 / SymbolicUtils v4 a symbolic matrix unwraps to `BasicSymbolic{SymReal}` —
+# the same concrete type as a symbolic scalar — so a `@register_symbolic` macro
+# would emit a `f(::BasicSymbolic{SymReal})` method that overwrites SymbolicUtils'
+# own scalar `f` and aborts precompilation. Build the matrix term directly off the
+# `Arr` wrapper instead via `SymbolicUtils.term`, which leaves the scalar methods
+# untouched. Some scalar atoms (e.g. `sqrt`) register a `promote_shape` that
+# rejects matrix shapes, which `term` would invoke unless the shape is supplied up
+# front, so pass `shape` explicitly.
+function matrix_atom(f, A::Symbolics.Arr)
+    a = Symbolics.unwrap(A)
+    return Symbolics.wrap(
+        SymbolicUtils.term(
+            f, a; type = SymbolicUtils.symtype(a), shape = Symbolics.shape(a)
+        )
+    )
+end
+function matrix_atom(f, A::AbstractMatrix{<:Num})
+    a = Symbolics.unwrap.(A)
+    return Symbolics.wrap(
+        SymbolicUtils.term(f, a; type = Matrix{Real}, shape = map(Base.OneTo, size(A)))
+    )
+end
+
+# A symbolic matrix–matrix product returns a bare `BasicSymbolic`; re-wrap the
+# 2-argument `Arr` product so it reaches the matrix-atom methods below (which
+# dispatch on `Arr`), restoring the array shape. This 2-arg method is more
+# specific than Symbolics' variadic `*`, so it does not overwrite it.
+function Base.:*(x::Symbolics.Arr{<:Any, 2}, y::Symbolics.Arr{<:Any, 2})
+    return Symbolics.wrap(Symbolics.unwrap(x) * Symbolics.unwrap(y))
+end
+
+# SymbolicUtils v4 gives `log` a matrix-permissive `promote_shape` but leaves
+# `sqrt` scalar-only, so rewriting a matrix `sqrt` term (e.g. during the
+# analysis walk's `maketerm`) throws "Invalid shapes for sqrt". Add a matrix
+# rule for `sqrt`; dispatching on the concrete `ShapeVecT` (rather than the
+# `ShapeT` union SymbolicUtils uses) makes this strictly more specific, so it
+# extends rather than overwrites the existing method.
+function SymbolicUtils.promote_shape(::typeof(sqrt), sh::SymbolicUtils.ShapeVecT)
+    (length(sh) == 0 || length(sh) == 2) && return sh
+    return SymbolicUtils._throw_array(sqrt, sh)
+end
+
+Base.log(A::Symbolics.Arr) = matrix_atom(log, A)
+Base.log(A::Matrix{Num}) = matrix_atom(log, A)
 add_dcprule(log, array_domain(RealLine(), 2), Positive, Concave, Increasing)
 
-@register_symbolic LinearAlgebra.inv(A::Symbolics.Arr)
+LinearAlgebra.inv(A::Symbolics.Arr) = matrix_atom(inv, A)
 add_dcprule(inv, semidefinite_domain(), AnySign, Convex, Decreasing)
 
-@register_symbolic LinearAlgebra.sqrt(A::Symbolics.Arr)
+LinearAlgebra.sqrt(A::Symbolics.Arr) = matrix_atom(sqrt, A)
 add_dcprule(sqrt, semidefinite_domain(), Positive, Concave, Increasing)
 
 add_dcprule(
@@ -420,7 +467,7 @@ Returns the log of the normal cumulative distribution function of `x`.
     - `x::Real`: A Real.
 """
 function lognormcdf(x::Real)
-    return logcdf(Normal, x)
+    return logcdf(Normal(), x)
 end
 Symbolics.@register_symbolic lognormcdf(x::Real)
 add_dcprule(lognormcdf, RealLine(), Negative, Concave, Increasing)
@@ -434,6 +481,11 @@ add_dcprule(min, (RealLine(), RealLine()), AnySign, Concave, Increasing)
 
 # special cases which depend on arguments:
 function dcprule(::typeof(^), x::Symbolic, i)
+    # A literal exponent is wrapped as a constant `BasicSymbolic`, so
+    # `isinteger`/`isone`/comparisons below would operate on a symbolic and error;
+    # `Symbolics.value` unwraps it to the underlying number (identity for an
+    # already-numeric exponent).
+    i = Symbolics.value(i)
     args = (x, i)
     if isone(i)
         return makerule(RealLine(), AnySign, Affine, Increasing), args
@@ -512,7 +564,10 @@ add_dcprule(vec, array_domain(RealLine(), 2), AnySign, Affine, Increasing)
 add_dcprule(vcat, array_domain(array_domain(RealLine(), 1), 1), AnySign, Affine, Increasing)
 
 function dcprule(::typeof(broadcast), f, x...)
-    return dcprule(f, x...)
+    # The broadcasted function is wrapped as a constant symbolic (e.g.
+    # `broadcast(exp, z)` carries a symbolic `exp`); `Symbolics.value` recovers the
+    # underlying function (identity for a plain function).
+    return dcprule(Symbolics.value(f), x...)
 end
 hasdcprule(::typeof(broadcast)) = true
 
 
@@ -1,13 +1,9 @@
 function canonize(ex)
+    # Symbolics v7 flattens `x' * (Y * x)` / `(B' * X) * B` to a single
+    # scalar/matrix `*` term, so match the flattened `*` directly.
     rs = [
-        @rule (adjoint(~x) * (~Y * ~x))[1] => quad_form(~x, ~Y)
-        @rule ((adjoint(~B) * ~X) * ~B)[
-            Base.OneTo(size(~B, 2)), Base.OneTo(
-                size(
-                    ~B, 1
-                )
-            ),
-        ] => conjugation(~X, ~B)
+        @rule adjoint(~x) * ~Y * ~x => quad_form(~x, ~Y)
+        @rule adjoint(~B) * ~X * ~B => conjugation(~X, ~B)
     ]
     try
         rc = SymbolicUtils.Chain(rs)
 
@@ -48,7 +48,7 @@ function mul_gcurvature(args)
                 return GUnknownCurvature
             end
         else
-            constant_prod *= arg
+            constant_prod *= constval(arg)
         end
     end
     if non_constant_expr !== nothing
@@ -101,6 +101,8 @@ function find_gcurvature(ex)
     if iscall(ex)
         f, args = operation(ex), arguments(ex)
         knowngcurv = false
+        f_curvature = GUnknownCurvature
+        f_monotonicity = (GAnyMono,)
 
         if hasgdcprule(f) && !any(iscall.(args))
             rule, args = gdcprule(f, args...)
@@ -161,16 +163,17 @@ function find_gcurvature(ex)
                 end
             end
         elseif Symbol(f) == :*
-            if args[1] isa Number && args[1] > 0
+            a1 = constval(args[1])
+            if a1 isa Number && a1 > 0
                 return find_gcurvature(args[2])
-            elseif args[1] isa Number && args[1] < 0
+            elseif a1 isa Number && a1 < 0
                 argscurv = find_gcurvature(args[2])
                 if argscurv == GConvex
                     return GConcave
                 elseif argscurv == GConcave
                     return GConvex
                 else
-                    argscurv
+                    return argscurv
                 end
             else
                 @warn "Disciplined Programming does not support multiple non-constant arguments in multiplication"
@@ -242,6 +245,11 @@ function find_gcurvature(ex)
 end
 
 function propagate_gcurvature(ex, M::AbstractManifold)
+    # Operate on the raw symbolic: on Symbolics v7 walking a `Num`/`Arr` wrapper
+    # round-trips through wrap/unwrap and loses the gcurvature metadata that the
+    # final `getgcurvature` reads. `analyze` already unwraps; do the same here so
+    # the function is correct when called directly on a wrapped expression.
+    ex = Symbolics.unwrap(ex)
     r = [
         @rule *(~~x) => setgcurvature(~MATCH, mul_gcurvature(~~x))
         @rule +(~~x) => setgcurvature(~MATCH, add_gcurvature(~~x))
 
@@ -6,32 +6,39 @@
 
 using Manifolds
 using LinearAlgebra
-using Symbolics: Symbolic, @register_symbolic, unwrap, variables
-
-@register_symbolic Manifolds.distance(
-    M::Manifolds.Lorentz,
-    p::AbstractVector,
-    q::Union{Symbolics.Arr, AbstractVector}
-) false
+using Symbolics: @register_symbolic, unwrap, variables
+
+# See the SPD `distance` note in gdcp/spd.jl: build the term directly off the
+# symbolic point so the SPD and Lorentz methods do not collide.
+function Manifolds.distance(
+        M::Manifolds.Lorentz,
+        p::AbstractVector,
+        q::Symbolics.Arr
+    )
+    return Symbolics.wrap(
+        SymbolicUtils.term(Manifolds.distance, M, p, Symbolics.unwrap(q); type = Real)
+    )
+end
 add_gdcprule(Manifolds.distance, Manifolds.Lorentz, Positive, GConvex, GAnyMono)
 
 """
-    lorentz_log_barrier(a, p)
+    lorentz_log_barrier(p)
 
-Computes the log-barrier function for the Lorentz model: `-log(-1 - <a, p>_L)`.
+Computes the log-barrier function for the Lorentz model: `-log(-1 - <a, p>_L)`,
+with the fixed vector `a = (0, ..., 0, 1)` in R^(d+1).
 
 # Arguments
 
-    - `a`: The vector (0, ..., 0, 1) in R^(d+1).
     - `p`: A point on the Lorentz manifold.
 """
 function lorentz_log_barrier(p::AbstractVector)
-    # Lorentzian inner product: a⋅p_L = a1*p1 + ... + a_d*p_d - a_{d+1}*p_{d+1}
-    inner_prod = a[end] * p[end]
-    return -log(-1 + inner_prod)
+    # a = (0, ..., 0, 1), so the Lorentzian inner product
+    # <a, p>_L = a1*p1 + ... + a_d*p_d - a_{d+1}*p_{d+1} reduces to -p[end].
+    # The barrier is -log(-1 - <a, p>_L) = -log(-1 + p[end]).
+    return -log(-1 + p[end])
 end
 
-@register_symbolic lorentz_log_barrier(p::Union{Symbolics.Arr, AbstractVector})
+@register_symbolic lorentz_log_barrier(p::Vector{Num})
 add_gdcprule(lorentz_log_barrier, Manifolds.Lorentz, Positive, GConvex, GIncreasing)
 
 """
@@ -69,7 +76,7 @@ end
 
 @register_symbolic lorentz_homogeneous_quadratic(
     A::AbstractMatrix,
-    p::Union{Symbolics.Arr, AbstractVector}
+    p::Vector{Num}
 )
 add_gdcprule(lorentz_homogeneous_quadratic, Manifolds.Lorentz, Positive, GConvex, GAnyMono)
 
@@ -102,7 +109,7 @@ end
 
 @register_symbolic lorentz_homogeneous_diagonal(
     a::AbstractVector,
-    p::Union{Symbolics.Arr, AbstractVector}
+    p::Vector{Num}
 )
 add_gdcprule(lorentz_homogeneous_diagonal, Manifolds.Lorentz, Positive, GConvex, GAnyMono)
 
@@ -170,7 +177,18 @@ function lorentz_least_squares(X::AbstractMatrix, y::AbstractVector, p::Abstract
     return lorentz_nonhomogeneous_quadratic(A, b, c, p)
 end
 
-@register_symbolic lorentz_least_squares(X::Matrix{Num}, y::Vector{Num}, p::Vector{Num})
+# A `@register_symbolic` three-array registration explodes into a combinatorial,
+# mutually-ambiguous set of wrapper methods on Symbolics v7 (Aqua flags them; see
+# the `sdivergence` note in gdcp/spd.jl). Build the `lorentz_least_squares(X, y, p)`
+# term directly off the symbolic point `p` instead (the gDCP pass only needs
+# `operation(ex) == lorentz_least_squares`).
+function lorentz_least_squares(X::AbstractMatrix, y::AbstractVector, p::Symbolics.Arr)
+    return Symbolics.wrap(
+        SymbolicUtils.term(
+            lorentz_least_squares, X, y, Symbolics.unwrap(p); type = Real
+        )
+    )
+end
 add_gdcprule(lorentz_least_squares, Manifolds.Lorentz, Positive, GConvex, AnyMono)
 
 """
@@ -203,7 +221,7 @@ end
 
 @register_symbolic lorentz_transform(
     O::AbstractMatrix,
-    p::Union{Symbolics.Arr, AbstractVector}
+    p::Vector{Num}
 )
 # Not adding a rule since this preserves geodesic convexity but doesn't have a specific curvature