Fixes for SumOfSquares

src/algebra_elts.jl

@@ -18,6 +18,12 @@ end
 (a::AlgebraElement)(x) = coeffs(a)[basis(a)[x]]
 Base.setindex!(a::AlgebraElement, v, idx) = a.coeffs[basis(a)[idx]] = v
 
+function nonzero_pairs(a::AlgebraElement)
+    return Base.Generator(nonzero_pairs(coeffs(a))) do (k, v)
+        return basis(a)[k], v
+    end
+end
+
 # AlgebraElement specific functions
 
 function supp(a::AlgebraElement)

src/arithmetic.jl

@@ -9,6 +9,10 @@ function _preallocate_output(op, args::Vararg{Any,N}) where {N}
     return similar(args[1], T)
 end
 
+function MA.promote_operation(::typeof(similar), ::Type{<:Vector}, ::Type{T}) where {T}
+    return Vector{T}
+end
+
 # module structure:
 
 Base.:*(X::AlgebraElement, a::Number) = a * X
@@ -25,12 +29,21 @@ function Base.:div(X::AlgebraElement, a::Number)
     return MA.operate_to!(_preallocate_output(div, X, a), div, X, a)
 end
 
+function MA.promote_operation(
+    op::Union{typeof(+),typeof(-)},
+    ::Type{AlgebraElement{A,T,VT}},
+    ::Type{AlgebraElement{A,S,VS}},
+) where {A,T,VT,S,VS}
+    U = MA.promote_operation(op, T, S)
+    return AlgebraElement{A,U,MA.promote_operation(similar, VT, U)}
+end
 function Base.:+(X::AlgebraElement, Y::AlgebraElement)
     return MA.operate_to!(_preallocate_output(+, X, Y), +, X, Y)
 end
 function Base.:-(X::AlgebraElement, Y::AlgebraElement)
     return MA.operate_to!(_preallocate_output(-, X, Y), -, X, Y)
 end
+
 function Base.:*(X::AlgebraElement, Y::AlgebraElement)
     return MA.operate_to!(_preallocate_output(*, X, Y), *, X, Y)
 end
@@ -92,7 +105,7 @@ function MA.operate_to!(
     X::AlgebraElement,
     Y::AlgebraElement,
 )
-    @assert parent(res) === parent(X) === parent(Y)
+    @assert parent(res) == parent(X) == parent(Y)
     MA.operate_to!(coeffs(res), -, coeffs(X), coeffs(Y))
     return res
 end
@@ -138,3 +151,8 @@ function unsafe_push!(a::Vector, k, v)
     a[k] = MA.add!!(a[k], v)
     return a
 end
+
+function unsafe_push!(a::AlgebraElement, k, v)
+    unsafe_push!(coeffs(a), basis(a)[k], v)
+    return a
+end

src/bases.jl

@@ -36,8 +36,8 @@ Implicit bases are not stored in memory and can be potentially infinite.
 """
 abstract type ImplicitBasis{T,I} <: AbstractBasis{T,I} end
 
-function zero_coeffs(::Type{S}, ::ImplicitBasis{T}) where {S,T}
-    return SparseCoefficients(T[], S[])
+function zero_coeffs(::Type{S}, ::ImplicitBasis{T,I}) where {S,T,I}
+    return SparseCoefficients(I[], S[])
 end
 
 """
@@ -77,3 +77,28 @@ function coeffs!(res, cfs, source::AbstractBasis, target::AbstractBasis)
     end
     return res
 end
+
+"""
+    adjoint_coeffs(cfs, source, target)
+Return `A' * cfs` where `A` is the linear map applied by
+`coeffs`.
+"""
+function adjoint_coeffs(cfs, source::AbstractBasis, target::AbstractBasis)
+    source === target && return cfs
+    source == target && return cfs
+    res = zero_coeffs(valtype(cfs), source)
+    return adjoint_coeffs!(res, cfs, source, target)
+end
+
+function adjoint_coeffs!(res, cfs, source::AbstractBasis, target::AbstractBasis)
+    MA.operate!(zero, res)
+    for (k, v) in nonzero_pairs(cfs)
+        x = target[k]
+        # If `x` is not in `source` then the corresponding row in `A` is zero
+        # so the column in `A'` is zero hence we can ignore it.
+        if x in source
+            res[source[x]] += v
+        end
+    end
+    return res
+end

src/sparse_coeffs.jl

@@ -42,6 +42,14 @@ function Base.zero(sc::SparseCoefficients)
     return SparseCoefficients(empty(keys(sc)), empty(values(sc)))
 end
 
+function MA.promote_operation(
+    ::typeof(similar),
+    ::Type{SparseCoefficients{K,V,Vk,Vv}},
+    ::Type{T},
+) where {K,V,Vk,Vv,T}
+    return SparseCoefficients{K,T,Vk,MA.promote_operation(similar, Vv, T)}
+end
+
 function Base.similar(s::SparseCoefficients, ::Type{T} = valtype(s)) where {T}
     return SparseCoefficients(similar(s.basis_elements), similar(s.values, T))
 end

src/types.jl

@@ -25,20 +25,25 @@ struct StarAlgebra{O,T,B<:AbstractBasis{T}} <: AbstractStarAlgebra{O,T}
 end
 
 basis(A::StarAlgebra) = A.basis
+MA.promote_operation(::typeof(basis), ::Type{StarAlgebra{O,T,B}}) where {O,T,B} = B
 object(A::StarAlgebra) = A.object
 
+function algebra end
+
 struct AlgebraElement{A,T,V} <: MA.AbstractMutable
     coeffs::V
     parent::A
 end
 
 Base.parent(a::AlgebraElement) = a.parent
-Base.eltype(a::AlgebraElement) = valtype(coeffs(a))
+Base.eltype(a::AlgebraElement) = eltype(typeof(a))
+Base.eltype(::Type{<:AlgebraElement{A,T}}) where {A,T} = T
 coeffs(a::AlgebraElement) = a.coeffs
 function coeffs(x::AlgebraElement, b::AbstractBasis)
     return coeffs(coeffs(x), basis(x), b)
 end
 basis(a::AlgebraElement) = basis(parent(a))
+MA.promote_operation(::typeof(basis), ::Type{<:AlgebraElement{A}}) where {A} = MA.promote_operation(basis, A)
 
 function AlgebraElement(coeffs, A::AbstractStarAlgebra)
     _sanity_checks(coeffs, A)