~cgeoga/HalfSpectral.jl

abea9db714c645f288717bc253fb3d5f511575ac — Chris Geoga 3 years ago 72278d2
Adding an optional 'longmemory' flag for making sure that the frequency
grid of the FFT doesn't contain exactly zero.
1 files changed, 8 insertions(+), 2 deletions(-)

M src/HalfSpectral.jl
M src/HalfSpectral.jl => src/HalfSpectral.jl +8 -2
@@ 37,6 37,7 @@ module HalfSpectral
    mulfn :: H
    nobuild_ix :: NTuple{M, Int64}
    padsz::Int64
    longmemory::Bool
  end

  function timegrid(X; tol=0.075, warn_tol=2*tol)::TimeGridResult


@@ 56,19 57,24 @@ module HalfSpectral
                     addfn::G=ZeroFunction(), 
                     mulfn::H=IdFunction(),
                     nob=NTuple{0, Int64}(),
                     padsz::Int64=7) where{T,F,G,H}
                     padsz::Int64=7,
                     longmemory=false) where{T,F,G,H}
    dict = Dict{Tuple{Int64, T, T}, Float64}()
    lp   = length(p)
    nnb  = length(nob)
    out  = TimeFFTKernel{T,F,lp,G,H,nnb}(tg, integrand, dict, tuple(p...), 
                                         X, addfn, mulfn, nob, padsz)
                                         X, addfn, mulfn, nob, padsz, longmemory)
    build!(out, p)
    return out
  end

  function build!(K::TimeFFTKernel{T,F,N,G,H,M}, p) where{T,F,N,G,H,M}
    ftn   = nextprod([2,3,5,7], K.padsz*K.V.Xt[end])
    (K.longmemory && iseven(ftn)) && (ftn += 1) 
    wgd   = collect(range(-0.5, 0.5, length=ftn+1)[1:ftn])
    if K.longmemory
      @assert !(in(0.0, wgd)) "debug, this should be impossible."
    end
    ix    = sort(unique(map(z->abs(z[2]-z[1]), Iterators.product(K.V.Xt, K.V.Xt))))
    plan  = plan_ifft(wgd)
    for (x,y) in Iterators.product(K.X, K.X)