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> using FourierSeriesEvaluators
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 FourierSeriesEvaluators, 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> using FourierSeriesEvaluators
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, FourierSeriesEvaluators
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.FourierSeries — Type
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.

HermitianFourierSeries

In physics problems with periodicity, often is it possible to Fourier interpolate Hermitian operators, such as in Wannier interpolation. The HermitianFourierSeries optimizes the interpolation by using the Hermitian symmetry of the coefficients of such an interpolant, where $C_{i} = C_{-i}^{\dagger}$, resulting in a roughly 2x speedup and saving half the memory.

julia> using FourierSeriesEvaluators, OffsetArrays, StaticArrays
julia> c = OffsetVector(rand(SMatrix{2,2,ComplexF64,4}, 7), -3:3);
julia> c = c + reverse(adjoint.(c)); # symmetrize the coefficients
julia> f = FourierSeries(c, period=1);
julia> hf = HermitianFourierSeries(f)1-dimensional and (1,)-periodic HermitianFourierSeries

As shown above, a HermitianFourierSeries can be constructed from an existing FourierSeries whose coefficients satisfy the symmetry property.

julia> x = rand();
julia> hf(x)2×2 StaticArrays.SHermitianCompact{2, ComplexF64, 3} with indices SOneTo(2)×SOneTo(2):
 2.08197+0.0im         1.90007-0.859546im
 1.90007+0.859546im  -0.650441+0.0im

Observe that the output type is SHermitianCompact, which is a special type that exploits the symmetry.

FourierSeriesEvaluators.HermitianFourierSeries — Type
HermitianFourierSeries(f::FourierSeries)

Construct an efficient evaluator for an operator-valued FourierSeries whose coefficients have Hermitian symmetry, i.e. $C_{i} = C_{-i}^{\dagger}$. This evaluator guarantees that its output ishermitian, and is roughly 2x faster than a FourierSeries and uses half of the memory.

Inplace series are not currently supported

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> using FourierSeriesEvaluators
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.ManyFourierSeries — Type
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> using FourierSeriesEvaluators
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.DerivativeSeries — Type
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> using FourierSeriesEvaluators
julia> s = FourierSeries(rand(17,17,17), period=1);
julia> x = ntuple(_->rand(), ndims(s));
julia> ws = FourierSeriesEvaluators.workspace_allocate(s, Tuple(x));
julia> @time s(x);  0.007170 seconds (13.94 k allocations: 1009.523 KiB, 99.13% compilation time)
julia> @time ws(x);  0.008379 seconds (28.94 k allocations: 1.989 MiB, 99.26% compilation time)

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.workspace_allocate — Method
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.