From 5fbf5adb9374cf2403890dde39d8f3e433216c4d Mon Sep 17 00:00:00 2001 From: Lorenzo Van Munoz <66997677+lxvm@users.noreply.github.com> Date: Fri, 21 Jun 2024 17:25:43 -0400 Subject: [PATCH] Hermitian series (#3) * initial commit * fast evaluate * add tests * documentation --- Project.toml | 11 +- docs/Project.toml | 2 + docs/src/examples.md | 167 +++++++------- docs/src/reference.md | 3 +- ext/FourierSeriesEvaluatorsStaticArraysExt.jl | 7 + src/FourierSeriesEvaluators.jl | 18 +- src/hermitian.jl | 212 ++++++++++++++++++ test/hermitian.jl | 172 ++++++++++++++ test/runtests.jl | 2 + 9 files changed, 501 insertions(+), 93 deletions(-) create mode 100644 ext/FourierSeriesEvaluatorsStaticArraysExt.jl create mode 100644 src/hermitian.jl create mode 100644 test/hermitian.jl diff --git a/Project.toml b/Project.toml index a850d06..8cdd2c8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,16 @@ name = "FourierSeriesEvaluators" uuid = "2a892dea-6eef-4bb5-9d1c-de966c9f6db5" authors = ["lxvm "] -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" diff --git a/docs/Project.toml b/docs/Project.toml index dfa65cd..3e278c7 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/examples.md b/docs/src/examples.md index ae96d6c..fd849f9 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -8,9 +8,9 @@ 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))`. @@ -18,36 +18,26 @@ We can then evaluate `f` at a random point using its function-like interface 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 diff --git a/docs/src/reference.md b/docs/src/reference.md index f8a1b72..15dfd8b 100644 --- a/docs/src/reference.md +++ b/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"] ``` \ No newline at end of file diff --git a/ext/FourierSeriesEvaluatorsStaticArraysExt.jl b/ext/FourierSeriesEvaluatorsStaticArraysExt.jl new file mode 100644 index 0000000..9c70729 --- /dev/null +++ b/ext/FourierSeriesEvaluatorsStaticArraysExt.jl @@ -0,0 +1,7 @@ +module FourierSeriesEvaluatorsStaticArraysExt +import FourierSeriesEvaluators: _hermitianpart +using StaticArrays: SHermitianCompact, StaticMatrix + +_hermitianpart(A::StaticMatrix) = SHermitianCompact(A) + +end diff --git a/src/FourierSeriesEvaluators.jl b/src/FourierSeriesEvaluators.jl index c4bdeec..2dc9445 100644 --- a/src/FourierSeriesEvaluators.jl +++ b/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 diff --git a/src/hermitian.jl b/src/hermitian.jl new file mode 100644 index 0000000..5cc18e2 --- /dev/null +++ b/src/hermitian.jl @@ -0,0 +1,212 @@ +# here we make use of the fact that for coefficients that satisfy A[i] = A[-i]' produce +# hermitian-valued Fourier series since z*A[i] + z'*A[-i] = (z*A[i]) + (z*A[i])' + +_hermitianpart(a::Complex) = real(a) +_hermitianpart(A::AbstractMatrix) = Hermitian(A) + +""" + hermitian_fourier_evaluate(C::AbstractVector, x, [k=1, a=0]) + +Evaluates a 1-D Hermitian Fourier series `C`, which is a compact representation +of a series with coefficients `C[i] = C[-i]'`. The other arguments are the same +as [`fourier_evaluate`](@ref). The precise formula of what is evaluated is given +by the following, where the initial element ``C_{0}`` is equal to `C[begin]`: +```math +r = 0^a C_{0} + \\sum_{n=1}^{N-1} C_{n} (ik n)^{a} \\exp(ik n x)) + \\text{h.c.} +``` +""" +@fastmath function hermitian_fourier_evaluate(C::AbstractVector, x, k=inv(oneunit(x)), a=Val(0)) + u = k*x # unitless + @assert u isa Real "Hermitian series can only be evaluated at real arguments" + s = size(C, 1) + i = firstindex(C) # Find the index offset + e = complex(one(u)) # initial offset phase + + # apply global Fourier multiplier (im*k)^a (k^a could change type) + z = _iszero(a) ? e : + _isone(a) ? e*im*k : + _istwo(a) ? e*(im*k)^2 : e*_pow(im*k, a) + + # initial Fourier coefficient z*c0^a (won't change type) + b = _iszero(a) ? z : 0z + @inbounds r = _hermitianpart(b * C[i]) # unroll first loop iteration + s == 1 && return r # fast return for singleton dimensions + + dz = cis(1u) # forward-backward-winding phase steps + for n in 1:s-1 # symmetric winding loop + z *= dz # forward--winding phase + i += 1 # forward--winding index + # winding coefficients + b = _iszero(a) ? z : + _isone(a) ? z*n : + _istwo(a) ? z*n^2 : z*_pow(n, a) + + @inbounds B = b * C[i] + r += _hermitianpart(B) + _hermitianpart(B') + end + return r +end +#= +@fastmath function hermitian_fourier_contract!(r::AbstractArray, C::AbstractArray, x, k=inv(oneunit(x)), a=Val(0)) + ndims(C) != ndims(r)+1 && throw(ArgumentError("array dimensions incompatible")) + ax = axes(r) + axC = axes(C) + ax == axC[1:end-1] || throw(BoundsError("array axes incompatible (check size or indexing)")) + + preI = CartesianIndices(ax) + + u = k*x # unitless + @assert u isa Real "Hermitian series can only be evaluated at real arguments" + s = size(C, ndims(C)) + i = firstindex(C, ndims(C)) # Find the index offset + e = complex(one(u)) # initial offset phase + + # apply global Fourier multiplier (im*k)^a (k^a could change type) + z = _iszero(a) ? e : + _isone(a) ? e*im*k : + _istwo(a) ? e*(im*k)^2 : e*_pow(im*k, a) + + # initial Fourier coefficient z*c0^a (won't change type) + b = _iszero(a) ? z : 0z + for prefix in preI + @inbounds r[prefix] = _hermitianpart(b * C[prefix, i]) + end + s == 1 && return r # fast return for singleton dimensions + + dz = cis(1u) # forward-backward-winding phase steps + for n in 1:s-1 # symmetric winding loop + z *= dz # forward--winding phase + i += 1 # forward--winding index + # winding coefficients + b = _iszero(a) ? z : + _isone(a) ? z*n : + _istwo(a) ? z*n^2 : z*_pow(n, a) + + for prefix in preI + @inbounds B = b * C[prefix, i] + @inbounds r[prefix] += _hermitianpart(B) + _hermitianpart(B') + end + end + return r +end + +""" + hermitian_fourier_allocate(C, x, k, a) + +Allocate an array of the correct type for contracting the Hermitian Fourier +series along the last axis +""" +function hermitian_fourier_allocate(C::AbstractArray{T,N}, x, k, a) where {T,N} + ax = axes(C) + prefix = CartesianIndices(ax[1:N-1]) + v = view(C, prefix, firstindex(C, N)) + y = zero(T)*cis(k*x) + r = _iszero(a) ? y : + _isone(a) ? y*im*k : + _istwo(a) ? y*(im*k)^2 : y*_pow(im*k, a) + return similar(v, typeof(_hermitianpart(r))) +end +=# +struct HermitianFourierSeries{N,iip,C<:AbstractArray,A<:NTuple{N,Any},T,F} <: AbstractFourierSeries{N,T,iip} + c::C + a::A + t::NTuple{N,T} + f::NTuple{N,F} +end + +function _check_coeffs_hermsym(c::AbstractArray, o::NTuple) + ax = axes(c) + _ax = ax[ndims(c)-length(o)+1:ndims(c)] + ax_ = ax[1:ndims(c)-length(o)] + for i in CartesianIndices(ax_) + _c = @view c[i,CartesianIndices(_ax)] + # TODO single, not double, check + all(c -> c[begin] == c[end]', zip(_c, Iterators.reverse(_c))) || return false + end + return true +end +function _check_coeffs_centered(c::AbstractArray, o::NTuple) + @assert ndims(c) >= length(o) + all(zip(axes(c)[end-length(o)+1:end],o)) do (ax, o) + len = length(ax) + isodd(len) && iszero(ax[div(len,2)+begin]+o) + end +end + +#= +Plan for inplace evaluation +If `f` is an inplace evaluator, the Hermitian evaluator may evaluate a different result: +- when `ndims(f)+1 = ndims(f.c)`, the Hermitian evaluator will do the same as + `f`, which is to apply the same summation to each element, like a batch evaluator +- when `ndims(f)+2 = ndims(f.c)`, and `eltype(f.c) isa Number` the Hermitian + evaluator will treat the inner dimensions as the indices of an operator. The output +- when `ndims(f)+3 = ndims(f.c)`, and `eltype(f.c) isa Number` the Hermitian + evaluator will treat the two inner dimensions as the indices of an operator, + and the third innermost dimension as a batching dimension +- Otherwise, the behavior of the Hermitian evaluators is not defined (it will error) +=# +""" + HermitianFourierSeries(f::FourierSeries) + +Construct an efficient evaluator for an operator-valued [`FourierSeries`](@ref) +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 +""" +function HermitianFourierSeries(f::FourierSeries) + isinplace(f) && throw(ArgumentError("inplace series are not supported")) + _check_coeffs_hermsym(f.c, f.o) || throw(ArgumentError("Fourier coefficients do not have Hermitian symmetry")) + _check_coeffs_centered(f.c, f.o) + ax = axes(f.c) + _ax = ax[ndims(f.c)-length(f.o)+2:ndims(f.c)] + ax_ = ax[1:ndims(f.c)-length(f.o)] + _ax_ = ax[ndims(f.c)-length(f.o)+1] + c = f.c[CartesianIndices(ax_),_ax_[begin+div(length(_ax_),2):end],CartesianIndices(_ax)] # copy + return HermitianFourierSeries{ndims(f),isinplace(f),typeof(c),typeof(f.a),eltype(f.t),eltype(f.f)}(c, f.a, f.t, f.f) +end + +function evaluate!(c, s::HermitianFourierSeries{1}, x) + f = freq2rad(s.f[1]) + if isinplace(s) + error("not implemented") + # return hermitian_fourier_contract!(c, s.c, x, f, s.a[1]) + else + return hermitian_fourier_evaluate(s.c, x, f, s.a[1]) + end +end + +function allocate(s::HermitianFourierSeries, x, ::Val{d}) where {d} + if isinplace(s) || ndims(s) > 1 + M = ndims(s.c)-ndims(s) + return fourier_allocate(s.c, x, freq2rad(s.f[d]), s.a[d], Val(M+d)) + elseif isinplace(s) && ndims(s) == 1 + error("not implemented") + # return hermitian_fourier_allocate(s.c, x, freq2rad(s.f[d]), s.a[d]) + else + return nothing + end +end + +function contract!(c, s::HermitianFourierSeries, x, dim::Val{d}) where {d} + if d == 1 + error("not implemented") + # hermitian_fourier_contract!(c, s.c, x, freq2rad(s.f[d]), s.a[d]) + else + M = ndims(s.c)-ndims(s) + fourier_contract!(c, s.c, x, freq2rad(s.f[d]), s.a[d], -div(size(s.c, M+d),2)-1, Val(M+d)) + end + t = deleteat_(s.t, dim) + f = deleteat_(s.f, dim) + a = deleteat_(s.a, dim) + return HermitianFourierSeries{ndims(s)-1,isinplace(s),typeof(c),typeof(a),eltype(t),eltype(f)}(c, a, t, f) +end + +period(s::HermitianFourierSeries) = s.t +frequency(s::HermitianFourierSeries) = s.f + +function nextderivative(s::HermitianFourierSeries, dim) + a = raise_multiplier(s.a, dim) + return HermitianFourierSeries{ndims(s),isinplace(s),typeof(s.c),typeof(a),eltype(s.t),eltype(s.f)}(s.c, a, s.t, s.f) +end \ No newline at end of file diff --git a/test/hermitian.jl b/test/hermitian.jl new file mode 100644 index 0000000..21dd121 --- /dev/null +++ b/test/hermitian.jl @@ -0,0 +1,172 @@ +using Test +using LinearAlgebra + +using StaticArrays + +using FourierSeriesEvaluators +using FourierSeriesEvaluators: _check_coeffs_hermsym, _check_coeffs_centered +using FourierSeriesEvaluators: workspace_allocate, workspace_evaluate + +@testset "hermitian checks" begin + n = 3 + c = rand(SMatrix{3,3,ComplexF64,9}, 2n+1,2n+1) + @test !_check_coeffs_hermsym(c, (nothing, nothing)) + c = c + reverse(adjoint.(c)) + @test _check_coeffs_hermsym(c, (nothing, nothing)) + @test !_check_coeffs_centered(c, (-n, -n)) + @test _check_coeffs_centered(c, (-n-1, -n-1)) + #= + ninner=4 + cc = Array{ComplexF64}(undef, ninner, 2n+1,2n+1) + @test !_check_coeffs_hermsym(cc, (nothing, nothing)) + for i in axes(cc,1) + h = rand(ComplexF64, 2n+1, 2n+1) + cc[i,:,:] .= h + reverse(adjoint.(h)) + end + @test _check_coeffs_hermsym(cc, (nothing, nothing)) + @test !_check_coeffs_centered(cc, (-n, -n)) + @test _check_coeffs_centered(cc, (-n-1, -n-1)) + =# +end + +@testset "evaluate" for d in 1:2 + n = 5 + ntest = 10 + T = Float64 + idx = Iterators.product(repeat([1:(2n+1)], d)...) + for a in [ + [rand(Complex{T}) for _ in idx], + # [rand(Complex{T}, 2, 2) for _ in idx], + # [rand(Complex{T}, 3, 3) for _ in idx], + [rand(SMatrix{2,2,Complex{T},4}) for _ in idx], + [rand(SMatrix{3,3,Complex{T},9}) for _ in idx], + ] + c = a + reverse(adjoint.(a)) + f = FourierSeries(c; period=1.0, offset=-n-1) + if first(c) isa AbstractMatrix && d == 1 + f2 = FourierSeries(UpperTriangular.(c); period=1.0, offset=-n-1) + end + hf = HermitianFourierSeries(f) + for _ in 1:ntest + x = rand(T, d) + hfx = @inferred hf(x) + @test ishermitian(hfx) + @test f(x) ≈ hfx + if first(c) isa AbstractMatrix && d == 1 + @test hfx ≈ Hermitian(f2(x)) + end + end + end +end + +@testset "inplace" for d in 1:2 + n = 5 + ntest = 10 + ninner = 4 + T = Float64 + idx = Iterators.product(ninner, repeat([1:(2n+1)], d)...) + for a in [ + [rand(Complex{T}) for _ in idx], + # [rand(Complex{T}, 3, 3) for _ in idx], + [rand(SMatrix{2,2,Complex{T},4}) for _ in idx], + [rand(SMatrix{3,3,Complex{T},9}) for _ in idx], + ] + c = Array{eltype(a)}(undef, ninner, 2n+1,2n+1) + for i in axes(c,1) + c[i,:,:] .= a[i,:,:] + reverse(adjoint.(a[i,:,:])) + end + f = FourierSeries(c, d; period=1.0, offset=-n-1) + @test_throws ArgumentError hf = HermitianFourierSeries(f) + #= + for _ in 1:ntest + x = rand(T, d) + hfx = @inferred hf(x) + @test all(ishermitian, hfx) + @test f(x) ≈ hfx + end + =# + end +end + +@testset "derivative" begin + order = 2 + d = 2 + n = 5 + ntest = 3 + idx = Iterators.product(repeat([1:(2n+1)], d)...) + T = Float64 + c = [rand(SMatrix{2,2,Complex{T},4}) for _ in idx] + c = c + reverse(adjoint.(c)) + f = FourierSeries(c; period=1.0, offset=-n-1) + hf = HermitianFourierSeries(f) + df = DerivativeSeries{order}(f) + dhf = DerivativeSeries{order}(hf) + for i in 1:ntest + x = rand(d) + dfx = @inferred df(x) + dhfx = @inferred dhf(x) + for (n,(dfxi, dhfxi)) in enumerate(zip(dfx, dhfx)) + if n == (0+1) + @test ishermitian(dhfxi) + @test dfxi ≈ dhfxi + elseif n == (1+1) + for i in 1:d + @test ishermitian(dhfxi[i]) + @test dfxi[i] ≈ dhfxi[i] + end + elseif n == (2+1) + for j in 1:d + for i in 1:d-j+1 + @test ishermitian(dhfxi[j][i]) + @test dfxi[j][i] ≈ dhfxi[j][i] + end + end + elseif n == (3+1) + for k in 1:d + for j in 1:d-j+1 + for i in 1:d-j-k+2 + @test ishermitian(dhfxi[k][j][i]) + @test dfxi[k][j][i] ≈ dhfxi[k][j][i] + end + end + end + else + error("test not implemented for order") + end + end + end +end + +@testset "workspace" begin + d = 3 + n = 3 + T = SMatrix{2,2,ComplexF64,4} + C = rand(T, ntuple(_ -> n, d)...) + #= + # inplace + for nvar in 1:d-1 + periods = rand(nvar) + CC = Array{eltype(C)}(undef, size(C)...) + id = CartesianIndices(axes(C)[d-nvar+1:end]) + for i in CartesianIndices(axes(C)[1:d-nvar]) + CC[i,id] .= C[i,id] + reverse(adjoint.(C[i,id])) + end + fip = FourierSeries(CC, nvar, period=periods) + hfip = HermitianFourierSeries(fip) + x = tuple(rand(nvar)...) + for sip in (hfip, ManyFourierSeries(hfip,hfip,period=periods), JacobianSeries(hfip)) + wsip = workspace_allocate(sip, x) + @test workspace_evaluate(wsip, x) == sip(x) + end + end + =# + C = C + reverse(adjoint.(C)) + periods = rand(d) + f = FourierSeries(C, period=periods) + hf = HermitianFourierSeries(f) + x = tuple(rand(d)...) + for s in (hf, ManyFourierSeries(hf,hf,period=periods), JacobianSeries(hf)) + ws = workspace_allocate(s, x) + @test workspace_evaluate(ws, x) == s(x) + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 2047511..acfa281 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -206,3 +206,5 @@ include("fourier_reference.jl") end end end + +@testset "hermitian" begin include("hermitian.jl") end