Bernstein polynomial type for quadratic clipping

Question

I would like to implement a Bernstein polynomial class in Julia for learning purposes. The goal in the end is to implement the quadratic clipping algorithm from this paper. I am completely new to Julia coming from Python + Numpy + Numba, so I want to make sure that I program in the canonical Julia way when using Julia.

The naming of variables follows the aforementioned paper.

using Polynomials

type Bernstein
    n::Integer
    i::Integer
    α::Real
    β::Real
end
Bernstein(n, i) = Bernstein(n, i, 0., 1.)

Base.string(B::Bernstein) = "Bernstein $(get_p(B))"
Base.show(io::IO, B::Bernstein) = print(io, string(B))
Base.display(io::IO, B::Bernstein) = print(io, string(B))

function get_p(n, i, α=0, β=1)
    @assert β >= α "Invalid interval"
    return (binomial(n, i)
            * poly([α for _ in 1:i])
            * poly([β for _ in 1:(n - i)])
            * (-1)^(n - i)
            / (β - α)^n)
end
get_p(B::Bernstein) = get_p(B.n, B.i, B.α, B.β)

function ∫(p::Poly, α=NaN, β=NaN)
    integrated = polyint(p)
    if isequal(α, NaN) || isequal(β, NaN)
        return integrated
    else
        return integrated(β) - integrated(α)
    end
end

function inner_product(p::Poly, q::Poly, α=0., β=1.)
    @assert β >= α "Invalid interval"
    ∫(p * q, α, β)
end

function inner_product(B::Bernstein, Q::Bernstein)
    @assert B.α >= Q.α "Invalid interval"
    @assert B.β >= Q.β "Invalid interval"
    m, i, n, j, α, β = B.n, B.i, Q.n, Q.i, B.α, B.β
    (β - α)binomial(m, i)binomial(n, j) / ((m + n + 1)binomial(m + n, i + j))
end

# QUESTION: Is it possible to automatically account for commutativity?
(inner_product(p::Poly, B::Bernstein)
 = inner_product(B::Bernstein, p::Poly)
 = inner_product(B.p, p::Poly))

function Base.LinAlg.norm(q::Poly, α=0., β=1.)
    @assert β >= α "Invalid interval"
    √inner_product(q, q) / (β - α)
end

function Base.LinAlg.norm(B::Bernstein)
    √inner_product(B, B) / (B.β - B.α)
end

B = Bernstein(3, 2)
Q = Bernstein(4, 3)

inner_product(B, Q)
norm(B)

Answer 1

The below code uses Julia 1.0. Full code as a Pkg3 project is available here.

Preliminaries

There's enough code for it to get its own module. This also makes testing easier.

module Bernstein

using Polynomials

import Base: convert, promote_rule, show
import LinearAlgebra: dot, norm
import Polynomials: poly

Since we're in a module, we must define the exported things. Since most of the things is just overloaded methods, only the type remains.

export BernsteinPoly

You do this check so often, we can factor it out. Alternative: using a separate type for intervals.

macro checkinterval(α, β)
    :(@assert $(esc(α)) <= $(esc(β)) "Invalid interval")
end

Parametrization by a constrained type is better than fields of an abstract type. But in the case of n and i, concrete Ints make more sense, since they represent natural numbers always (this makes conversions below easier). And I'm calling it BernsteinPoly, since it is a variant of Poly.

struct BernsteinPoly{T<:Number}
    n::Int
    i::Int
    α::T
    β::T
end

BernsteinPoly(n, i) = BernsteinPoly(n, i, 0.0, 1.0)

Just implementing show using print is the recommended way for custom pretty-printing. Here, we can reuse Polynomials.printpoly for nicer and more consistent formatting:

function show(io::IO, b::BernsteinPoly)
    print(io, "BernsteinPoly(")
    Polynomials.printpoly(io, poly(b))
    print(io, ")")
end

Conversions

get_p was essentially the conversion from BernsteinPoly to Poly. We can replace this by the implementation of proper convert methods. The default parameters are unnecessary, since already occuring defaulted in the BernsteinPoly constructor.

function convert(::Type{Poly{S}}, b::BernsteinPoly) where {S<:Number}
    n, i, α, β = b.n, b.i, convert(S, b.α), convert(S, b.β)
    @checkinterval α β

    return (binomial(n, i)
            * poly(fill(α, i))     # `fill` instead of list comprehension
            * poly(fill(β, n - i))
            * (-1)^(n - i)
            / (β - α)^n)
end

This allows you to just say convert(Poly, b), automatically reusing the inner type of b

convert(::Type{Poly}, b::BernsteinPoly{T}) where {T<:Number} =
    convert(Poly{T}, b)

While we're at it: conversion between different BernsteinPoly values

convert(::Type{BernsteinPoly{S}}, b::BernsteinPoly) where {S<:Number} =
    BernsteinPoly(b.n, b.i, convert(S, b.α), convert(S, b.β))

If we're handling different representations, we sometimes need to determine the "most general" form, which is called promotion:

promote_rule(a::Type{BernsteinPoly{S}}, b::Type{Poly{T}}) where {S<:Number, T<:Number} =
    Poly{promote_type(S, T)}

promote_rule(::Type{BernsteinPoly{S}}, ::Type{BernsteinPoly{T}}) where {S<:Number, T<:Number} = 
    BernsteinPoly{promote_type(S, T)}

Also add a method to the poly smart constuctor, which is now trivial:

poly(b::BernsteinPoly) = convert(Poly, b)

Linear Algebra

Now to the linear algebra part. We could add methods to dot and norm from LinearAlgebra:

function dot(p::Poly{T}, q::Poly{T}, α = zero(T), β = one(T)) where {T<:Number}
    @checkinterval α β
    polyint(p * q, α, β)
end

And the norm induced by that inner product:

function norm(q::Poly{T}, α = zero(T), β = one(T)) where {T<:Number}
    @checkinterval α β
    √dot(q, q, α, β) / (β - α)
end

But such "overloads from outside" are frowned upon. As of writing this, there's a norm method in Polynomials, but not dot. There's an issue about that; basically, dot is not unique, hence we need to specify the intervals each time here.

For Bernstein polynomials, on the other hand, the inner product is defined uniquely, if I understood correctly.

function dot(b::BernsteinPoly{T}, q::BernsteinPoly{T}) where {T<:Number}
    @checkinterval b.α q.α
    @checkinterval b.β q.β
    m, i, n, j, α, β = b.n, b.i, q.n, q.i, b.α, b.β
    (β - α) * binomial(m, i) * binomial(n, j) / ((m + n + 1) * binomial(m + n, i + j))
end

And the induced norm, as before.

function norm(b::BernsteinPoly)
    √dot(b, b) / (b.β - b.α)
end

As the last remaining question, how to do cross-type inner products. If dot is defined for Poly as written above, we can use dot(promote(p, q)..., α, β). But it's difficult to get that working as a method. I tried

dot(p::Union{Poly{T}, BernsteinPoly{T}}, q::Union{Poly{T}, BernsteinPoly{T}},
     α = zero(T), β = one(T)) where {T<:Number} =
     dot(promote(p, q)..., α, β)

but that doesn't work. It's probably not recommended to do that anyway. Use promote explicitely where necessary.

Example output

julia> B = BernsteinPoly(3, 2)
BernsteinPoly(3.0*x^2 - 3.0*x^3)

julia> Q = BernsteinPoly(4, 3)
BernsteinPoly(4.0*x^3 - 4.0*x^4)

julia> using LinearAlgebra

julia> dot(B, Q)
0.07142857142857142

julia> norm(B)
0.29277002188455997


end # module