Examples

The following are several examples of how to use the Fourier series evaluators exported by FourierSeriesEvaluators together with their documentation.

FourierSeries

The FourierSeries constructor accepts an array of coefficients and a mandatory keyword argument, period. To create a 3-dimensional series with random coefficients we can use the following

julia> f = FourierSeries(rand(4,4,4), period=1.0)
4×4×4 and (1.0, 1.0, 1.0)-periodic FourierSeries with Float64 coefficients, (0, 0, 0) derivative, (0, 0, 0) offset

We can then evaluate f at a random point using its function-like interface f(rand(3)).

To control the indices of the coefficients, which determine the phase factors according to a formula given later, we can either use arrays with offset indices, or provide the offsets ourselves

julia> using OffsetArrays

julia> c = OffsetVector(rand(7), -3:3);

julia> x = rand();

julia> FourierSeries(c, period=1)(x) == FourierSeries(parent(c), period=1, offset=-4)(x)
true

If we want to take the derivative of a Fourier series such as the derivative of sine, which is cosine, we can first construct and validate our version of sin

julia> sine = FourierSeries([-1, 0, 1]/2im, period=2pi, offset=-2)
3-element and (6.283185307179586,)-periodic FourierSeries with ComplexF64 coefficients, (0,) derivative, (-2,) offset

julia> x = rand();

julia> sine(x) ≈ sin(x)
true

and then reuse the arguments to the series and provide the order of derivative with the deriv keyword, which is typically integer but can even be fractional

julia> cosine = FourierSeries([-1, 0, 1]/2im, period=2pi, offset=-2, deriv=1)
3-element and (6.283185307179586,)-periodic FourierSeries with ComplexF64 coefficients, (1,) derivative, (-2,) offset

julia> cosine(x) ≈ cos(x)
true

For multidimensional series, a scalar value of deriv means it is applied to derivatives of all variables, whereas a tuple or vector of orders will select the order of derivative applied to each variable. For automatically generating all orders of derivatives see the section on DerivativeSeries.

Lastly, if we wish to evaluate multiple Fourier series at the same time, we may either use array-valued coefficients, such as with StaticArrays.jl, or for very large arrays we may pass a second positional argument to FourierSeries to provide the number of variables in the FourierSeries corresponding to the outermost axes of the coefficient array. Below we have an example constructing equivalent evaluators using both styles.

julia> c = rand(5,7);

julia> using StaticArrays

julia> sc = reinterpret(reshape, SVector{5,eltype(c)}, c);

julia> x = rand();

julia> FourierSeries(sc, period=1.0)(x) ≈ FourierSeries(c, 1, period=1.0)(x)
true

When using the form with the positional argument, note that inplace array operations are used to avoid allocations, so be careful not to overwrite the result if reusing it.

FourierSeriesEvaluators.FourierSeriesType
FourierSeries(coeffs::AbstractArray, [N]; period, offset=0, deriv=0)

Construct a Fourier series whose coefficients are given by the coefficient array array coeffs whose elements should support addition and scalar multiplication, This object represents the Fourier series

$$$f(\vec{x}) = \sum_{\vec{n} \in \mathcal I} C_{\vec{n}} \left( \prod_{i} \left( 2\pi f_{i} (n_{i} + o_{i}) \sqrt{-1} \right)^{a_{i}} \exp\left(2\pi f_{i} x_{i} (n_{i} + o_{i}) \sqrt{-1} \right) \right)$$$

Here, the indices $\mathcal I$ are the CartesianIndices of coeffs, $f_{i} = 1/t_{i}$ is the frequency/inverse period, $a_{i}$ is the order of derivative, and $o_{i}$ is the index offset. Also, the keywords, which can either be a single value applied to all dimensions or a tuple/vector describing each dimension mean

• period: The periodicity of the Fourier series, which must be a real number
• offset: An offset in the phase index, which must be integer
• deriv: The degree of differentiation, implemented as a Fourier multiplier. Can be any number a such that x^a is well-defined, or a Val(a). a::Integer performs best.

If inplace evaluation or evaluation of multiple series is desired, the optional argument N when set fixes the number of variables of the Fourier series, which may be less than or equal to ndims(coeffs), and the series is evaluated inplace, returning the innermost, continguous ndims(coeffs)-N axes evaluated at the variables corresponding to the remaining outer axes.

ManyFourierSeries

A third way to evaluate multiple series at the same time, possibly with different types of coefficients, is with a ManyFourierSeries, which behaves like a tuple of FourierSeries. We can construct one from multiple series and a period (which the element series are rescaled by so that all periods match).

julia> f1 = FourierSeries(rand(3,3), period=3);

julia> f2 = FourierSeries(rand(3,3), period=2);

julia> x = rand(2);

julia> ms = ManyFourierSeries(f1, f2, period=1)
2-dimensional and (1, 1)-periodic ManyFourierSeries with 2 series

julia> ms(x)[1] == f1(3x)
true

julia> ms(x)[2] == f2(2x)
true

There are situations in which inplace FourierSeries may be more efficient than ManyFourierSeries, since the former calculates phase factors fewer times than the later.

FourierSeriesEvaluators.ManyFourierSeriesType
ManyFourierSeries(fs::AbstractFourierSeries{N,T,iip}...; period) where {N,T,iip}

