Make exp/measure/rational implementation independent

Author: blegat

src/MultivariatePolynomials.jl

@@ -6,9 +6,23 @@ import Base: *, +, -, /, ^, ==,
     promote_rule, convert, show, isless, size, getindex,
     one, zero, transpose, isapprox, @pure, dot, copy
 
-include("utils.jl")
-include("types.jl")
-include("operators.jl")
+abstract type AbstractPolynomialLike{T} end
+abstract type AbstractTermLike{T} <: AbstractPolynomialLike{T} end
+abstract type AbstractMonomialLike <: AbstractTermLike{Int} end
+
+abstract type AbstractVariable <: AbstractMonomialLike end
+abstract type AbstractMonomial <: AbstractMonomialLike end
+abstract type AbstractTerm{T} <: AbstractTermLike{T} end
+abstract type AbstractPolynomial{T} <: AbstractPolynomialLike{T} end
+
+const APL{T} = AbstractPolynomialLike{T}
+
+include("measure.jl")
+include("exp.jl")
+
+#include("utils.jl")
+#include("types.jl")
+#include("operators.jl")
 #include("show.jl")
 #include("substitution.jl")
 

src/exp.jl

@@ -1,9 +1,9 @@
 export expectation
 
-function _dot{C}(m::Measure{C}, p::TermContainer{C}, f)
+function _dot(m::Measure, p::APL, f)
     i = 1
     s = 0
-    for t in p
+    for t in terms(p)
         while i <= length(m.x) && t.x != m.x[i]
             i += 1
         end
@@ -15,11 +15,8 @@ function _dot{C}(m::Measure{C}, p::TermContainer{C}, f)
     end
     s
 end
-dot(m::Measure, p::TermContainer) = _dot(m, p, (*))
-dot(p::TermContainer, m::Measure) = _dot(m, p, (a, b) -> b * a)
+Base.dot(m::Measure, p::APL) = _dot(m, p, (*))
+Base.dot(p::APL, m::Measure) = _dot(m, p, (a, b) -> b * a)
 
-dot(m::Measure, p::PolyType) = dot(m, TermContainer(p))
-dot(p::PolyType, m::Measure) = dot(TermContainer(p), m)
-
-expectation(m::Measure, p::PolyType) = dot(m, p)
-expectation(p::PolyType, m::Measure) = dot(p, m)
+expectation(m::Measure, p::APL) = dot(m, p)
+expectation(p::APL, m::Measure) = dot(p, m)

src/measure.jl

@@ -1,39 +1,30 @@
 export Measure, zeta, ζ
 
-type Moment{C, T}
+type Moment{T, MT <: AbstractMonomial}
     α::T
-    x::Monomial{C}
+    x::MT
 end
 
 # If a monomial is not in x, it does not mean that the moment is zero, it means that it is unknown/undefined
-type Measure{C, T}
+type Measure{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}}
     a::Vector{T}
-    x::MonomialVector{C}
+    x::MVT
 
-    function Measure{C, T}(a::Vector{T}, x::MonomialVector{C}) where {C, T}
+    function Measure{T, MT, MVT}(a::Vector{T}, x::MVT) where {T, MT, MVT}
         if length(a) != length(x)
             error("There should be as many coefficient than monomials")
         end
         new(a, x)
     end
 end
 
-Measure{C, T}(a::Vector{T}, x::MonomialVector{C}) = Measure{C, T}(a, x)
-function (::Type{Measure{C}}){C}(a::Vector, x::Vector)
-    if length(a) != length(x)
-        error("There should be as many coefficient than monomials")
-    end
-    σ, X = sortmonovec(PolyVar{C}, x)
-    Measure(a[σ], X)
-end
-Measure{T<:VectorOfPolyType{true}}(a::Vector, X::Vector{T}) = Measure{true}(a, X)
-Measure{T<:VectorOfPolyType{false}}(a::Vector, X::Vector{T}) = Measure{false}(a, X)
+Measure{T, MT <: AbstractMonomial}(a::Vector{T}, x::AbstractVector{MT}) = Measure{T, MT, monovectype(x)}(monovec(a, x)...)
 
-function ζ{C, T}(v::Vector{T}, x::MonomialVector{C}, varorder::Vector{PolyVar{C}})
+function ζ{T, MT <: AbstractMonomial, PVT <: AbstractVariable}(v::Vector{T}, x::AbstractVector{MT}, varorder::AbstractVector{PVT})
     Measure(T[m(v, varorder) for m in x], x)
 end
 
