Add support for orthonormal basis

src/dependence.jl

@@ -219,6 +219,21 @@ function Base.convert(
     return BasisDependence{StaircaseDependence}(d, d.basis)
 end
 
+_full_basis(basis::MB.AbstractPolynomialBasis) = basis, eachindex(basis)
+
+function _full_basis(basis::MB.AbstractMonomialBasis)
+    full = MB.maxdegree_basis(
+        typeof(basis),
+        MP.variables(basis),
+        MP.maxdegree(basis.monomials),
+    )
+    index_map = Int[
+        something(_index(basis, full[i]), 0)
+        for i in eachindex(full)
+    ]
+    return full, index_map
+end
+
 """
     function BasisDependence{StaircaseDependence}(
         is_dependent::Function,
@@ -236,68 +251,91 @@ Foundations of Computational Mathematics 8 (2008): 607-647.
 """
 function BasisDependence{StaircaseDependence}(
     r,
-    basis::MB.MonomialBasis{M},
-) where {M}
-    if isempty(basis.monomials)
+    basis::MB.AbstractPolynomialBasis,
+)
+    if isempty(basis)
         return BasisDependence(basis, StaircaseDependence[])
     end
-    function dependence(mono)
-        i = _index(basis, mono)
-        return if isnothing(i)
+    vars = MP.variables(basis)
+    full_basis, index_map = _full_basis(basis)
+    function dependence(full_index)
+        i = index_map[full_index]
+        return if iszero(i)
             TRIVIAL
         else
             is_dependent!(r, i) ? DEPENDENT : INDEPENDENT
         end
     end
-    vars = MP.variables(basis)
-    full_basis =
-        MB.maxdegree_basis(typeof(basis), vars, MP.maxdegree(basis.monomials))
-    d = StaircaseDependence[]
     # This sieve of [LLR08, Algorithm 1] is a performance improvement but not only.
     # It also ensures that the standard monomials have the "staircase structure".
-    function is_corner_multiple(mono, indices, dependence)
-        for i in eachindex(dependence)
-            if is_dependent(dependence[i]) &&
-               MP.divides(full_basis.monomials[indices[i]], mono)
+    corners = Int[]
+    function is_corner_multiple(mono)
+        for corner in corners
+            if MP.divides(full_basis[corner], mono) # TODO tol
                 return true
             end
         end
         return false
     end
-    keep = Int[]
     # Compute standard monomials and corners
-    for (i, mono) in enumerate(full_basis.monomials)
-        if !is_corner_multiple(mono, keep, d)
-            push!(keep, i)
-            push!(d, StaircaseDependence(true, dependence(mono)))
+    keep = falses(length(full_basis))
+    d = Vector{StaircaseDependence}(undef, length(full_basis))
+    for (i, mono) in enumerate(full_basis)
+        if !is_corner_multiple(mono)
+            keep[i] = true
+            dep = dependence(i)
+            d[i] = StaircaseDependence(true, dep)
+            if is_dependent(dep)
+                push!(corners, i)
+            end
         end
     end
-    full_basis = typeof(full_basis)(full_basis.monomials[keep])
-    new_basis = MB.MonomialBasis(
-        eltype(basis.monomials)[
-            full_basis.monomials[i] * shift for
-            i in eachindex(d) if !is_dependent(d[i]) for shift in vars
-        ],
-    )
-    full_basis, I1, I2 = MB.merge_bases(full_basis, new_basis)
-    deps = Vector{StaircaseDependence}(undef, length(full_basis.monomials))
-    for (i, mono) in enumerate(full_basis.monomials)
-        if iszero(I1[i])
-            @assert !iszero(I2[i])
-            if is_corner_multiple(mono, 1:(i-1), view(deps, 1:(i-1)))
-                std = false
-            else
-                # If it was not seen before, it means it is outside the basis
-                # so it is trivial standard
-                @assert isnothing(_index(basis, mono))
-                std = true
+    new_basis = eltype(MB.generators(full_basis))[]
+    new_deps = StaircaseDependence[]
+    for (i, mono) in enumerate(full_basis)
+        if keep[i] && is_standard(d[i])
+            for shift in vars
+                shifted = mono * shift
+                j = _index(full_basis, shifted)
+                if isnothing(j)
+                    push!(new_basis, shifted)
+                    push!(new_deps, StaircaseDependence(
+                        is_corner_multiple(shifted),
+                        TRIVIAL,
+                    ))
+                else
+                    if keep[j]
+                        continue
+                    end
+                    keep[j] = true
+                    d[j] = StaircaseDependence(false, dependence(j))
+                end
             end
-            deps[i] = StaircaseDependence(std, dependence(mono))
+        end
+    end
+    I = findall(keep)
+    bd = BasisDependence(full_basis[I], d[I])
+    if isempty(new_basis)
+        return bd
+    else
+        return MB.merge_bases(
+            bd,
+            BasisDependence(MB.MonomialBasis(new_basis), new_deps),
+        )
+    end
+end
+
+function MB.merge_bases(b1::BasisDependence, b2::BasisDependence)
+    basis, I1, I2 = MB.merge_bases(b1.basis, b2.basis)
+    deps = Vector{StaircaseDependence}(undef, length(basis))
+    for i in eachindex(basis)
+        if iszero(I1[i])
+            deps[i] = b2.dependence[I2[i]]
         else
-            deps[i] = d[I1[i]]
+            deps[i] = b1.dependence[I1[i]]
         end
     end
-    return BasisDependence(full_basis, deps)
+    return BasisDependence(basis, deps)
 end
 
 function _exponents(monos, i)

test/null.jl

@@ -18,13 +18,16 @@ function test_partial_commutative_fix(x, y)
         0 0 1 0 # x * y
         1 0 0 0 # x^2
     ]
-    basis = MB.MonomialBasis(MP.monomials((x, y), 0:2))
-    null = MM.MacaulayNullspace(matrix, basis, 1e-8)
-    D = MM.StaircaseDependence
-    solver = MM.StaircaseSolver{Float64}(max_partial_iterations = 1)
-    shift_solver = MM.ShiftNullspace{D}(solver)
-    sol = MM.solve(null, shift_solver)
-    return testelements(sol, [[1, 1], [-1, 1]], 1e-8)
+    monos = MP.monomials((x, y), 0:2)
+    for basis in [MB.MonomialBasis(monos), MB.OrthonormalCoefficientsBasis(monos)]
+        null = MM.MacaulayNullspace(matrix, basis, 1e-8)
+        D = MM.StaircaseDependence
+        solver = MM.StaircaseSolver{Float64}(max_partial_iterations = 1)
+        shift_solver = MM.ShiftNullspace{D}(solver)
+        sol = MM.solve(null, shift_solver)
+        testelements(sol, [[1, 1], [-1, 1]], 1e-8)
+    end
+    return
 end
 
 function test_dependent_border(x, y)