Hermitian series (#3)

* initial commit

* fast evaluate

* add tests

* documentation

Project.toml

@@ -1,7 +1,16 @@
 name = "FourierSeriesEvaluators"
 uuid = "2a892dea-6eef-4bb5-9d1c-de966c9f6db5"
 authors = ["lxvm <lorenzo@vanmunoz.com>"]
-version = "1.0.1"
+version = "1.1.0"
+
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[weakdeps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[extensions]
+FourierSeriesEvaluatorsStaticArraysExt = "StaticArrays"
 
 [compat]
 Aqua = "0.8"

docs/Project.toml

@@ -1,2 +1,4 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

docs/src/examples.md

@@ -8,46 +8,36 @@ exported by `FourierSeriesEvaluators` together with their documentation.
 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
+```@repl
+using FourierSeriesEvaluators
+f = FourierSeries(rand(4,4,4), period=1.0)
 ```
 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
+```@repl
+using FourierSeriesEvaluators, OffsetArrays
+c = OffsetVector(rand(7), -3:3);
+x = rand();
+FourierSeries(c, period=1)(x) == FourierSeries(parent(c), period=1, offset=-4)(x)
 ```
 
 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
+```@repl sincos
+using FourierSeriesEvaluators
+sine = FourierSeries([-1, 0, 1]/2im, period=2pi, offset=-2)
+x = rand();
+sine(x) ≈ sin(x)
 ```
 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
+```@repl sincos
+cosine = FourierSeries([-1, 0, 1]/2im, period=2pi, offset=-2, deriv=1)
+cosine(x) ≈ cos(x)
 ```
 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
@@ -62,17 +52,12 @@ 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
+```@repl
+c = rand(5,7);
+using StaticArrays, FourierSeriesEvaluators
+sc = reinterpret(reshape, SVector{5,eltype(c)}, c);
+x = rand();
+FourierSeries(sc, period=1.0)(x) ≈ FourierSeries(c, 1, period=1.0)(x)
 ```
 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
@@ -82,27 +67,52 @@ result if reusing it.
 FourierSeriesEvaluators.FourierSeries
 ```
 
-## [`ManyFourierSeries`](@ref)
+## [`HermitianFourierSeries`](@ref)
 
-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).
+In physics problems with periodicity, often is it possible to Fourier
+interpolate Hermitian operators, such as in [`Wannier
+interpolation`](https://en.wikipedia.org/wiki/Wannier_function). 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.
+
+```@repl herm
+using FourierSeriesEvaluators, OffsetArrays, StaticArrays
+c = OffsetVector(rand(SMatrix{2,2,ComplexF64,4}, 7), -3:3);
+c = c + reverse(adjoint.(c)); # symmetrize the coefficients
+f = FourierSeries(c, period=1);
+hf = HermitianFourierSeries(f)
 ```
-julia> f1 = FourierSeries(rand(3,3), period=3);
 
-julia> f2 = FourierSeries(rand(3,3), period=2);
+As shown above, a `HermitianFourierSeries` can be constructed from an existing
+`FourierSeries` whose coefficients satisfy the symmetry property.
 
-julia> x = rand(2);
+```@repl herm
+x = rand();
+hf(x)
+```
 
-julia> ms = ManyFourierSeries(f1, f2, period=1)
-2-dimensional and (1, 1)-periodic ManyFourierSeries with 2 series
+Observe that the output type is `SHermitianCompact`, which is a special type
+that exploits the symmetry.
 
-julia> ms(x)[1] == f1(3x)
-true
+```@docs
+FourierSeriesEvaluators.HermitianFourierSeries
+```
+
+## [`ManyFourierSeries`](@ref)
 
-julia> ms(x)[2] == f2(2x)
-true
+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).
+```@repl
+using FourierSeriesEvaluators
+f1 = FourierSeries(rand(3,3), period=3);
+f2 = FourierSeries(rand(3,3), period=2);
+x = rand(2);
+ms = ManyFourierSeries(f1, f2, period=1)
+ms(x)[1] == f1(3x)
+ms(x)[2] == f2(2x)
 ```
 There are situations in which inplace `FourierSeries` may be more efficient than
 `ManyFourierSeries`, since the former calculates phase factors fewer times than
@@ -118,23 +128,16 @@ 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,))
+```@repl deriv
+using FourierSeriesEvaluators
+sine = FourierSeries([-1, 0, 1]/2im, period=2pi, offset=-2)
+ds = DerivativeSeries{1}(sine)
+ds(0)
 ```
 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,),),),))
+```@repl deriv
+d4s = DerivativeSeries{4}(sine)
+d4s(pi/3)
 ```
 When applied to multidimensional series, all mixed partial derivatives are
 computed exactly once and their layout in the tuple of results is explained
@@ -152,30 +155,22 @@ 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)
+```@repl workspace
+using FourierSeriesEvaluators
+s = FourierSeries(rand(17,17,17), period=1);
+x = ntuple(_->rand(), ndims(s));
+ws = FourierSeriesEvaluators.workspace_allocate(s, Tuple(x));
+@time s(x);
+@time ws(x);
 ```
 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
+```@repl workspace
+ws3 = FourierSeriesEvaluators.workspace_allocate(s, Tuple(x), (2,2,2));
+ws2 = FourierSeriesEvaluators.workspace_contract!(ws3, x[3], 1);
+ws1 = FourierSeriesEvaluators.workspace_contract!(ws2, x[2], 2);
+FourierSeriesEvaluators.workspace_evaluate!(ws1, x[1], 1) == s(x)
 ```
 Note that the 3rd argument of `workspace_allocate`, `workspace_contract!`, and
 `workspace_evaluate!` either specifies the number of nested workspaces to create

docs/src/reference.md

@@ -5,5 +5,6 @@ The following functions are exported by the `FourierSeriesEvaluators` package
 
 ```@autodocs
 Modules = [FourierSeriesEvaluators]
-Pages   = ["fourier_kernel.jl"]
+Order = [:function]
+Pages   = ["fourier_kernel.jl", "hermitian.jl"]
 ```

ext/FourierSeriesEvaluatorsStaticArraysExt.jl

@@ -0,0 +1,7 @@
+module FourierSeriesEvaluatorsStaticArraysExt
+import FourierSeriesEvaluators: _hermitianpart
+using StaticArrays: SHermitianCompact, StaticMatrix
+
+_hermitianpart(A::StaticMatrix) = SHermitianCompact(A)
+
+end

src/FourierSeriesEvaluators.jl

@@ -7,13 +7,15 @@ small compared to the number of evaluation points.
 
 As a first example, we evaluate cosine from its Fourier series coefficients:
 
-    julia> using FourierSeriesEvaluators
+```
+julia> using FourierSeriesEvaluators
 
-    julia> cosine = FourierSeries([0.5, 0.0, 0.5], period=2pi, offset=-2)
-    3-element and (6.283185307179586,)-periodic FourierSeries with Float64 coefficients, (0,) derivative, (-2,) offset
+julia> cosine = FourierSeries([0.5, 0.0, 0.5], period=2pi, offset=-2)
+3-element and (6.283185307179586,)-periodic FourierSeries with Float64 coefficients, (0,) derivative, (-2,) offset
 
-    julia> cosine(pi)
-    -1.0 + 0.0im
+julia> cosine(pi)
+-1.0 + 0.0im
+```
 
 Notice that we can create a series object and interpolate with a function-like interface.
 For more examples, see the documentation.
@@ -23,6 +25,7 @@ For more examples, see the documentation.
 Generic Fourier series can be implemented with the [`AbstractFourierSeries`](@ref) interface
 and the following implementations are provided as building blocks for others:
 - [`FourierSeries`](@ref): (inplace) evaluation at real/complex arguments
+- [`HermitianFourierSeries`](@ref): series coefficients with Hermitian symmetry
 - [`ManyFourierSeries`](@ref): evaluation of multiple series
 - [`DerivativeSeries`](@ref): automatic evaluation of derivatives of series to any order
 
@@ -41,6 +44,8 @@ These routines have the following features
 """
 module FourierSeriesEvaluators
 
+using LinearAlgebra: Hermitian
+
 export AbstractFourierSeries
 include("definitions.jl")
 
@@ -53,4 +58,7 @@ include("fourier_series.jl")
 export FourierWorkspace
 include("workspace.jl")
 
+export hermitian_fourier_evaluate, HermitianFourierSeries
+include("hermitian.jl")
+
 end