# 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.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.

`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.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> 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`

.

`FourierSeriesEvaluators.JacobianSeries`

— Type`JacobianSeries`

Alias for a `DerivativeSeries`

of order 1

`FourierSeriesEvaluators.HessianSeries`

— Type`HessianSeries`

Alias for a `DerivativeSeries`

of order 2

`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.FourierWorkspace`

— Type`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_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.

`FourierSeriesEvaluators.workspace_evaluate!`

— Function`workspace_evaluate!(ws, x, [i=1])`

Return the 1-d series evaluated at the variable `x`

, using cache sector `i`

.