# Bernstein polynomial type for quadratic clipping

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)  ## 1 Answer 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