problems getting boundary conditions to work...

[zahachtah]

I am probably missing something obvious but trying to keep close to the documentation:

N=20
M=10
V=3
z=sort(randn(N).*V .+M)
u=gaussian(z,10,3)
B=interpolate(z,u, BSplineOrder(4))
bB = RecombinedBSplineBasis(Derivative(0), B)

gives error:

MethodError: no method matching BSplineKit.BSplines.BSplineBasis(::BSplineKit.SplineInterpolations.SplineInterpolation{BSplineKit.Splines.Spline{Float64, BSplineKit.BSplines.BSplineBasis{4, Float64, Vector{Float64}}, Vector{Float64}}, LinearAlgebra.LU{Float64, BSplineKit.Collocation.CollocationMatrix{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, Vector{Float64}})
Closest candidates are:
BSplineKit.BSplines.BSplineBasis(!Matched::Integer, !Matched::Any...; kwargs...) at /Users/raukhur/.julia/packages/BSplineKit/9U37V/src/BSplines/basis.jl:76
BSplineKit.BSplines.BSplineBasis(!Matched::BSplineKit.BSplines.BSplineOrder{k}, !Matched::AbstractVector{T}; augment) where {k, T} at /Users/raukhur/.julia/packages/BSplineKit/9U37V/src/BSplines/basis.jl:54
var"#RecombinedBSplineBasis#9"(::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, ::Type{BSplineKit.Recombinations.RecombinedBSplineBasis}, ::Tuple{BSplineKit.DifferentialOps.Derivative{1}}, ::BSplineKit.SplineInterpolations.SplineInterpolation{BSplineKit.Splines.Spline{Float64, BSplineKit.BSplines.BSplineBasis{4, Float64, Vector{Float64}}, Vector{Float64}}, LinearAlgebra.LU{Float64, BSplineKit.Collocation.CollocationMatrix{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, Vector{Float64}})@bases.jl:247
BSplineKit.Recombinations.RecombinedBSplineBasis(::Tuple{BSplineKit.DifferentialOps.Derivative{1}}, ::BSplineKit.SplineInterpolations.SplineInterpolation{BSplineKit.Splines.Spline{Float64, BSplineKit.BSplines.BSplineBasis{4, Float64, Vector{Float64}}, Vector{Float64}}, LinearAlgebra.LU{Float64, BSplineKit.Collocation.CollocationMatrix{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, Vector{Float64}})@bases.jl:247
BSplineKit.Recombinations.RecombinedBSplineBasis(::BSplineKit.DifferentialOps.Derivative{1}, ::BSplineKit.SplineInterpolations.SplineInterpolation{BSplineKit.Splines.Spline{Float64, BSplineKit.BSplines.BSplineBasis{4, Float64, Vector{Float64}}, Vector{Float64}}, LinearAlgebra.LU{Float64, BSplineKit.Collocation.CollocationMatrix{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, Vector{Float64}})@bases.jl:212
top-level scope@Local: 8[inlined]

[jipolanco]

Hi, I'm not sure I understand what you're trying to do. Could you provide some more details?

• If I understand correctly, you're trying to construct an interpolator that satisfies homogeneous Dirichlet boundary conditions?
• Where is gaussian defined, and what does it return?

In any case, your code fails because your B is not a B-spline basis, it is a spline interpolation. The following would work, although I'm not sure if that's what you're looking for:

bB = RecombinedBSplineBasis(Derivative(0), basis(B))

[zahachtah]

Here's what I am trying to do:

I get a distribution of species abundances, u, along a trait gradient, z. z is somewhat randomly distributed. I want to get a function zi such that I can get its approximate estimated abundance at any point, not just at each value in the array z. the gaussian function is just a way to create one distribution, but in my simulations I will get many different kinds of distributions.

 gaussian(z,Z,V)=exp.(-((z.-Z).^2 ./V))

What I would preferably like to get is a function approximation that asymptotically goes towards zero (if the endpoints abundances are not zero), or stays at zero if the end points abundances are zero. I'd like to avoid the bumps at both ends like in this figure:

I have to admit I am no engineer and am not familiar with the boundary nomenclature.

current code:

        gaussian(z,Z,V)=exp.(-((z.-Z).^2 ./V))

    N=20
    M=10
    V=3
    z=sort(randn(N).*V .+M)
    u=gaussian(z,10,3)
    zi=interpolate(z,u, BSplineOrder(4))

    scatter(z,u) #using makie
    zz=range(0.0,stop=25.0,length=100)
    lines!(zz,zi.(zz))
    current_figure()

I tried to use

        bB = RecombinedBSplineBasis(Derivative(0), basis(zi))
    scatter(z,u) #using makie
    zz=range(0.0,stop=25.0,length=100)
    lines!(zz,bB.(zz))
    current_figure()

but then when trying to plot it I get error that bB is not callable. Sorry, I am sure its somewhere in the documentation but after trying and reading I have still not solved it. thanks for any hint

[jipolanco]

Thanks, now I understand better what you're trying to do.

Unfortunately there's no easy way to do this in BSplineKit, but here is an idea that may help. One thing you could try is to avoid using the high-level interpolate function (which only works for very simple cases), and to start by building a (recombined) B-spline basis with boundary conditions and knots that you specify.

One way to do this is to define your B-spline breakpoints by extending your data locations z with an extra point on the left and on the right. Then, the boundary conditions will be enforced on these two points. In your example, you could do something like this:

offset = sqrt(V)
breaks = vcat(z[begin] - offset, z, z[end] + offset)

The offset is of course a free parameter that one could play with. Here, in your Gaussian example, it made sense to me to make it equal (or proportional) to the standard deviation of the distribution, but other values may make more sense in general. Also, note that the breakpoints don't really need to match the data points, so there are other (infinite) possible choices of breakpoints.

Then, you would define a B-spline basis and then a recombined B-spline basis satisfying certain boundary conditions (BCs):

B = BSplineBasis(BSplineOrder(4), breaks)
R = RecombinedBSplineBasis((Derivative(0), Derivative(1)), B)

Note that here I'm imposing both the value and the derivative of the function to be zero, which again can make sense for the asymptotic behaviour of a distribution. But actually, in this example, setting two BCs per boundary is actually necessary to reduce the number of degrees of freedom by 4, as we want the number of unknown basis coefficients to be equal to the number of data points. In other words:

length(z)  # = 20
length(B)  # = 24 (too many degrees of freedom)
length(R)  # = 20 (this is what we want)

Then, you can manually perform an interpolation on the recombined basis by building a collocation matrix and solving a linear equation to get the spline coefficients (see this post of mine on the Julia Discourse for more details):

C = collocation_matrix(R, z)  # this is a banded 20×20 matrix
coefs = lu!(C) \ u            # these are the spline coefficients (length 20)
S = Spline(R, coefs) 

---

Here is the full example and the resulting plot:

using BSplineKit
using GLMakie
using LinearAlgebra
using Random

gaussian(z,Z,V)=exp.(-((z.-Z).^2 ./V))

N = 20
M = 10
V = 3
rng = MersenneTwister(48)
z = sort(randn(rng, N) .* V .+ M)

u=gaussian(z,10,3)
offset = sqrt(V)
breaks = vcat(z[begin] - offset, z, z[end] + offset)
B = BSplineBasis(BSplineOrder(4), breaks)
R = RecombinedBSplineBasis((Derivative(0), Derivative(1), ), B)
@show length(B), length(R), length(z)

C = collocation_matrix(R, z)
coefs = lu!(C) \ u
S = Spline(R, coefs)

fig = Figure()
ax = Axis(fig[1, 1])
scatter!(ax, z, u)
scatter!(ax, [breaks[begin], breaks[end]], [0, 0])
lines!(ax, zz, S.(zz))
fig

The orange points are the breakpoints that were added, where both BCs are satisfied.

Note that this approach is not perfect if you're trying to approximate a probability distribution. In particular, the resulting function may have negative values. For instance, in the above example, if one evaluates S(3), one gets something close to -1.2e-6. One way to avoid this would be to make sure that all spline coefficients (the coefs vector) are non-negative, which would imply that the spline itself is non-negative. This may be done by, instead of solving the linear system coefs = lu!(C) \ u, solving a constrained optimisation problem (with the constraint coefs .>= 0). Or simpler, if you don't need an optimal approximation, you could do clamp!.(coefs, 0, Inf) to enforce this condition.

[zahachtah]

This is amazing! Many thanks for this thorough answer! I realise I would not have been able to solve this myself, but this solution may indeed work nicely for my purpose! The approximation does not have to be optimal, so the clamp solution could indeed work.

Again, many thanks for taking time to explain and show how to do this!

/Jon