A => .gitignore +5 -0
@@ 1,5 @@
+*.csv
+*.png
+!simulationplot.png
+plotsim.gp
+cmocean_balance.pal
A => Manifest.toml +106 -0
@@ 1,106 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[AbstractFFTs]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716"
+uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+version = "0.5.0"
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[FFTW]]
+deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
+git-tree-sha1 = "14536c95939aadcee44014728a459d2fe3ca9acf"
+uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+version = "1.2.2"
+
+[[FFTW_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b"
+uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
+version = "3.3.9+5"
+
+[[IntelOpenMP_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6"
+uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
+version = "2018.0.3+0"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[LibGit2]]
+deps = ["Printf"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[MKL_jll]]
+deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
+git-tree-sha1 = "0ce9a7fa68c70cf83c49d05d2c04d91b47404b08"
+uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
+version = "2020.1.216+0"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Pkg]]
+deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Reexport]]
+deps = ["Pkg"]
+git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "0.2.0"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
A => Project.toml +9 -0
@@ 1,9 @@
+name = "HalfSpectral"
+uuid = "28fda99a-6fd1-4789-996a-91dd14997ef0"
+authors = ["Chris Geoga <cgeoga@protonmail.com>"]
+version = "0.1.0"
+
+[deps]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
A => README.md +47 -0
@@ 1,47 @@
+
+# HalfSpectral.jl
+
+A convenient interface for working with half-spectral covariance functions in
+the space-time domain. This is part of the software companion to *Flexible
+nonstationary spatio-temporal modeling of high-frequency monitoring data* (arxiv
+link coming soon).
+
+See the paper for a discussion of the covariance function itself, and see the
+example file `./example/basicexample.jl` for a quick-start guide. With just a
+few lines of code, you can specify a complicated covariance function
+corresponding to fields like this:
+
+<!---
+
+--->
+<p align="center">
+ <img src="https://git.sr.ht/~cgeoga/HalfSpectral.jl/blob/master/simulationplot.png">
+</p>
+
+For the code used to fit models in the paper, see `./examples/paperscripts/`. I
+do not plan to modify these files if this code evolves, so if you want to run
+those scripts, please use the appropriate tagged release.
+
+For estimation, please see my other software for estimation,
+[GPMaxlik.jl](https://git.sr.ht/~cgeoga/GPMaxlik.jl). The examples for setting
+up and performing optimization in that repository are much simpler than the
+examples in `./examples/paperscripts/`.
+
+**If you use this software (HalfSpectral.jl), please cite the paper and not
+the software package itself.**
+
+# Installation
+
+This package is not registered. Please install by opening Julia, pressing the ]
+key to enter the Pkg REPL, and entering
+````{julia}
+Pkg> add https://git.sr.ht/~cgeoga/HalfSpectral.jl
+````
+
+# Contact
+
+Please reach out to me at <christopher.geoga@rutgers.edu> or <cgeoga@anl.gov> if
+you have any issues.
+
+
A => example/basicexample.jl +57 -0
@@ 1,57 @@
+
+using HalfSpectral, LinearAlgebra
+
+# Marginal spectral density, with the twist of making a few parameters depend on
+# the spatial index.
+function Sf(w,x,p)
+ (scale, rate, smoothness) = p[1:3]
+ rate += log(x)
+ smoothness *= (1+x/10)
+ scale*(rate*sinpi(w)^2 + 1)^(-smoothness-1/2)
+end
+
+# A coherence function that does happen to be stationary. See the paperscripts
+# model for code implementing the Paciorek-Schervish nonstationary coherence.
+function Cw(w,x,y,p)
+ gamxy = p[4]/(1 + (sinpi(w)/0.3)^2)^4
+ exp(-abs(x-y)/gamxy)
+end
+
+# A simple phase function.
+phase(w,x,y,p) = exp(2.0*pi*im*p[5]*sin(sinpi(w))*(x-y))
+
+# The integrand in the form that HalfSpectral requests.
+integrand(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)*phase(w,x,y,p)
+
+# Create the kernel for spatial locations 1:40 and 200 unit time measurements.
+const tg = timegrid(1:200)
+const xpts = 1:40
+const params = (50.0,100.0,0.5,5.0,0.5) # sample parameter values.
+kernel = tftkernel(integrand, tg, xpts, params)
+
+# Build a second time to get a timing example.
+println("Building the kernel takes this long")
+@time kernel = tftkernel(integrand, tg, xpts, params)
+
+# Simulate a realization of the field with this function. The points here should
+# be ordered as a tuple (time, spatial coordinate), and the spatial coordinate
+# does not need to be a real number. You can treat the kernel like any normal
+# function, calling it with either two points x and y like kernel(x,y) or
+# additionally with a parameter collection, and if the parameter collection
+# disagrees with the internal one it will rebuild the dictionary of function
+# values. In general, I suggest including the parameter vector for code safety.
+# But to avoid runtime checks, you can use the unchecked kernel(x,y).
+pts = vec(reverse.(collect(Iterators.product(xpts, tg.Xt))))
+K = Symmetric([kernel(x,y) for x in pts, y in pts])
+sim = reshape(cholesky(K).L*randn(length(pts)), length(xpts), length(tg.Xt))
+
+#= Optionally, save the simulation file:
+using DelimitedFiles
+writedlm("sim.csv", sim, ',')
+=#
+
+#= Optionally, visualize it:
+using Plots
+heatmap(sim, xlabel="time", ylabel="altitude")
+=#
+
A => simulationplot.png +0 -0
A => src/HalfSpectral.jl +119 -0
@@ 1,119 @@
+
+module HalfSpectral
+
+ export timegrid, tftkernel
+
+ using LinearAlgebra, Statistics, FFTW
+
+ abstract type TimeGridResult end
+
+ struct ZeroFunction <: Function
+ end
+ (zrf::ZeroFunction)(x,y,p) = 0.0
+
+ struct IdFunction <: Function
+ end
+ (idf::IdFunction)(x,p) = 1.0
+
+ struct TimeGrid <: TimeGridResult
+ Tv :: Vector{Float64}
+ Xt :: Vector{Int64}
+ dt :: Float64
+ tol :: Float64
+ end
+
+ struct TimeGridFailure <: TimeGridResult
+ Tv ::Vector{Float64}
+ tol ::Float64
+ end
+
+ mutable struct TimeFFTKernel{T,F<:Function,N,G<:Function,H<:Function,M} <: Function
+ V :: TimeGrid
+ Fn :: F
+ D :: Dict{Tuple{Int64, T, T}, Float64}
+ p :: NTuple{N, Float64}
+ X :: AbstractVector{T}
+ addfn :: G
+ mulfn :: H
+ nobuild_ix :: NTuple{M, Int64}
+ padsz::Int64
+ end
+
+ function timegrid(X; tol=0.075, warn_tol=2*tol)::TimeGridResult
+ dt, dtmax = extrema(abs.(diff(X)))
+ d_reldiff = rem(dtmax-dt, dt)/dt
+ if d_reldiff <= warn_tol
+ d_reldiff > tol && (@info "Only the acceptable tol was acheived.")
+ Xd = (X .- X[1])./dt
+ Xd_int = Int64.(round.(Xd))
+ minimum(diff(Xd_int)) > 0 || throw(error("time grid construction failed."))
+ return TimeGrid(X, Xd_int, dt, d_reldiff)
+ end
+ return TimeGridFailure(X, d_reldiff)
+ end
+
+ function tftkernel(integrand::F, tg::TimeGrid, X::AbstractVector{T}, p;
+ addfn::G=ZeroFunction(),
+ mulfn::H=IdFunction(),
+ nob=NTuple{0, Int64}(),
+ padsz::Int64=7) 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)
+ 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])
+ wgd = collect(range(-0.5, 0.5, length=ftn+1)[1:ftn])
+ 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)
+ covxy = real(plan*(fftshift([K.Fn(w,x,y,p) for w in wgd])))[ix.+1]
+ newv = Dict(zip(zip(ix, fill(x, length(ix)), fill(y, length(ix))), covxy))
+ merge!(K.D, newv)
+ end
+ K.p = convert(NTuple{N, Float64}, Tuple(p))
+ return nothing
+ end
+
+ # Avoid the runtime cost of checking the parameters AND the post-function:
+ function (K::TimeFFTKernel{T,F,N,ZeroFunction,IdFunction,0})(x, y) where{T,F,N}
+ opt1 = (abs(x[1]-y[1]), x[2], y[2])
+ haskey(K.D, opt1) && return K.D[opt1]
+ throw(error("No computed covariance for these two points: $x, $y"))
+ end
+
+ # Avoid the runtime cost of checking the parameters, but call the post
+ # computation function:
+ function (K::TimeFFTKernel{T,F,N,G,H,M})(x, y) where{T,F,N,G,H,M}
+ opt1 = (abs(x[1]-y[1]), x[2], y[2])
+ haskey(K.D, opt1) && return K.mulfn(x,y,K.p)*K.D[opt1] + K.addfn(x,y,K.p)
+ throw(error("No computed covariance for these two points: $x, $y"))
+ end
+
+ # Safer checked version:
+ function (K::TimeFFTKernel{T,F,N,M})(x, y, p) where{T,F,N,M}
+ update_p_flag = false
+ for (j, (Kpj, pj)) in enumerate(zip(K.p, p))
+ if Kpj != pj
+ if !in(j, K.nobuild_ix)
+ build!(K, p)
+ update_p_flag = false
+ break
+ else
+ update_p_flag = true
+ end
+ end
+ end
+ if update_p_flag
+ K.p = tuple(p...)
+ end
+ K(x, y)
+ end
+
+end
+
