Sampling points in 2D

[kuanxu]

For Framefun approximation in 2D:

1. Are the default sampling points equispaced points in a cartesian grid? Taking the 'disk with a hole' in README.md as an example, how can I know/retrieve the actual sampling points?

2. Is it possible to construct a 2D Framefun approximant by giving the ctor a set of sampling points and discrete data of a same size (as function values at these points)?

3. Is the construction adaptive in 2D?

[daanhb]

1. You can ask any dictionary or any platform what its sampling_grid is. This function takes an optional oversamplingfactor argument. For example:

   
   julia> using BasisFunctions, FrameFun, IntervalSets

julia> frame = ExtensionFrame(0..1.0, ChebyshevT(10)) Dictionary 𝔼(C)

𝔼 : Extension frame, from 0.0..1.0 to -1.0..1.0 (Chebyshev) C : Chebyshev polynomials (first kind) ↳ length = 10 ↳ Float64 -> Float64 ↳ support = -1.0..1.0 (Chebyshev)

julia> sampling_grid(frame) 20-element GridArrays.MaskedGrid{ChebyshevTNodes{Float64},BitArray{1},Int64,Interval{:closed,:closed,Float64},Float64}: 0.039259815759068326 0.11753739745783758 0.1950903220161282 0.2714404498650744 0.3461170570774927 0.4186597375374278 ...

This code creates a frame by restricting Chebyshev polynomials of the first kind from `[-1,1]` to `[0,1]`. Its sampling points are a subset of ChebyshevTNodes: Chebyshev points are generated on `[-1,1]` with a length such that their restriction to `[0,1]` is twice the number of degrees of freedom in the frame (i.e. the oversamplingfactor here is 2).

2. We don't have a convenient constructor for this yet, it is possible if you dig a little deeper. You can easily construct the linear system that is needed to find the approximation, but you can not easily provide the constructor with sampled data. It expects a function.
Any suggestion for an interface? This should be easy to add, and it would come in handy for some other things we do.

3. The `Fun` method is adaptive if you supply it a platform, rather than a fixed basis. This works in 2D, yes. A platform is essentially an object that maps a parameter (such as N) to a basis (such as Chebyshev polynomials of degree N). The platform can also overrule any defaults (such as which sampling grids are used to compute an approximation). It is very flexible, but I'm afraid largely undocumented. There are some examples in the `Adaptivity` notebook in the examples folder.

The easiest way to get a platform is to construct any basis and then use `ModelPlatform(basis)`. The platform will then attempt to resize the basis to other lengths in the context of adaptivity and it will (try to) use sensible defaults.

[daanhb]

Also, apologies for the late reply! I'm not actively monitoring the issues here. Please let me know if there is something you try to do that is not working, we'll make it work. Ask here or drop me an e-mail ;-)

[kuanxu]

@daanhb Thanks for the hints! They're most helpful! In fact, one of my undergrad students is now working on a project of computational optics for which we are looking at an example problem with the following code:

using BasisFunctions, Plots, DomainSets, FrameFun, IntervalSets
gr();
using StaticArrays
C = disk(1.0)
B = Fourier(31,-1.3,1.3) ⊗ Fourier(31,-1.3,1.3)
f = (x,y)->exp(x+y)
F = Fun(f,B,C)

P = heatmap(F,plot_ext=true,layout=2,aspect_ratio=1)
plot!(F,f,plot_ext=true,subplot=2,aspect_ratio=1)

1. sampling_grid worked for us as we tried to plot the sampling points in the disk by plot(sampling_grid(F.expansion.dictionary)) (see the figure below). Thanks!

2. In fact, we do rely on 2D evenly-spaced grids in this project. But we don't have a function handle to pass to the constructor like in the code above. The only thing we have is discrete sampling data (evenly-spaced data inside the circle as shown in the figure above) for the framefun constructor since the data were obtained from experiments. It'll be great if the constructor allows us to pass a set of data instead of a function handle. Is there anyone in the framefun team can make this function possible in the near future?