-type MatMeasure{C, T}
+type MatMeasure{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}}
     Q::Vector{T}
-    x::MonomialVector{C}
+    x::MVT
 end

src/rational.jl

@@ -1,59 +1,49 @@
 export RationalPoly
 import Base.+, Base.-, Base.*, Base./
 
-immutable RationalPoly{C, S, T} <: PolyType{C}
-    num::TermContainer{C, S}
-    den::TermContainer{C, T}
+immutable RationalPoly{NT <: APL, DT <: APL}
+    num::NT
+    den::DT
 end
-iscomm{C, S, T}(r::Type{RationalPoly{C, S, T}}) = C
 
-Base.convert{C, S, T}(::Type{RationalPoly{C, S, T}}, q::RationalPoly{C, S, T}) = q
-Base.convert{C, S, T, U, V}(::Type{RationalPoly{C, S, T}}, q::RationalPoly{C, U, V}) = TermContainer{C, S}(q.num) / TermContainer{C, T}(q.den)
-function Base.convert{C, S, T}(::Type{RationalPoly{C, S, T}}, p::TermContainer{C, S})
-    p / one(TermContainer{C, T})
+Base.convert{NT, DT}(::Type{RationalPoly{NT, DT}}, q::RationalPoly{NT, DT}) = q
+Base.convert{NTout, DTout, NTin, DTin}(::Type{RationalPoly{NTout, DTout}}, q::RationalPoly{NTin, DTin}) = convert(NTout, q.num) / convert(DTout, q.den)
+function Base.convert{NT, DT}(::Type{RationalPoly{NT, DT}}, p::NT)
+    p / one(DT)
 end
-function Base.convert{C, S, T}(::Type{RationalPoly{C, S, T}}, p::TermContainer)
-    convert(RationalPoly{C, S, T}, TermContainer{C, S}(p))
-end
-function Base.convert{C, S, T}(::Type{RationalPoly{C, S, T}}, p)
-    Base.convert(RationalPoly{C, S, T}, TermContainer{C, S}(p))
+function Base.convert{NT, DT}(::Type{RationalPoly{NT, DT}}, p)
+    convert(RationalPoly{NT, DT}, convert(NT, p))
 end
 
-(/)(r::RationalPoly, p::TermContainer) = r.num / (r.den * p)
-function (/){C, S, T}(num::TermContainer{C, S}, den::TermContainer{C, T})
-    RationalPoly{C, S, T}(num, den)
+(/)(r::RationalPoly, p) = r.num / (r.den * p)
+function (/){NT <: APL, DT <: APL}(num::NT, den::DT)
+    RationalPoly{NT, DT}(num, den)
 end
-function (/){C}(num, den::PolyType{C})
-    TermContainer{C}(num) / den
+function (/)(num, den::APL)
+    term(num) / den
 end
-(/){C}(num::PolyType{C}, den::PolyType{C}) = TermContainer{C}(num) / TermContainer{C}(den)
-
 # Polynomial divided by coefficient is a polynomial not a rational polynomial
-(/){C}(num::PolyType{C}, den) = num * (1 / den)
+# (1/den) * num would not be correct in case of noncommutative coefficients
+(/)(num::APL, den) = num * (1 / den)
 
 function (+)(r::RationalPoly, s::RationalPoly)
     (r.num*s.den + r.den*s.num) / (r.den * s.den)
 end
-function (+)(p::TermContainer, r::RationalPoly)
+function (+)(p, r::RationalPoly)
     (p*r.den + r.num) / r.den
 end
-(+)(r::RationalPoly, p::Polynomial) = p + r
-(+)(r::RationalPoly, t::Term) = t + r
+(+)(r::RationalPoly, p) = p + r
 function (-)(r::RationalPoly, s::RationalPoly)
     (r.num*s.den - r.den*s.num) / (r.den * s.den)
 end
-(-)(p::PolyType, s::RationalPoly) = (p * s.den - s.num) / s.den
-(-)(s::RationalPoly, p::PolyType) = (s.num - p * s.den) / s.den
+(-)(p, s::RationalPoly) = (p * s.den - s.num) / s.den
+(-)(s::RationalPoly, p) = (s.num - p * s.den) / s.den
 
 (*)(r::RationalPoly, s::RationalPoly) = (r.num*s.num) / (r.den*s.den)
