Cfrtjryk

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)

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