3. Cool. Except the adaptivity which is not quite documented, can we say the the underlying approximation algorithms/methods in 2D framefun are pretty much based on the paper Function Approximation on Arbitrary Domains Using Fourier Extension Frames?

4. In the code above, we specify the number of freedom in each direction as 31. However, we found the actual discretization uses 65 points instead and this change happens in function deduce_samplingparameter as shown below in a screenshot. Why does this happen?

5. We know we can evaluate the framefun approximant F at a point inside the circle like F(0.5, 0.3). How can we do the evaluation of F at multiple points in a single batch? Can you give us a hint? We try quite a few possibilities, but failed.

Sorry for this long list of questions! Thanks a lot!

[daanhb]

1. Good to hear. Note that the sampling grid in this case is a MaskedGrid, which knows that it is a restriction of a cartesian product periodic equispaced grid (FourierGrid). You can do sampling_grid(dictionary(F)), but the sampling_grid function actually accepts the same arguments as the Fun constructor without the f. This way you can know for sure in which grid f will be sampled if you invoke Fun. This may differ from asking the Fun afterwards, because the Fun does not store the actual grid that was used during construction. You can also override the grid to anything you like with an optional grid argument.

In your example:

julia> g = sampling_grid(B, C)
1976-element GridArrays.MaskedGrid{...

julia> size(g)
(1976,)

julia> supergrid(g)
65×65 ProductGrid{Tuple{MappedGrid{FourierGrid{Float64},AffineMap{Float64,Float64,Float64},Float64,1},MappedGrid{FourierGrid{Float64},AffineMap{Float64,Float64,Float64},Float64,1}},SArray{Tuple{2},Float64,1,2},2}
...

The sampling grid acts like a vector with length 1976, but it is a subset of a 65 by 65 product of periodic equispaced grids (FourierGrid).

Everything is configurable, but sensible defaults are chosen normally so you don't need to be aware of how things can be changed. Still, if you want to provide your own grid, you can use the GridStyle as follows:

julia> b = ChebyshevT(20)

julia> sampling_grid(b)
20-element ChebyshevTNodes{Float64}:
 -0.996917333733128
 -0.9723699203976766
 -0.9238795325112867
...

julia> sampling_grid(b, samplingstyle=GridStyle(), grid=EquispacedGrid(1000, -1.0..1.0))
1000-element EquispacedGrid{Float64}:
 -1.0
 -0.997997997997998
 -0.9959959959959961
 -0.9939939939939939
...

This means that if you invoke Fun(exp, b, samplingstyle=GridStyle(), grid=EquispacedGrid(1000, -1.0..1.0)), you won't be using the 20 point Chebyshev nodes to approximate the exponential function, but f will be sampled in 1000 equispaced points on [-1,1] and a least squares approximation is computed.

This may come in handy if you only know samples of a function. You can force the Fun constructor to only use your grid.

[daanhb]

2. What we would need to provide is a way to pass function values rather than a function to Fun. I think this will most easily be achieved by accepting an array as the f argument to Fun. We never assume anything about f other than it can be called and evaluated. We could easily dispatch on arrays, and if the size of the array matches the size of the sampling grid we can safely assume that the user intended to provide data rather than a function. This would be useful to us in another ongoing project, so we'll get to this soon.

[daanhb]

3. Yes, for Fourier extension that reference is the latest one. Next week we will put a new paper on the arxiv about the work of @vincentcp on spline-based extension approximations. It is titled "Efficient function approximation on general bounded domains using splines on a cartesian grid". This paper is done.

We are also finalizing a paper generalizing the Fourier extension algorithm to what we've called the AZ algorithm, which works for much more general frames. The paper is called "The AZ algorithm for least squares systems with a known incomplete generalized inverse", co-authored by Vincent Coppé, Roel Matthysen, Marcus Webb and myself. This will supersede the previous Fourier extension papers. I'll e-mail you a draft. This paper is nearly done.

[daanhb]

4. This is a good question. By default, the FrameFun code will aim to compute an approximation based on oversampling by a factor of 2. You typically supply a number of degrees of freedom N (by providing a fixed dictionary like in your example) so we know how many sampling points we need: we aim by default for M = 2N. However, it is not easy to generate a grid with exactly that many sampling points on a disk.

In contrast, it is easy to generate cartesian product grids on the bounding box. Say you have an L-by-L grid of points on the box. The restriction of that grid to the disk results in a subset of h(L) points, where h is some non-linear function. In the code, we call L the sampling parameter. The routine deduce_sampling_parameter solves the problem h(L)=M for L. Internally, we then store the parameter L. That means we can easily generate the grid on the fly: construct the L-by-L grid on the bounding box and keep only those points that belong to the domain at hand.

If you supply your own grid, then you won't need this. You can also give the value for L to the Fun constructor and avoid the optimization problem:

julia> b = ExtensionFrame(0..1.0, ChebyshevT(20))
Dictionary 𝔼(C)

𝔼   :   Extension frame, from 0.0..1.0 to -1.0..1.0 (Chebyshev)
C   :   Chebyshev polynomials (first kind)
            ↳ length = 20
            ↳ Float64 -> Float64
            ↳ support = -1.0..1.0 (Chebyshev)

julia> length(supergrid(sampling_grid(b)))
80

julia> sampling_grid(b, verbose=true)
Sampling parameter: oversamplingfactor=2, length=20, M=40
Sampling parameter: best match for M = 40 is L = 80
Sampling parameter: using L = 80
Sampling parameter: using L = 80
40-element GridArrays.MaskedGrid{ChebyshevTNodes{Float64},BitArray{1},Int64,Interval{:closed,:closed,Float64},Float64}:
...

julia> sampling_grid(b, verbose=true, L=150)
Sampling parameter: using L = 150
Sampling parameter: using L = 150
75-element GridArrays.MaskedGrid{ChebyshevTNodes{Float64},BitArray{1},Int64,Interval{:closed,:closed,Float64},Float64}:
 0.010471784116245747
 0.03141075907812828
 0.05233595624294384
...

Note the verbose option :-) It comes in handy if you want to figure out what is actually happening! In the example above, the sampling grid is a length M=40 subgrid on 0..1 of L=80 Chebyshev nodes on [-1,1]. We changed L to 150 by passing it to the constructor.