-(*)(p::TermContainer, r::RationalPoly) = p == r.den ? r.num : (p * r.num) / r.den
-(*)(r::RationalPoly, p::Polynomial) = p == r.den ? r.num : (r.num * p) / r.den
-(*)(r::RationalPoly, t::Term)          = t == r.den ? r.num : (r.num * t) / r.den
-(*)(p::PolyType, r::RationalPoly) = TermContainer(p) * r
-(*)(r::RationalPoly, p::Monomial) = r * TermContainer(p)
-(*)(r::RationalPoly, p::PolyVar)  = r * TermContainer(p)
-(*){C}(α, r::RationalPoly{C}) = TermContainer{C}(α) * r
-(*){C}(r::RationalPoly{C}, α) = r * TermContainer{C}(α)
+# Not type stable, currently it is a hack for SumOfSquares/test/SOSdemo2.jl:line 22
+# We should take gcd between numerator and denominator instead and in sosdemo2, we should cast to polynomial manually
+(*)(p, r::RationalPoly)       = p == r.den ? r.num : (p * r.num) / r.den
+(*)(r::RationalPoly, p)       = p == r.den ? r.num : (r.num * p) / r.den
 
-zero(r::RationalPoly) = zero(r.num)
-zero{C, S, T}(::Type{RationalPoly{C, S, T}}) = zero(Polynomial{C, S})
+zero{NT}(::RationalPoly{NT}) = zero(NT)
+zero{NT, DT}(::Type{RationalPoly{NT, DT}}) = zero(NT)

test/alltests.jl

@@ -0,0 +1,20 @@
+#include("mono.jl")
+#include("ncmono.jl")
+#include("poly.jl")
+#include("rational.jl")
+include("measure.jl")
+include("exp.jl")
+#include("promote.jl")
+#include("comp.jl")
+#include("nccomp.jl")
+#include("alg.jl")
+#include("ncalg.jl")
+#include("diff.jl")
+#include("subs.jl")
+#include("algebraicset.jl")
+#include("hash.jl")
+
+include("show.jl")
+
+#include("example1.jl")
+#include("example2.jl")

test/measure.jl

@@ -1,7 +1,7 @@
 @testset "Measure" begin
     @polyvar x y
     @test_throws ErrorException Measure([1, 2], [x, x*y, y])
-    @test_throws ErrorException Measure([1, 2, 3, 4], MonomialVector([x, x*y, y]))
+    @test_throws ErrorException Measure([1, 2, 3, 4], [x, x*y, y])
     m = Measure([1, 0, 2, 3], [x^2*y^2, y*x^2, x*y*x^2, x*y^2])
     @test m.a == [2, 1, 0, 3]
 end

test/mono.jl

@@ -1,11 +1,11 @@
 @testset "PolyVar and Monomial tests" begin
-#   @testset "polyvar macro index set" begin
-#       n = 3
-#       @polyvar x[1:n] y z[1:n-1]
-#       @test length(x) == 3
-#       @test length(z) == 2
-#       @test x[1] > x[2] > x[3] > y > z[1] > z[2]
-#   end
+    @testset "polyvar macro index set" begin
+        n = 3
+        @polyvar x[1:n] y z[1:n-1]
+        @test length(x) == 3
+        @test length(z) == 2
+        @test x[1] > x[2] > x[3] > y > z[1] > z[2]
+    end
     @testset "PolyVar" begin
         @polyvar x
         @test copy(x) == x
@@ -21,19 +21,19 @@
         @inferred zero(x^2)
         @test one(x^2) == 1
         @inferred one(x^2)
-#       @polyvar y[1:7]
-#       m = y[1] * y[3] * y[5] * y[7]
-#       @test issorted(vars(y[2] * m), rev=true)
-#       @test issorted(vars(m * y[4]), rev=true)
-#       @test issorted(vars(y[6] * m), rev=true)
+        @polyvar y[1:7]
+        m = y[1] * y[3] * y[5] * y[7]
+        @test issorted(vars(y[2] * m), rev=true)
+        @test issorted(vars(m * y[4]), rev=true)
+        @test issorted(vars(y[6] * m), rev=true)
+    end
+    @testset "MonomialVector" begin
+        @polyvar x y
+        X = [x^2,x*y,y^2]
+        for (i, m) in enumerate(monomials([x, y], 2))
+            @test m == X[i]
+        end
     end
-#   @testset "MonomialVector" begin
-#       @polyvar x y
-#       X = [x^2,x*y,y^2]
-#       for (i, m) in enumerate(monomials([x, y], 2))
-#           @test m == X[i]
-#       end
-#   end
 end
 module newmodule
     using Base.Test