Represents a tuple of Fourier series of the same dimension and contracts them all simultaneously. All the series are required to be either inplace or not.

DerivativeSeries

To systematically compute all orders of derivatives of a Fourier series, we can wrap a series with a DerivativeSeries{O}, where O::Integer specifies the order of derivative to evaluate athe series and take all of its derivatives up to and including degree O.

julia> sine = FourierSeries([-1, 0, 1]/2im, period=2pi, offset=-2)
3-element and (6.283185307179586,)-periodic FourierSeries with ComplexF64 coefficients, (0,) derivative, (-2,) offset

julia> ds = DerivativeSeries{1}(sine)
1-dimensional and (6.283185307179586,)-periodic DerivativeSeries of order 1

julia> ds(0)
(0.0 + 0.0im, (1.0 + 0.0im,))

We can use this to show that the fourth derivative of sine is itself

julia> d4s = DerivativeSeries{4}(sine)
1-dimensional and (6.283185307179586,)-periodic DerivativeSeries of order 4

julia> d4s(pi/3)
(0.8660254037844386 + 0.0im, (0.5000000000000001 + 0.0im,), ((-0.8660254037844386 + 0.0im,),), (((-0.5000000000000001 + 0.0im,),),), ((((0.8660254037844386 + 0.0im,),),),))

When applied to multidimensional series, all mixed partial derivatives are computed exactly once and their layout in the tuple of results is explained below.

FourierSeriesEvaluators.DerivativeSeriesType
DerivativeSeries{O}(f::AbstractFourierSeries)

Construct an evaluator of a Fourier series and all of its derivatives up to order O, which must be a positive integer. O=1 gives the gradient, O=2 gives the Hessian, and so on. The derivatives are returned in order as a tuple (f(x), df(x), d2f(x), ..., dOf(x)) where the entry of order O is given by:

• O=0: f
• O=1: (dfdx1, ..., dfdxN)
• O=2: ((d2fdx1dx1, ..., d2fdx1dxN), ..., (d2fdxNdxN,))
• O=3: (((d3fdx1dx1dx1, ..., d3fdx1dx1dxN), ..., (d3fdx1dxNdxN,)), ..., ((d3fdxNdxNdxN,),))

and so on. The fewest number of contractions are made to compute all derivatives. As can be seen from the pattern above, the O-th derivative with partial derivatives i = [a_1 ≤ ... ≤ a_N] is stored in ds(x)[O+1][i[1]][i[2]]...[i[N]]. These indices are given by the simplical generalization of triangular numbers. For examples of how to index into the solution see the unit tests.

For this routine to work, f must implement nextderivative.

FourierWorkspace

By default, evaluating a Fourier series using the function-like interface allocates several intermediate arrays in order to achieve the fastest-possible evaluation with as few as possible calculations of phase factors. These arrays can be preallocated in a FourierWorkspace

julia> s = FourierSeries(rand(17,17,17), period=1)
17×17×17 and (1, 1, 1)-periodic FourierSeries with Float64 coefficients, (0, 0, 0) derivative, (0, 0, 0) offset

julia> ws = FourierSeriesEvaluators.workspace_allocate(s, Tuple(x));

julia> @time s(x);
0.000044 seconds (3 allocations: 5.047 KiB)

julia> @time ws(x);
0.000063 seconds (1 allocation: 32 bytes)

We can also allocate nested workspaces that can be used independently to evaluate the same series at many points in a hierarchical or product grid in parallel workloads.

julia> ws3 = FourierSeriesEvaluators.workspace_allocate(s, Tuple(x), (2,2,2));

julia> ws2 = FourierSeriesEvaluators.workspace_contract!(ws3, x[3], 1);

julia> ws1 = FourierSeriesEvaluators.workspace_contract!(ws2, x[2], 2);

julia> FourierSeriesEvaluators.workspace_evaluate!(ws1, x[1], 1) == s(x)
true

Note that the 3rd argument of workspace_allocate, workspace_contract!, and workspace_evaluate! either specifies the number of nested workspaces to create or selects the workspace to use.

FourierSeriesEvaluators.FourierWorkspaceType
FourierWorkspace

A workspace for storing a Fourier series and the intermediate arrays used for evaluations. Given a ws::FourierWorkspace, you can evaluate it at a point x with ws(x).

All functionality is implemented by the workspace_allocate, workspace_contract!, and workspace_evaluate! routines, which allow multiple caches within each dimension of evaluation to enable parallel workloads.

FourierSeriesEvaluators.workspace_allocateMethod
workspace_allocate(s::AbstractFourierSeries{N}, x::NTuple{N}, [len::NTuple{N}=ntuple(one,N)])

Allocates a FourierWorkspace for the Fourier series s that can be used to evaluate the series multiple times without allocating on-the-fly. The len argument can indicate how many copies of workspace should be made for each variable for downstream use in parallel workloads.

The workspace is constructed recursively starting from the outer dimension and moving towards the inner dimension so as to access memory contiguously. Thus, the outer dimension has len[N] workspace copies and each of these has len[N-1] workspaces for the next variable. In total there are prod(len) leaf-level caches to use for parallel workloads.

FourierSeriesEvaluators.workspace_contract!Function
workspace_contract!(ws, x, [i=1])

Returns a workspace with the series contracted at variable x in the outer dimension. The index i selects which workspace in the cache to assign the new data.