In your original example, the sampling parameter is a tuple:

julia> samplingparameter(B, C, verbose=true)
Sampling parameter: oversamplingfactor=2, length=961, M=1922
Sampling parameter: best match for M = 1922 is L = (65, 65)
Sampling parameter: using L = (65, 65)
(65, 65)

[daanhb]

5. The elegant solution would be to use broadcast: F.([0.5, 0.6, 0.8], [0.7, 0.8, 0.9]).This is currently not always the fastest way.

You can also construct the evaluation operator in some grid:

x = rand(128)
y = rand(128)
A = BasisFunctions.evaluation_operator(dictionary(F), ScatteredGrid(collect(zip(x,y))))

You can multiply A by a vector of coefficients, and the result will be the evaluation of the corresponding expansion in the given grid. I've used a random grid above. If the dictionary recognizes the grid, a more efficient algorithm may be used (such as FFT when evaluating Fourier series in equispaced points).

[kuanxu]

1. Thanks for the elaboration! It's crystal clear to us now.
2. That's a very good news - we can then (brazenly) take advantage of this new feature for our project. :D
3. So the current Framefun code is based on Function Approximation on Arbitrary Domains Using Fourier Extension Frames, correct? Will the code in this repo be replaced by the one corresponding to the two coming papers in the near future? My student needs to make sure that he cites the right paper when gives the readers a pointer to the software.
4. I see. Just to confirm - the oversampling rate 2 is for the problem as a whole, not for only sampling grid in one direction, right?
5. Thanks for the hint! Very rusty on Julia these days - totally forgot the fusion.

[kuanxu]

Forgot to say - we'd like to see the two new papers when they're available! Thanks!

[daanhb]

The new interface where Fun accepts an array of samples instead of a function is now available in master (b8ae094). You simply have to make sure that the array matches the size of the sampling operator.

The oversampling factor is indeed (approximately) 2 as a whole, not per dimension:

julia> length(sampling_grid(ChebyshevT(10)^2, disk())) / 100
2.2

julia> length(sampling_grid(ChebyshevT(10)^2, disk(), oversamplingfactor=3.5)) / 100
3.62

The oversampling factor is only approximately reached, because not all oversampling ratios can be achieved exactly from a family of grids. This depends on the shape of the domain. We err on the side of safety (more oversampling).

The code follows the conventions of the upcoming AZ paper, which you cannot cite yet :-) But citing Roel's paper is fine (thanks!), because for Fourier extensions we essentially use his algorithm.

[kuanxu]

@daanhb You're lightning fast! I'll ask the student to try the master branch and let you know if we have further questions. Thanks!!

It'll take a little while until the student starts to focus on writing the thing up. Will confirm with you then for correction citations. For the time being, we'll take Roel's paper as the reference.

[kuanxu]

@daanhb We have tried the new feature, but failed. Can you show us the correction usage by the example of disk that we used above?

[daanhb]

I had the following in mind:

using BasisFunctions, Plots, DomainSets, FrameFun, IntervalSets
using StaticArrays

C = disk(1.0)
B = Fourier(31,-1.3,1.3) ⊗ Fourier(31,-1.3,1.3)
f = (x,y)->exp(x+y)
F = Fun(f,B,C)

g = sampling_grid(B, C)
fvect(x) = f(x[1],x[2])
Z = fvect.(g)
F2 = Fun(Z, B, C)

abs(F(0.3,0.4)-F2(0.3,0.4))

I obtain the grid that Fun would like to use with sampling_grid. I've had to create a function fvect that takes a vector in order to be able to write Z = fvect.(g), since that automatically creates an array with the same size as g. I then pass this data to Fun. The final line gives me something on the order of 1e-9.

Can you confirm that the above works?

[kuanxu]

We updated all the dependent packages and still failed:

Are we missing anything?

[daanhb]

That's an informative message that rings a bell. I have probably defined the size of a GenericOperator in a local branch that hasn't been pushed yet. I'll look into it.

[daanhb]

That was indeed the case, I had uncommitted changes in my local development branch of BasisFunctions.jl which I did not think of, sorry about that. I'm running the test suite now before I commit them.

[kuanxu]

@daanhb Thanks very much!

[daanhb]

Unfortunately the development branch is in, wel..., development :-) and it breaks the tests currently. I think I've added the relevant definition of the size of a generic operator in one dimension to master. Can you try again?

[kuanxu]

@daanhb Which branch shall we try? Development or master?

[daanhb]

I've pushed to the master branch of BasisFunctions.jl (not FrameFun), so you'll have to update that package.

[kuanxu]

@daanhb Cool. Its works perfectly now! Thanks!!! :D

[daanhb]

Great! Also, please see the very recent paper by @vincentcp on the arXiv. This paper uses B-splines instead of Fourier series on a bounding box, which makes the approximation algorithm much more efficient. The corresponding code is in the package BSplineExtensions.jl.

[kuanxu]

@daanhb Thanks for the pointer to the paper. Looks very nice. Can I say that we need to consult this paper only when we work with non-evenly spaced grids? With evenly-spaced paper, it is Roel's paper that is most relevant, right?

[daanhb]

Roel's paper describes the Fourier extension solver, which is what you are using so it is the relevant paper. The new paper describes approximations based on splines, but still using equispaced grids. We did not consider unequispaced grids yet. The main difference is that B-splines are compactly supported, which makes the algorithm more efficient, but they offer only algebraic convergence rather than spectral convergence.

[kuanxu]

@daanhb Cool. So the spline-based methods are only shipped in BSplineExtensions.jl and will not replace what is currently in Framefun, correct?

[daanhb]

That is correct: the package doesn't replace what is in FrameFun, it just adds functionality and implements the necessary components to make the computation of an approximation more efficient specifically for splines.

[kuanxu]

Awesome! Will play with BSplineExtensions.jl to see the speed.