text:
using IterativeSolvers
using LinearMaps
using LinearAlgebra
using Test
using Random

# Already defined in another file
#=
include("laplace_matrix.jl")

struct JacobiPrec{TD}
    diagonal::TD
end

Base.ldiv!(y, P::JacobiPrec, x) = y .= x ./ P.diagonal
=#

function max_err(R)
    r = zeros(real(eltype(R)), size(R, 2))
    for j in 1:length(r)
        for i in 1:size(R, 1)
            r[j] += conj(R[i,j])*R[i,j]
        end
        r[j] = sqrt(r[j])
    end
    maximum(r)
end

@testset "Locally Optimal Block Preconditioned Conjugate Gradient" begin
    srand(1234323)
    @testset "Single eigenvalue" begin
        n = 10
        @testset "Small full system" begin
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        b = rand(T, n, 1)
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, largest, b; tol=tol, maxiter=Inf, log=false)
                        λ, X = r.λ, r.X
                        @test norm(A*X - X*λ) ≤ tol

                        # If you start from the exact solution, you should converge immediately
                        r = lobpcg(A, largest, X; tol=10tol, log=true)
                        @test length(r.trace) == 1
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        b = rand(T, n, 1)
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, B, largest, b; tol=tol, maxiter=Inf, log=true)
                        λ, X = r.λ, r.X
                        @test max_err(A*X - B*X*λ) ≤ tol

                        # If you start from the exact solution, you should converge immediately
                        r = lobpcg(A, B, largest, X; tol=10tol, log=true)
                        @test length(r.trace) == 1
                    end
                end
            end
        end
        @testset "Sparse Laplacian" begin
            A = laplace_matrix(Float64, 20, 2)
            rhs = randn(size(A, 2), 1)
            rmul!(rhs, inv(norm(rhs)))
            tol = IterativeSolvers.default_tolerance(Float64)
            @testset "Matrix" begin
                @testset "largest = $largest" for largest in (true, false)
                    r = lobpcg(A, largest, rhs; tol=tol, maxiter=Inf)
                    λ, X = r.λ, r.X
                    @test norm(A*X - X*λ) ≤ tol
                end
            end
        end
        @testset "Zero initial solution" begin
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        b = zeros(T, n, 1)
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, largest, b; tol=tol, maxiter=Inf, log=false)
                        λ, X = r.λ, r.X
                        @test norm(A*X - X*λ) ≤ tol
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        b = zeros(T, n, 1)
                        tol = IterativeSolvers.default_tolerance(T)

                        r = lobpcg(A, B, largest, b; tol=tol, maxiter=Inf, log=true)
                        λ, X = r.λ, r.X
                        @test max_err(A*X - B*X*λ) ≤ tol
                    end
                end
            end
        end
        @testset "No initial solution" begin
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        tol = IterativeSolvers.default_tolerance(T)

                        r = lobpcg(A, largest, 1; tol=tol, maxiter=Inf, log=false)
                        λ, X = r.λ, r.X
                        @test norm(A*X - X*λ) ≤ tol
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        tol = IterativeSolvers.default_tolerance(T)

                        r = lobpcg(A, B, largest, 1; tol=tol, maxiter=Inf, log=true)
                        λ, X = r.λ, r.X
                        @test max_err(A*X - B*X*λ) ≤ tol
                    end
                end
            end
        end
        @testset "Inplace" begin
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        tol = IterativeSolvers.default_tolerance(T)
                        b = rand(T, n, 1)
                        itr = LOBPCGIterator(A, largest, b)

                        r = lobpcg!(itr; tol=tol, maxiter=Inf, log=false)
                        λ, X = r.λ, r.X
                        @test norm(A*X - X*λ) ≤ tol
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        b = rand(T, n, 1)
                        tol = IterativeSolvers.default_tolerance(T)
                        itr = LOBPCGIterator(A, B, largest, b)

                        r = lobpcg!(itr; tol=tol, maxiter=Inf, log=true)
                        λ, X = r.λ, r.X
                        @test max_err(A*X - B*X*λ) ≤ tol
                    end
                end
            end
        end
        @testset "Jacobi preconditioner" begin
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        tol = IterativeSolvers.default_tolerance(T)
                        P = JacobiPrec(diag(A))
                        r = lobpcg(A, largest, 1; P=P, tol=tol, maxiter=Inf, log=false)
                        λ, X = r.λ, r.X
                        @test norm(A*X - X*λ) ≤ tol
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        P = JacobiPrec(diag(A))
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        tol = IterativeSolvers.default_tolerance(T)

                        r = lobpcg(A, B, largest, 1; P=P, tol=tol, maxiter=Inf, log=true)
                        λ, X = r.λ, r.X
                        @test max_err(A*X - B*X*λ) ≤ tol
                    end
                end
            end
        end
        @testset "Constraint" begin
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, largest, 1; tol=tol, maxiter=Inf, log=false)
                        λ1, X1 = r.λ, r.X
                        r = lobpcg(A, largest, 1; C=copy(r.X), tol=tol, maxiter=Inf, log=false)
                        λ2, X2 = r.λ, r.X
                        @test norm(A*X2 - X2*λ2) ≤ tol
                        @test isapprox(real((adjoint(X1)*X2)[1,1]), 0, atol=2*n*tol)
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, B, largest, 1; tol=tol, maxiter=Inf, log=false)
                        λ1, X1 = r.λ, r.X
                        r = lobpcg(A, B, largest, 1; C=copy(r.X), tol=tol, maxiter=Inf, log=false)
                        λ2, X2 = r.λ, r.X
                        @test norm(A*X2 - B*X2*λ2) ≤ tol
                        @test isapprox(real((adjoint(X1)*(B*X2))[1,1]), 0, atol=2*n*tol)
                    end
                end
            end
        end
    end
    @testset "Two eigenvalues" begin
        @testset "Small full system" begin
            n = 10
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        b = rand(T, n, 2)
                        tol = IterativeSolvers.default_tolerance(T)

                        r  = lobpcg(A, largest, b; tol=tol, maxiter=Inf, log=false)
                        λ, X = r.λ, r.X
                        @test max_err(A*X - X*Matrix(Diagonal(λ))) ≤ tol

                        # If you start from the exact solution, you should converge immediately
                        r = lobpcg(A, largest, X; tol=10tol, log=true)
                        @test length(r.trace) == 1
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        b = rand(T, n, 2)
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, B, largest, b; tol=tol, maxiter=Inf, log=true)
                        λ, X = r.λ, r.X
                        @test max_err(A*X - B*X*Matrix(Diagonal(λ))) ≤ tol

                        # If you start from the exact solution, you should converge immediately
                        r = lobpcg(A, B, largest, X; tol=10tol, log=true)
                        @test length(r.trace) == 1
                    end
                end
            end
        end
    end
    @testset "nev = 3, block size = $block_size" for block_size in (1, 2)
        n = 10
        @testset "Simple eigenvalue problem" begin
            @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                @testset "largest = $largest" for largest in (true, false)
                    A = rand(T, n, n)
                    A = A' + A + 20I
                    tol = IterativeSolvers.default_tolerance(T)
                    X0 = rand(T, n, block_size)
                    r = lobpcg(A, largest, X0, 3, tol=tol, maxiter=Inf, log=true)
                    λ, X = r.λ, r.X
                    @test max_err(A*X - X*Matrix(Diagonal(λ))) ≤ tol
                    @test all(isapprox.(adjoint(X)*X, Matrix{T}(I, 3, 3), atol=2*n*tol))
                end
            end
        end
        @testset "Generalized eigenvalue problem" begin
            @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                @testset "largest = $largest" for largest in (true, false)
                    A = rand(T, n, n)
                    A = A' + A + 20I
                    B = rand(T, n, n)
                    B = B' + B + 20I
                    tol = IterativeSolvers.default_tolerance(T)
                    X0 = rand(T, n, block_size)
                    r = lobpcg(A, B, largest, X0, 3, tol=tol, maxiter=Inf, log=true)
                    λ, X = r.λ, r.X
                    @test max_err(A*X - B*X*Matrix(Diagonal(λ))) ≤ tol
                    @test all(isapprox.(adjoint(X)*(B*X), Matrix{T}(I, 3, 3), atol=2*n*tol))
                end
            end
        end
        @testset "Constraint" begin
            @testset "Simple eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, largest, 1; tol=tol, maxiter=Inf, log=false)
                        λ1, X1 = r.λ, r.X

                        X0 = rand(T, n, block_size)
                        r = lobpcg(A, largest, X0, 3, C=copy(r.X), tol=tol, maxiter=Inf, log=true)
                        λ2, X2 = r.λ, r.X
                        @test max_err(A*X2 - X2*Matrix(Diagonal(λ2))) ≤ tol
                        @test all(isapprox.(adjoint(X2)*X2, Matrix{T}(I, 3, 3), atol=2*n*tol))
                        @test all(isapprox.(real(adjoint(X1)*X2), 0, atol=2*n*tol))
                    end
                end
            end
            @testset "Generalized eigenvalue problem" begin
                @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
                    @testset "largest = $largest" for largest in (true, false)
                        A = rand(T, n, n)
                        A = A' + A + 20I
                        B = rand(T, n, n)
                        B = B' + B + 20I
                        tol = IterativeSolvers.default_tolerance(T)
                        r = lobpcg(A, B, largest, 1; tol=tol, maxiter=Inf, log=false)
                        λ1, X1 = r.λ, r.X

                        X0 = rand(T, n, block_size)
                        r = lobpcg(A, B, largest, X0, 2, C=copy(r.X), tol=tol, maxiter=Inf, log=true)
                        λ2, X2 = r.λ, r.X
                        @test max_err(A*X2 - B*X2*Matrix(Diagonal(λ2))) ≤ tol
                        @test all(isapprox.(adjoint(X2)*(B*X2), Matrix{T}(I, 2, 2), atol=2*n*tol))
                        @test all(isapprox.(real(adjoint(X1)*(B*X2)), 0, atol=2*n*tol))
                    end
                end
            end
        end
    end
end


text:
# ============================================
#  Barycenter example
# ============================================

# --------------------------------------------
#  Load Packages, set paths and load the data.
# --------------------------------------------

using DataFrames
using JuMP
using Clp
using PyCall
using PyPlot

# ---- Set directory paths to code and the paper
figurepath = joinpath(pwd(),"figures")
codepath   = joinpath(pwd(),"scripts")

# ---- Set directory paths to shapefile for mappitn california
# citation: from http://www.census.gov/cgi-bin/geo/shapefiles2010/main
shapename = "tl_2010_06_state10"
shapefile = joinpath(pwd(),"mapshapefiles/$shapename/$shapename")

# -----------  Load the data
# citation: from https://en.wikipedia.org/wiki/Climate_of_California
citydf = DataFrame(
	city = [
		:Bakersfield,
		:Eureka,
		:Fresno,
		:LosAngeles,
		:Sacramento,
		:SanBernardino,
		:SanFrancisco,
		:SanJose,
		:SouthLakeTahoe ],
	lat_y      = [35.3667,   40.8019,   36.7500,   34.0500,   38.5556,	 34.1000,   37.7833,	37.3382,	38.9400],
	lon_x      = [-119.0167, -124.1636, -119.7667, -118.2500, -121.4689, -117.3000, -122.4167,	-121.8863,	-119.9769],
	population = [363_630,   26_913,    509_924,   3_884_000, 479_686,	 2_088_000, 837_442,	998_537,	21_387],
	jan = [57, 56, 55, 71, 54, 69, 58, 61, 42],
	feb = [64, 57, 62, 72, 61, 71, 61, 64, 43],
	mar = [70, 57, 68, 75, 66, 75, 63, 67, 48],
	apr = [76, 58, 75, 78, 72, 80, 64, 73, 53],
	may = [84, 61, 84, 80, 80, 85, 66, 77, 63],
	jun = [92, 63, 92, 83, 87, 92, 68, 82, 72],
	jul = [98, 64, 99, 90, 92, 10, 68, 84, 80],
	aug = [96, 65, 97, 92, 92, 103,69, 84, 80],
	sep = [91, 64, 91, 91, 88, 100, 71, 81, 73],
	oct = [80, 62, 79, 84, 78, 86, 70, 76, 62],
	nov = [66, 58, 65, 78, 64, 75, 64, 65, 50],
	dec = [57, 55, 55, 71, 54, 68, 58, 61, 42],
)

# ---------------------------------------------------
# Transform temp data and package into a custom type
# ---------------------------------------------------
# This is basically a named array type.
# It is designed to contain the images in the columns (indexed by month)
immutable city_x_month_image <: AbstractMatrix{Float64}
	imagematrix::Array{Float64,2}
	city_nms ::Array{Symbol,1}
	month_nms::Array{Symbol,1}
end

# here is a constructor
# This allows one to
function city_x_month_image(df::DataFrame, months::Symbol...)
	nrw = length(df[:city])
	ncl = length(months)
	transform_normalize(tempu, pop)  = (z=pop.*((tempu.-72).^2); return z./sum(z))
	imageslocal = zeros(Float64, (nrw, ncl))
	for cl = 1:ncl
		imageslocal[:,cl] = transform_normalize(citydf[months[cl]], citydf[:population])
	end
	return city_x_month_image(imageslocal, df[:city], [months...])
end
city_x_month_image(df::DataFrame) = city_x_month_image(df::DataFrame, :jan,:feb,:mar,:apr,:may,:jun,:jul,:aug,:sep,:oct,:nov,:dec)

Base.size(im::city_x_month_image)   = size(im.imagematrix)
Base.length(im::city_x_month_image) = length(im.imagematrix)
Base.getindex(im::city_x_month_image, i::Symbol, j::Symbol) = im.imagematrix[find(im.city_nms .== i)[1], find(im.month_nms .== j)[1]]
Base.getindex(im::city_x_month_image, i::Symbol, j) = im.imagematrix[find(im.city_nms .== i)[1], j]
Base.getindex(im::city_x_month_image, i, j::Symbol) = im.imagematrix[i, find(im.month_nms .== j)[1]]
Base.getindex(im::city_x_month_image, i, j) = im.imagematrix[i, j]
months(im::city_x_month_image) = im.month_nms
cities(im::city_x_month_image) = im.city_nms

# make an instance of this type
images = city_x_month_image(citydf, :dec, :jan, :feb, :mar, :jun, :jul, :aug, :sep)
# now we can index as follows: images[:Sacramento, :jan]

# ---------------------------------------------------
#   Get barycenter support
# ---------------------------------------------------

# The following object iterates over the number of ways to choose 8 locations out of 20.
function get_support(citydf, images)
	ncities      = length(cities(images))
	nmonths      = length(months(images))
	ndividers    = ncities - 1          # 8
	nones        = nmonths + ndividers  # 20
	itercomb     = combinations(1:nones, ndividers)

	bary_support = zeros(Float64, (2, length(itercomb))) # columns give the support coordinates
	for (index, value) in enumerate(itercomb)
		tmp        = vcat(value[1], diff(sort(value))) .- 1 # this gives the counts for the first 8  partitions
		citycounts = vcat(tmp, nmonths - sum(tmp))
		for i = 1:ncities
			bary_support[:, index] .+= citycounts[i].*[citydf[i, :lon_x], citydf[i, :lat_y]]
		end
		bary_support[:, index] .*= (1.0 / nmonths)
	end
	return bary_support
end
@time bary_support = get_support(citydf, images)

# ------------------------------------------------------------
# Define and solve the Wasserstein Barycenter linear program
# -----------------------------------------------------------

# ------- initialize distance matrix and variable size
N = size(images, 2)        # number of images
n = size(images, 1)        # each image size
b = size(bary_support, 2)  # barycenter size
zk = (citydf[:lon_x] .- bary_support[1,:]).^2  + (citydf[:lat_y] .- bary_support[2,:]).^2

# ------- define the solver
m = Model(solver = ClpSolver(SolveType=5))

# ------ set variables and constraints
@defVar(m, p[1:b] >= 0 )
@defVar(m, mk[1:n, 1:b, 1:N] >= 0 )
for i = 1:n, k = 1:N
	@addConstraint(m, sum{mk[i, j, k], j = 1:b} == images[i, k])
end
for j = 1:b, k = 1:N
	@addConstraint(m, sum{mk[i, j, k], i = 1:n} == p[j])
end

# ------ set objective
@setObjective(m, Min, sum{zk[i, j] * mk[i, j, k], i=1:n, j=1:b, k=1:N})

# ------- solve it
@time status = solve(m)

# ----- get the solution and the support
pval  = getValue(p[:])
mkval = Array(Float64, n, b, N)
for kntr = 1:N
	mkval[:,:,kntr] = getValue(mk[:,:,kntr])
end
psupport = find(pval)

# ------------------------------------------------------------
# Generate figure1* (the individual month demands)
# ------------------------------------------------------------

@pyimport mpl_toolkits.basemap as basemap
Basemap =  basemap.Basemap
lat_0, lon_0 =37.0, -119.3
map = Basemap(
	projection="merc",
	lat_0=lat_0,
	lon_0=lon_0,
    resolution = "h",
    area_thresh=1000,
    llcrnrlon=lon_0 - 5.2 ,
    llcrnrlat=lat_0 - 5.2,
    urcrnrlon=lon_0 + 5.2 ,
    urcrnrlat=lat_0 + 5.2
    )

# ---- plot individual months
function monthplots(setmonth::Symbol)
	map[:readshapefile](shapefile, shapename, linewidth=1.5)
	map[:scatter](citydf[:lon_x], citydf[:lat_y], latlon = true, s=5000*images[:,find(months(images) .== setmonth)], c=[0,0,0], alpha=0.5)
	map[:scatter](citydf[:lon_x], citydf[:lat_y], latlon = true,
		s=50,
		linewidths=2,
		marker="x",
		c=[1,0,0],
		alpha=1)
	axis("off")
end

monthplots(:feb)
savefig(joinpath(figurepath, "figure1a.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

monthplots(:mar)
savefig(joinpath(figurepath, "figure1b.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

monthplots(:jun)
savefig(joinpath(figurepath, "figure1c.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

monthplots(:jul)
savefig(joinpath(figurepath, "figure1d.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

# -------------------------------------------------------------------
# Generate figure2* (the barycenter and the possible support points)
# -------------------------------------------------------------------

# ---- barycenter  plot
map[:readshapefile](shapefile, shapename, linewidth=1.5)
map[:scatter](bary_support[1,psupport][:], bary_support[2,psupport][:],  latlon = true,
	s=5000*pval[psupport],
	c=[0,0,0],
	alpha=0.5)
map[:scatter](citydf[:lon_x], citydf[:lat_y], latlon = true,
	s=50,
	linewidths=2,
	marker="x",
	c=[1,0,0],
	alpha=1)
axis("off")
savefig(joinpath(figurepath, "figure1e.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

# ---- barycenter support,
map[:readshapefile](shapefile, shapename, linewidth=1.5)
map[:scatter](bary_support[1,:][:], bary_support[2,:][:], latlon = true, s = 10,  marker=".", c=[1,1,1], alpha=0.2)
map[:scatter](citydf[:lon_x], citydf[:lat_y], latlon = true,
	s=50,
	linewidths=2,
	marker="x",
	c=[1,0,0],
	alpha=1)
axis("off")
savefig(joinpath(figurepath, "figure1f.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

# -------------------------------------------------------------------
# Generate figure3* (the transportation maps)
# -------------------------------------------------------------------

# barycntr_trnsprt_map:
#	- Each row corresponds to a barycenter support point
#   - Columns will correspond to month and contain the name of the servicing city.
barycntr_trnsprt_map = DataFrame(bary_supp_indx = psupport)

for iim = 1:length(months(images))   # iim indexes the months
	findcity = Array(Int64, length(psupport)) # initialize the container that  this should contain entries 1,2,...,9 which correspond to the rows of images.
	for kr = 1:length(psupport)
		_, findcity[kr] = findmax(mkval[:, psupport[kr], iim])
	end
	barycntr_trnsprt_map[months(images)[iim]] = cities(images)[findcity]
end
@assert barycntr_trnsprt_map[:bary_supp_indx] == psupport

# cit_srvic_fun:
#   - A function cit_srvic_map(:month, :city) => the barycenter support points which serve that city for that month
cit_srvic_fun(mnth, cit) =  (barycntr_trnsprt_map[cit .== barycntr_trnsprt_map[mnth], :bary_supp_indx])
for mnth in months(images)
	@assert length(psupport) == sum([length(cit_srvic_fun(mnth,city)) for city in cities(images)])
end

# ---- transportation into :sink_city on month :month
function makequiver(monthindex::Int)
	println(months(images)[monthindex])
	map[:readshapefile](shapefile, shapename, linewidth=1.5)
	for sink_city in [:SanFrancisco, :Sacramento, :LosAngeles, :SanBernardino]
		month      = months(images)[monthindex] # pick a month [:mar, :sep]
		b_sply_inx = cit_srvic_fun(month, sink_city)
		suppl_x, suppl_y  = map[:__call__](bary_support[1, b_sply_inx][:], bary_support[2, b_sply_inx][:])
		sinkl_x, sinkl_y  = map[:__call__](citydf[find(citydf[:city].==sink_city)[1], :lon_x], citydf[find(citydf[:city].==sink_city)[1],:lat_y])
		map[:scatter](suppl_x, suppl_y,
			s=5000*pval[b_sply_inx],
			c=[0,0,0],
			alpha=0.15)
		function smallt(u,v)
			nrm = √(u.^2 + v.^2)
			return 30_000 .* u ./ nrm, 30_000 .* v ./ nrm
		end
		sub1, sub2 = smallt(sinkl_x .- suppl_x, sinkl_y .- suppl_y)
		quiver(suppl_x, suppl_y, sinkl_x .- suppl_x .- sub1, sinkl_y .- suppl_y .- sub2,
			angles="xy",
			scale_units="xy",
			scale = 1,
			width = 0.004,
			headwidth = 4.9,
			headlength =  3.0,
			headaxislength =  3.5,
			minshaft =  2.5, #1.5
			color = [0,0,0],
			alpha=0.9)
		map[:scatter](sinkl_x, sinkl_y,
			s=50,
			linewidths=2,
			marker="x",
			c=[1,0,0],
			alpha=1)
		axis("off")
	end
end

makequiver(4)
savefig(joinpath(figurepath, "figure1g.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

makequiver(6)
savefig(joinpath(figurepath, "figure1h.pdf"), dpi=300, bbox_inches="tight", transparent=true)
plt[:close]()

ext: jl


text:
__precompile__(true)
module ContinuousPOMDPTreeSearchExperiments

importall POMDPs
using BasicPOMCP
using POMCPOW
using Parameters
using LightDarkPOMDPs
using StaticArrays
using POMDPToolbox
using ParticleFilters
using ControlSystems
using Distributions
using LaserTag
# import DESPOT
using ARDESPOT
using QMDP
using MCTS
using DiscreteValueIteration
using RecipesBase
using SubHunt

using DataFrames
using DataArrays
using ProgressMeter
using CPUTime

# using VDPTag
using VDPTag2
const VDPTag = VDPTag2 # yikes!!

# import DESPOT: bounds, init_bounds
import ARDESPOT: bounds, init_bounds
import POMCPOW.n_children

export
    RadiusRandom,
    AdaptiveRadiusRandom,
    SimpleFeedback,
    SymmetricNormalResampler,
    MinPopResampler,
    ObsAdaptiveParticleFilter,
    ObsAdaptiveSRFilter,
    LightDarkLQRSolver,
    LightDarkLQRPolicy,
    ModeAugmentedBelief,
    InfoGatherHeur,
    InfoGatherUpdater,
    OneStepValue,
    LaserBounds,
    InevitableInit,
    RootToNextMLFirst,
    NoGapTag,
    ModeMDP,
    ModeMDPSolver,
    nogap_tag,
    VDPUpper,
    SampleRollout,
    MeanRewardBeliefMDP,

    Sim,
    # SimQueue,
    run_parallel,
    run,
    metadata_as_pairs,
    Staged,
    LQG1D,
    LinearFeedback,
    Kalman1D,

    SimpleLightDark,
    DSimpleLightDark,
    LDHeuristic,
    LDHSolver,
    LDSide,
    LDSidePolicy,

    PDPWSolver,

    GBMDPSolver

include("simple_lightdark.jl")
include("lqg_example.jl")

include("policies.jl")
include("updaters.jl")
include("belief_mdp.jl")
include("action_gen.jl")

include("sr_filter.jl")
include("bounds_and_heuristics.jl")

# include("simulations.jl")
include("integration.jl")
include("pomcpdpw.jl")

n_children(h::BeliefNode) = length(h.children)

immutable OneStepValue end
BasicPOMCP.estimate_value(o::OneStepValue, pomdp::POMDP, s, h, steps) = reward(pomdp, s)

Base.srand(p::QMDPPolicy, i::Int) = p

end # module


text:
# Model hyperparameters

"""
RAPIDS Docs for KMeans: https://docs.rapids.ai/api/cuml/stable/api.html#k-means-clustering

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)

model = KMeans()
mach = machine(model, X)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct KMeans <: MMI.Unsupervised
    handle = nothing
    n_clusters::Int64 = 8::(_ > 0)
    max_iter::Int64 = 300::(_ > 0)
    tol::Float64 = 1e-4::(_ > 0)
    verbose::Bool = false
    random_state::Int = 1::(_ > 0)
    init::String = "scalable-k-means++"::(_ in ("scalable-k-means++","k-means||", "random"))
    n_init::Int = 1::(_ > 0)
    oversampling_factor::Float64 = 2.0::(_ > 0)
    max_samples_per_batch::Int64 = 32768::(_ > 0)
end

"""
RAPIDS Docs for DBSCAN: https://docs.rapids.ai/api/cuml/stable/api.html#dbscan

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)

model = DBSCAN()
mach = machine(model, X)
fit!(mach)
preds = mach.report.labels #DBSCAN does not have a predict method
```
"""
MLJModelInterface.@mlj_model mutable struct DBSCAN <: MMI.Unsupervised
    handle = nothing
    eps::Float64 = 1e-4::(_ > 0)
    min_samples::Int = 1::(_ > 0)
    metric::String = "euclidean"::(_ in ("euclidean","precomputed"))
    verbose::Bool = false
    max_mbytes_per_batch = nothing
    calc_core_sample_indices::Bool = true
end

"""
RAPIDS Docs for AgglomerativeClustering: https://docs.rapids.ai/api/cuml/stable/api.html#agglomerative-clustering

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)

model = AgglomerativeClustering()
mach = machine(model, X)
fit!(mach)
preds = mach.report.labels #AgglomerativeClustering does not have a predict method
```
"""
MLJModelInterface.@mlj_model mutable struct AgglomerativeClustering <: MMI.Unsupervised
    handle = nothing
    verbose::Bool = false
    affinity::String = "euclidean"::(_ in ("euclidean", "l1", "l2", "manhattan", "cosine"))
    # linkage::String = "single"::(_ in ("single", ))
    # only single linkage is supported. Error when specifying linkage
    n_neighbors::Int = 15::(_ > 0)
    connectivity::String = "knn"::(_ in ("knn", "pairwise"))
end

"""
RAPIDS Docs for HDBSCAN: https://docs.rapids.ai/api/cuml/stable/api.html#hdbscan

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)

model = HDBSCAN()
mach = machine(model, X)
fit!(mach)
preds = mach.report.labels #AgglomerativeClustering does not have a predict method
```
"""
MLJModelInterface.@mlj_model mutable struct HDBSCAN <: MMI.Unsupervised
    handle = nothing
    alpha::Float64 = 1.0::(_ > 0)
    verbose::Bool = false
    min_cluster_size::Int = 5::(_ > 0)
    min_samples = nothing
    cluster_selection_epsilon::Float64 = 0.0::(_ >= 0)
    max_cluster_size::Int = 0::(_ >= 0)
    # TODO: Why are we getting an affinity error when specifying this parameter
    # metric = "minkowski"
    p::Int = 2::(_ > 0)
    cluster_selection_method::String = "eom"::(_ in ("eom", "leaf"))
    allow_single_cluster::Bool = false
    # requires `hdbscan` cpu python package
    # gen_min_span_tree::Bool = false
    # gen_condensed_tree::Bool = false
    # gen_single_linkage_tree_::Bool = false
end

# Multiple dispatch for initializing models
model_init(mlj_model::KMeans) = cuml.KMeans(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::DBSCAN) = cuml.DBSCAN(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::AgglomerativeClustering) = cuml.AgglomerativeClustering(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::HDBSCAN) = cuml.HDBSCAN(; mlj_to_kwargs(mlj_model)...)

const CUML_CLUSTERING = Union{KMeans, DBSCAN, AgglomerativeClustering, HDBSCAN}

# add metadata
MMI.load_path(::Type{<:KMeans}) = "$PKG.KMeans"
MMI.load_path(::Type{<:DBSCAN}) = "$PKG.DBSCAN"
MMI.load_path(::Type{<:AgglomerativeClustering}) = "$PKG.AgglomerativeClustering"
MMI.load_path(::Type{<:HDBSCAN}) = "$PKG.HDBSCAN"

MMI.input_scitype(::Type{<:CUML_CLUSTERING}) = Union{AbstractMatrix, Table(Continuous)}

MMI.docstring(::Type{<:KMeans}) = "cuML's KMeans: https://docs.rapids.ai/api/cuml/stable/api.html#k-means-clustering"
MMI.docstring(::Type{<:DBSCAN}) = "cuML's DBSCAN: https://docs.rapids.ai/api/cuml/stable/api.html#dbscan"
MMI.docstring(::Type{<:AgglomerativeClustering}) = "cuML's AgglomerativeClustering: https://docs.rapids.ai/api/cuml/stable/api.html#agglomerative-clustering"
MMI.docstring(::Type{<:HDBSCAN}) = "cuML's HDBSCAN: https://docs.rapids.ai/api/cuml/stable/api.html#hdbscan"

# fit methods
function MMI.fit(mlj_model::CUML_CLUSTERING, verbosity, X, w=nothing)
    # fit the model
    model = model_init(mlj_model)
    model.fit(prepare_input(X))
    fitresult = model

    # save result
    cache = nothing
    report = ()
    return (fitresult, cache, report)
end

# predict methods
function MMI.predict(mlj_model::KMeans, fitresult, Xnew)
    model  = fitresult
    py_preds = model.predict(prepare_input(Xnew))
    preds = pyconvert(Array, py_preds)

    return preds
end

function MMI.predict(mlj_model::DBSCAN, fitresult, Xnew)
    @error "DBSCAN does not support predictions on new observations."
end

function MMI.predict(mlj_model::AgglomerativeClustering, fitresult, Xnew)
    @error "AgglomerativeClustering does not support predictions on new observations."
end

function MMI.predict(mlj_model::HDBSCAN, fitresult, Xnew)
    @error "HDBSCAN does not support predictions on new observations."
end

# Clustering metadata
MMI.metadata_pkg.((KMeans, DBSCAN, AgglomerativeClustering, HDBSCAN),
    name = "cuML Clustering Methods",
    uuid = "2764e59e-7dd7-4b2d-a28d-ce06411bac13", # see your Project.toml
    url  = "https://github.com/tylerjthomas9/RAPIDS.jl",  # URL to your package repo
    julia = false,          # is it written entirely in Julia?
    license = "MIT",        # your package license
    is_wrapper = true,      # does it wrap around some other package?
)


text:
"""
The examples module demonstrates how Wallace can be used to solve a number of
example problems. Each method within this module is named after the problem
being solved (using snake case), and returns a composed algorithm definition
which can be executed using the `run!` method.
"""
module example
  using core, common, algorithm, logger, criterion, evaluator, replacement,
        crossover, mutation, selection, fitness, population, breeder,
        initialiser, species, representation, koza, utility

  # Load each of the examples.
  include("example/one_max.jl")
  include("example/rastrigin.jl")
  include("example/symbolic.jl")
  include("example/tsp.jl")
end


text:
v1 = ConstrainedResampling(
    TruncateStd(1.5),
    NoConstraint(),
    [i % 2 == 0 ? TruncateQuantiles(0.1, 0.1*i) : TruncateStd(0.1*i) for i = 2:10])

@test v1 isa ConstrainedResampling
@test length(v1) == 3
@test v1[1] |> typeof <: SamplingConstraint


text:
module Features
import Combinatorics: combinations

import JuliaRL.FeatureCreators

struct RandomFeatures{T, F<:AbstractFloat} <: FeatureCreators.AbstractFeatureCreator
    # feature_map::Array{Array{Float64, 1}, 1}
    feature_map::Dict{T, Array{F, 1}}
    size::Int64
    active::Int64
    function RandomFeatures{F}(state_list::Array{T, 1}, feature_size, active_features; rng=Random.GLOBAL_RNG) where {T, F<:AbstractFloat}
        comb = shuffle(rng, collect(combinations(1:feature_size, active_features)))
        features = [zeros(F, feature_size) for i in 1:chain_size]
        for (idx, state) in enumerate(state_list)
            features[state[comb[idx]]] .= 1.0f0
        end
        new(features, feature_size, active_features)
    end
end

@forward RandomFeatures.features Base.getindex

FeatureCreators.create_features(fc::RandomFeatures, s; kwargs...) = fc[s]
FeatureCreators.feature_size(fc::RandomFeatures) = fc.size

end


text:
using DataFrames
using HDF5, JLD

DATA_DIR = expanduser("~/Dropbox/Data/lastfm-dataset-360K/")
PROFILES_FILE = DATA_DIR * "usersha1-profile.tsv"
FULL_DATASET_FILE = DATA_DIR * "usersha1-artmbid-artname-plays.tsv"
SHORT_DATASET_FILE = DATA_DIR * "short.tsv"
ARTIST_FILE = DATA_DIR * "artists.jld"
USER_FILE = DATA_DIR * "users.jld"
MODEL_FILE = DATA_DIR * "model.jld"
TOPN = 1000
LOW_THRESHOLD = 10

# read full last.fm dataset (17M records)
function read_full_dataset()
    return readtable(FULL_DATASET_FILE, separator='\t', quotemark=Char[],
                     names=[:userid, :artmbid, :artname, :plays])
end

# read top 1K artist dataset (8M)
# unlike read_full_dataset(), this returns preprocessed data
# where artname and user id are replaced with corresponding
# continuous ids
function read_short_dataset()
    return readtable(SHORT_DATASET_FILE, separator='\t')
end

# read user profile data
function read_profiles()
    return readtable(PROFILES_FILE, separator='\t', quotemark=Char[],
                     names=[:userid, :gender, :age, :country, :regdate])
end

max_stars_repo_name: dfdx/lastfm-rbm

max_issues_repo_name: dfdx/lastfm-rbm

max_forks_repo_name: dfdx/lastfm-rbm


text:
module RasterIO
    using Compat
    import Compat.String
    using BinDeps
    const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    if isfile(depfile)
      include(depfile)
    else
      error("libgdal not properly installed. Please run Pkg.build(\"RasterIO\")")
    end

    include("gdal_api.jl")
    include("raster_types.jl")
    include("raster_functions.jl")

    GDALAllRegister() # register all known drivers

end


text:
# Autogenerated wrapper script for librttopo_jll for i686-linux-gnu
export librttopo

using GEOS_jll
JLLWrappers.@generate_wrapper_header("librttopo")
JLLWrappers.@declare_library_product(librttopo, "librttopo.so.1")
function __init__()
    JLLWrappers.@generate_init_header(GEOS_jll)
    JLLWrappers.@init_library_product(
        librttopo,
        "lib/librttopo.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export Spectral_GapR
using JOLI, SeisJOLI

###### Calculate the spectral gap ratio
using LinearAlgebra, Arpack
function calculate_SGR(Mask)
    #Input: the mask in the reconstructed domain
    #Output: the spectral gap ratio for the input mask
	Z = svds(Mask, nsv = 2);

    return Z[1].S[2]/Z[1].S[1], rank(Mask)
end

function Spectral_GapR(nR::Int64,locs::Vector)
    # Inputs:
    # nR:the number of sources
    # locs: subsampling positions
    # Output: Spectral ratio

    # construct mask
    nC = nR;
    Mask = zeros(nC,nR);

    #apply reciprocity
    Mask[76:225,locs] .= 1;
    Mask[locs,76:225] .= 1;

    # Transfer data from SR-> MO domain
    SR = joSRtoCMO(nC,nR,DDT=Float64);
    RecMask = reshape(SR*vec(Mask),nC,2*nR-1);

    #return Spectral ratio
    return calculate_SGR(RecMask)
end


text:
"""
    ExplicitRungeKuttaScheme(
        step_size::Float64,
        alpha::LinearAlgebra.LowerTriangular{Float64},
        beta::LinearAlgebra.LowerTriangular{Float64},
    )

Generic data container for ExplicitRungeKuttaScheme for use in `integrate_time`.

We are interested in numerically integrating ODEs of the form
    ``a'(t) = a(t) Q``

where ``a(t)`` is a vector of coefficients and ``Q`` is a matrix, given an initial condition ``a(0)``.

In this context, Explicit Runge-Kutta methods with ``s`` stages can be written as

``v^{(0)} = a(t),``

``v^{(l)} = \\sum_{k=0}^{l-1} \\alpha_{l k} v^{(k)} + h\\beta_{l k}  v^{(k)} Q,\\, l = 1,...,s,``

``a(t+h) =v^{(s)},``

where ``h`` is the step-size and ``\\alpha_{l k}`` and ``\\beta_{l k}`` are parameters which define the scheme.

# Arguments
-  `step_size::Float64`: the step size of the integration scheme.
-  `alpha::LinearAlgebra.LowerTriangular{Float64}`: coefficients as described above
-  `beta::LinearAlgebra.LowerTriangular{Float64}`: coefficients as described above
"""
struct ExplicitRungeKuttaScheme
    step_size::Float64
    alpha::LinearAlgebra.LowerTriangular{Float64}
    beta::LinearAlgebra.LowerTriangular{Float64}
    function ExplicitRungeKuttaScheme(step_size::Float64,
        alpha::LinearAlgebra.LowerTriangular{Float64},
        beta::LinearAlgebra.LowerTriangular{Float64})

        t1 = (step_size > 0)
        checksquare(alpha)
        l1 = size(alpha)
        l2 = size(beta)
        (!t1)&&throw(DomainError("step_size must be positive"))
        !(l1==l2)&&throw(DimensionMismatch("alpha, beta must have same size"))
        return new(step_size,alpha,beta)
    end
end

"""
    ForwardEuler(step_size::Float64) <: ExplicitRungeKuttaScheme

Defines an Euler integration scheme to be used in `integrate_time`.

# Arguments
-  `step_size::Float64`: the step size of the integration scheme.
"""
ForwardEuler(step_size::Float64) = ExplicitRungeKuttaScheme(
        step_size,
        LinearAlgebra.LowerTriangular([1.0][:,:]),
        LinearAlgebra.LowerTriangular([1.0][:,:])
    )

"""

    Heuns(step_size::Float64) <: ExplicitRungeKuttaScheme

Defines Heuns integration scheme to be used in `integrate_time`.

# Arguments
-  `step_size::Float64`: the step size of the integration scheme.
"""
Heuns(step_size::Float64) = ExplicitRungeKuttaScheme(
        step_size,
        LinearAlgebra.LowerTriangular([1.0 0.0; 0.5 0.5]),
        LinearAlgebra.LowerTriangular([1.0 0.0; 0.0 0.5])
    )

"""

    StableRK3(step_size::Float64) <: ExplicitRungeKuttaScheme

Defines a strong stability preserving Runge-Kutta integration scheme to be used in `integrate_time`.

# Arguments
-  `step_size::Float64`: the step size of the integration scheme.

Hesthaven, J. S. & Warburton, T. (2007), Nodal discontinuous Galerkin methods: algorithms, analysis, and applications, Springer Science & Business Media.
(Section 5.7)

"""
StableRK3(step_size::Float64) = ExplicitRungeKuttaScheme(
        step_size,
        LinearAlgebra.LowerTriangular([1.0 0.0 0.0; 0.75 0.25 0.0; 1/3 0.0 2/3]),
        LinearAlgebra.LowerTriangular([1.0 0.0 0.0; 0.0 0.25 0.0; 0.0 0.0 2/3])
    )

_α = [  1.0                0.0                 0.0                 0.0                 0.0                 ;
        0.44437049406734   0.55562950593266    0.0                 0.0                 0.0                 ;
        0.62010185138540   0.0                 0.37989814861460    0.0                 0.0                 ;
        0.17807995410773   0.0                 0.0                 0.82192004589227    0.0                 ;
        0.00683325884039   0.0                 0.51723167208978    0.12759831133288    0.34833675773694    ]
_β = [  0.39175222700392    0.0                 0.0                 0.0                 0.0                 ;
        0.0                 0.36841059262959    0.0                 0.0                 0.0                 ;
        0.0                 0.0                 0.25189177424738    0.0                 0.0                 ;
        0.0                 0.0                 0.0                 0.54497475021237    0.0                 ;
        0.0                 0.0                 0.0                 0.08460416338212    0.22600748319395    ]
"""
    StableRK4(step_size::Float64) = ExplicitRungeKuttaScheme(
            step_size,
            LinearAlgebra.LowerTriangular(_α),
            LinearAlgebra.LowerTriangular(_β)
        )

Defines a strong stability preserving Runge-Kutta integration scheme to be used in `integrate_time`.

# Arguments
-  `step_size::Float64`: the step size of the integration scheme.

Raymond J. Spiteri and Steven J. Ruuth. A new class of optimal high-order
strong-stability-preserving time discretization methods. SIAM J. Numer. Anal.,
40(2):469-491, 2002.

Hesthaven, J. S. & Warburton, T. (2007), Nodal discontinuous Galerkin methods: algorithms, analysis, and applications, Springer Science & Business Media.
(Section 5.7)
"""
StableRK4(step_size::Float64) = ExplicitRungeKuttaScheme(
        step_size,
        LinearAlgebra.LowerTriangular(_α),
        LinearAlgebra.LowerTriangular(_β)
    )

"""
Given `x0` and `D` apprximate `x0 exp(Dy)`.

    integrate_time(x0::Array{Float64,2}, D::AbstractArray{Float64,2},
        y::Float64, scheme::ExplicitRungeKuttaScheme [; limiter])

# Arguments
- `x0`: An initial vector
- `D`: A square matrix
- `y`: time to integrate up to
- `h`: ExplicitRungeKuttaScheme.
- `limiter`: optional named argument spacifying the slope-limiter to used. The
    limiter function will be applied at every stage and every time-step of the
    ExplicitRungeKuttaScheme.
"""
function integrate_time(x0::Array{Float64,2}, D::AbstractArray{Float64,2},
    y::Float64, scheme::ExplicitRungeKuttaScheme)

    checksquare(D)
    !(size(x0,2)==size(D,1))&&throw(DimensionMismatch("x0 must have length size(D,1)"))

    return _integrate(x0[:],D,y,scheme)
end
function integrate_time(x0::Array{Float64,1}, D::AbstractArray{Float64,2},
    y::Float64, scheme::ExplicitRungeKuttaScheme)

    checksquare(D)

    return _integrate(x0,D,y,scheme)
end

function integrate_time(x0::SFMDistribution, D::AbstractArray{Float64,2},
    y::Float64, scheme::ExplicitRungeKuttaScheme)
    return SFMDistribution(integrate_time(x0.coeffs,D,y,scheme),x0.dq)
end

function integrate_time(x0::SFMDistribution{DGMesh{T}}, D::AbstractArray{Float64,2},
    y::Float64, scheme::ExplicitRungeKuttaScheme; limiter::Limiter=GeneralisedMUSCL) where T

    checksquare(D)

    limiter_params = limiter.generate_params(x0.dq)
    limiter_function = x->limiter.fun(x,limiter_params...)

    return SFMDistribution(_integrate(x0.coeffs,D,y,scheme,limiter_function),x0.dq)
end

"""
    _integrate(x0::Array{Float64,1},
        D::Union{Array{Float64,2},SparseArrays.SparseMatrixCSC{Float64,Int}},
        y::Float64, scheme::ExplicitRungeKuttaScheme)

Use ExplicitRungeKuttaScheme method.
"""
function _integrate(x0::Array{Float64,1}, D::AbstractArray{Float64,2},
    y::Float64, scheme::ExplicitRungeKuttaScheme)

    return _integrate(x0, D, y, scheme, identity)
end

function _integrate(x0::Array{Float64,1}, D::AbstractArray{Float64,2},
    y::Float64, scheme::ExplicitRungeKuttaScheme, limit_function::Function)

    x = limit_function(x0)
    h = scheme.step_size
    l = size(scheme.alpha,1)
    α = scheme.alpha
    βh = scheme.beta*h
    v = Array{Float64,2}(undef,length(x0),l+1)
    vD = Array{Float64,2}(undef,length(x0),l+1)
    for t = h:h:y
        v[:,1] = x
        for i in 1:l
            vD[:,i] = transpose(v[:,i])*D
            initialised = false
            for j in 1:i
                if α[i,j]!=0.0
                    initialised ? (v[:,i+1]+=v[:,j]*α[i,j]) : (v[:,i+1]=v[:,j]*α[i,j]; initialised=true)
                end
                if βh[i,j]!=0.0
                    initialised ? (v[:,i+1]+=vD[:,j]*βh[i,j]) : (v[:,i+1]=vD[:,j]*βh[i,j]; initialised=true)
                end
            end
            v[:,i+1] = limit_function(v[:,i+1])
        end
        x = v[:,end]
    end
    return x
end


text:
# Control whether timings are enabled or not, by default yes
# Note: TimerOutputs is not thread-safe, so do not use `@timeit`
# or `@timing` in threaded regions unless you know what you are doing.
if get(ENV, "DFTK_TIMING", "1") == "1"
    timer_enabled() = true
else
    timer_enabled() = false
end

"""TimerOutput object used to store DFTK timings."""
const timer = TimerOutput()

"""
Shortened version of the `@timeit` macro from `TimerOutputs`,
which writes to the DFTK timer.
"""
macro timing(args...)
    if DFTK.timer_enabled()
        TimerOutputs.timer_expr(__module__, false, :($(DFTK.timer)), args...)
    else  # Disable taking timings
        :($(esc(last(args))))
    end
end


text:
__precompile__()

module EisensteinNumbers

import Base: +, -, *, /, \, real, imag, complex, Complex, isreal, conj, inv, abs2, show, zero, one, iszero

"Error when finding the inverse of zero."
const ZeroInverse = "inverse of zero"

"Error when the denominator in a quotient is zero."
const ZeroDenominator = "denominator is zero"

# Half of the square root of 3
Base.@irrational halfsqrt3 0.86602540378443864676 (sqrt(big(3)) / 2)

include("Eisenstein.jl")

export Eisenstein, Eisenstein16, Eisenstein32, Eisenstein64, Eisenstein128, Eisenstein256, EisensteinBigInt
export unreal, random, associates, ω, omega, asarray, halfsqrt3

end


text:
module Utils

export save_global!, print_banner, prefix_tuple_names, select_device!
export SimulationProgressMessenger
export fit_cubic
export InterpolatedProfile, InterpolatedProfileTimeSeries, ∂z, ∂t
export diagnose_buoyancy_flux

using Printf

using CUDA: CuDevice, device!

using Oceananigans.Simulations
using Oceananigans.Diagnostics
using Oceananigans.Utils

function select_device!(ndev)
    @show dev = CuDevice(ndev)
    device!(dev)
    return nothing
end

function save_global!(file, group, name)
    val = try
        eval(name)
    catch err
        elaboration = err isa UndefVarError ? "because $(err.var) is not defined in the global scope." : ""
        @warn "$name will not be saved $elaboration"
        nothing
    end

    if val !== nothing
        file["$group/$name"] = val
    end

    return nothing
end

function print_banner(simulation)

    model = simulation.model

    banner = """

    Simulation with

                 Nx, Ny: $(model.grid.Nx)
                     Nz: $(model.grid.Nz)
                 Lx, Ly: $(model.grid.Lx) meters
                     Lz: $(model.grid.Lz) meters
                 output: $(simulation.output_writers[:fields].filepath)
    """
        #averages output: $(simulation.output_writers[:averages].filepath)

    println(banner)

    return nothing
end

function prefix_tuple_names(prefix, tup)
    new_keys = (Symbol(prefix, :_, key) for key in keys(tup))
    return (; zip(new_keys, values(tup))...)
end

include("progress_messenger.jl")
include("fit_cubic.jl")
include("interpolated_profiles.jl")
include("diagnose_buoyancy_flux.jl")

end # module


text:
include("CauchyWeight.jl")
include("Geometry.jl")
include("evaluation.jl")
include("skewProductFun.jl")
include("convolutionProductFun.jl")

# GreensFun

export GreensFun

# In GreensFun, K must be <: BivariateFun, since a kernel
# could consist of a Vector of ProductFun's and LowRankFun's.
# If there are GreensFun's in the Vector, then we recurse.

immutable GreensFun{K<:BivariateFun,T} <: BivariateFun{T}
    kernels::Vector{K}
    function GreensFun(Kernels::Vector{K})
        @assert all(map(domain,Kernels).==domain(Kernels[1]))
        if any(K->K<:GreensFun,map(typeof,Kernels))
            return GreensFun(vcat(map(kernels,Kernels)...))
        end
        new(Kernels)
    end
end
GreensFun{K<:BivariateFun}(kernels::Vector{K}) = GreensFun{eltype(kernels),mapreduce(eltype,promote_type,kernels)}(kernels)

GreensFun{K<:BivariateFun}(F::K) = GreensFun(K[F])

Base.length(G::GreensFun) = length(G.kernels)
Base.transpose(G::GreensFun) = GreensFun(mapreduce(transpose,+,G.kernels))
Base.convert(::Type{GreensFun},F::Union{ProductFun,LowRankFun}) = GreensFun(F)
Base.rank(G::GreensFun) = error("Not all kernels are low rank approximations.")

domain(G::GreensFun) = domain(first(G.kernels))
@compat (G::GreensFun)(x,y)=evaluate(G,x,y)
evaluate(G::GreensFun,x,y) = mapreduce(f->evaluate(f,x,y),+,G.kernels)
kernels(B::BivariateFun) = B
kernels(G::GreensFun) = G.kernels

Base.rank{L<:LowRankFun}(G::GreensFun{L}) = mapreduce(rank,+,G.kernels)
slices{L<:LowRankFun}(G::GreensFun{L}) = mapreduce(x->x.A,vcat,G.kernels),mapreduce(x->x.B,vcat,G.kernels)
slices{L<:LowRankFun}(G::GreensFun{L},k::Int) = slices(G)[k]
LowRankIntegralOperator{L<:LowRankFun}(G::GreensFun{L}) = LowRankIntegralOperator(slices(G)...)

Base.promote_rule{K,T,K1,T1}(::Type{GreensFun{K,T}},::Type{GreensFun{K1,T1}}) = GreensFun{promote_type(K,K1),promote_type(T,T1)}

defaultgetindex(⨍::Operator,G::GreensFun) = mapreduce(f->getindex(⨍,f),+,G.kernels)

function getindex{F<:BivariateFun}(⨍::DefiniteLineIntegral,B::Matrix{F})
    m,n = size(B)
    wsp = domainspace(⨍)
    @assert m == length(wsp.spaces)
    ⨍j = DefiniteLineIntegral(wsp[1])
    ret = Array(Operator{promote_type(eltype(⨍j),map(eltype,B)...)},m,n)
    for j=1:n
        ⨍j = DefiniteLineIntegral(wsp[j])
        for i=1:m
            ret[i,j] = ⨍j[B[i,j]]
        end
    end
    ops=promotespaces(ret)
    InterlaceOperator(ops,
                      PiecewiseSpace(map(domainspace,ret[1,:])),
                      PiecewiseSpace(map(rangespace,ret[:,1])))
end

# Algebra with BivariateFun's

+(F::GreensFun,G::GreensFun) = GreensFun([F.kernels;G.kernels])
-(F::GreensFun,G::GreensFun) = GreensFun([F.kernels;-G.kernels])
+(G::GreensFun,B::BivariateFun) = GreensFun([G.kernels;kernels(B)])
-(G::GreensFun,B::BivariateFun) = GreensFun([G.kernels;-kernels(B)])
+(B::BivariateFun,G::GreensFun) = GreensFun([kernels(B);G.kernels])
-(B::BivariateFun,G::GreensFun) = GreensFun([kernels(B);-G.kernels])

# Custom operations on Arrays required to infer type of resulting Array{GreensFun}

for op in (:+,:-)
    @eval begin
        function $op{F<:BivariateFun,T1,G<:BivariateFun,T2}(A::Array{GreensFun{F,T1}},B::Array{GreensFun{G,T2}})
            C = similar(A, GreensFun{promote_type(F,G),promote_type(T1,T2)}, promote_shape(size(A),size(B)))
            for i=1:length(A)
                @inbounds C[i] = $op(A[i],B[i])
            end
            return C
        end
    end
end

## TODO: Get ProductFun & LowRankFun in different CauchyWeight spaces to promote to GreensFun.
#=
+{S<:UnivariateSpace,V<:UnivariateSpace,O1,O2}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,CauchyWeight{O2}}) = GreensFun([F,G])
+{S<:UnivariateSpace,V<:UnivariateSpace,O1,SS}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,SS}) = GreensFun([F,G])
+{S<:UnivariateSpace,V<:UnivariateSpace,SS<:AbstractProductSpace,O1}(F::ProductFun{S,V,SS},G::ProductFun{S,V,CauchyWeight{O1}}) = GreensFun([F,G])

-{S<:UnivariateSpace,V<:UnivariateSpace,O1,O2}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,CauchyWeight{O2}}) = GreensFun([F,-G])
-{S<:UnivariateSpace,V<:UnivariateSpace,O1,SS}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,SS}) = GreensFun([F,-G])
-{S<:UnivariateSpace,V<:UnivariateSpace,SS<:AbstractProductSpace,O1}(F::ProductFun{S,V,SS},G::ProductFun{S,V,CauchyWeight{O1}}) = GreensFun([F,-G])
=#

## Constructors

function GreensFun{SS<:AbstractProductSpace}(f::Function,ss::SS;method::Symbol=:lowrank,kwds...)
    if method == :standard
        F = ProductFun(f,ss,kwds...)
    elseif method == :convolution
        F = convolutionProductFun(f,ss,kwds...)
    elseif method == :unsplit
        # Approximate imaginary & smooth part.
        F1 = skewProductFun((x,y)->imag(f(x,y)),ss.space;kwds...)
        # Extract diagonal value.
        d = domain(ss)
        xm,ym = mean([first(d[1]),last(d[1])]),mean([first(d[2]),last(d[2])])
        F1m = F1(xm,ym)
        # Set this normalized part to be the singular part.
        F1 = ProductFun(-coefficients(F1)/F1m/2,ss)
        # Approximate real & smooth part after singular extraction.
        m,n = size(F1)
        if typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Chebyshev,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n)+1)
        elseif typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Laurent,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n))
        end
        F = [F1;F2]
    elseif method == :lowrank
        F = LowRankFun(f,ss;method=:standard,kwds...)
    elseif method == :Cholesky
        F = LowRankFun(f,ss;method=method,kwds...)
    end
    GreensFun(F)
end

function GreensFun{SS<:AbstractProductSpace}(f::Function,g::Function,ss::SS;method::Symbol=:unsplit,kwds...)
    if method == :unsplit
        # Approximate Riemann function of operator.
        G = skewProductFun(g,ss.space;kwds...)
        # Normalize and set to the singular part.
        F1 = ProductFun(-coefficients(G)/2,ss)
        # Approximate real & smooth part after singular extraction.
        m,n = size(F1)
        if typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Chebyshev,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n)+1)
        elseif typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Laurent,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n))
        end
        F = [F1;F2]
    end
    GreensFun(F)
end

# Array of GreensFun on TensorSpace of PiecewiseSpaces

function GreensFun{PWS1<:PiecewiseSpace,PWS2<:PiecewiseSpace}(f::Function,ss::AbstractProductSpace{Tuple{PWS1,PWS2}};method::Symbol=:lowrank,tolerance::Symbol=:absolute,kwds...)
    M,N = length(ss[1]),length(ss[2])
    @assert M == N
    G = Array(GreensFun,N,N)
    if method == :standard
        for i=1:N,j=1:N
            G[i,j] = GreensFun(f,ss[i,j];method=method,kwds...)
        end
    elseif method == :convolution
        for i=1:N,j=i:N
            G[i,j] = GreensFun(f,ss[i,j];method=method,kwds...)
        end
        for i=2:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    elseif method == :unsplit
        maxF = Array(Number,N)
        for i=1:N
          G[i,i] = GreensFun(f,ss[i,i];method=method,kwds...)
          maxF[i] = one(real(mapreduce(eltype,promote_type,G[i,i].kernels)))/2π
        end
        for i=1:N,j=i+1:N
            G[i,j] = GreensFun(f,ss[i,j].space;method=:lowrank,tolerance=(tolerance,max(maxF[i],maxF[j])),kwds...)
        end
        for i=2:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    elseif method == :lowrank
        for i=1:N,j=1:N
            G[i,j] = GreensFun(f,ss[i,j];method=method,kwds...)
        end
    elseif method == :Cholesky
        if tolerance == :relative
            for i=1:N
                G[i,i] = GreensFun(f,ss[i,i];method=method,kwds...)
                for j=i+1:N
                    G[i,j] = GreensFun(f,ss[i,j];method=:lowrank,kwds...)
                end
            end
        elseif tolerance == :absolute
            maxF = Array(Number,N)
            for i=1:N
                F,maxF[i] = LowRankFun(f,ss[i,i];method=method,retmax=true,kwds...)
                G[i,i] = GreensFun(F)
            end
            for i=1:N,j=i+1:N
                G[i,j] = GreensFun(f,ss[i,j];method=:lowrank,tolerance=(tolerance,max(maxF[i],maxF[j])),kwds...)
            end
        end
        for i=2:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    end
    G
end

function GreensFun{PWS1<:PiecewiseSpace,PWS2<:PiecewiseSpace}(f::Function,g::Function,ss::AbstractProductSpace{Tuple{PWS1,PWS2}};method::Symbol=:unsplit,tolerance::Symbol=:absolute,kwds...)
    M,N = length(ss[1]),length(ss[2])
    @assert M == N
    G = Array(GreensFun,N,N)
    if method == :unsplit
        maxF = Array(Number,N)
        for i=1:N
            G[i,i] = GreensFun(f,g,ss[i,i];method=method,kwds...)
            maxF[i] = one(real(mapreduce(eltype,promote_type,G[i,i].kernels)))/2π
        end
        for i=1:N
            for j=i+1:N
                G[i,j] = GreensFun(f,ss[i,j].space;method=:lowrank,tolerance=(tolerance,max(maxF[i],maxF[j])),kwds...)
            end
        end
        for i=1:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    end
    G
end

# HierarchicalMatrix of GreensFuns on TensorSpace of HierarchicalSpaces

converttoPiecewiseSpace(H::Space) = H
converttoPiecewiseSpace(H::HierarchicalSpace) = PiecewiseSpace(H)

function partition{HS1<:HierarchicalSpace,HS2<:HierarchicalSpace}(ss::ApproxFun.TensorSpace{Tuple{HS1,HS2}})
    ss11,ss12 = partition(ss[1])
    ss21,ss22 = partition(ss[2])
    (ss11⊗ss21,ss12⊗ss22),(converttoPiecewiseSpace(ss12)⊗converttoPiecewiseSpace(ss21),converttoPiecewiseSpace(ss11)⊗converttoPiecewiseSpace(ss22))
end

function partition{O,HS1<:HierarchicalSpace,HS2<:HierarchicalSpace}(ss::CauchyWeight{O,Tuple{HS1,HS2}})
    ss11,ss12 = partition(ss[1])
    ss21,ss22 = partition(ss[2])
    (CauchyWeight(ss11⊗ss21,O),CauchyWeight(ss12⊗ss22,O)),(CauchyWeight(converttoPiecewiseSpace(ss12)⊗converttoPiecewiseSpace(ss21),O),CauchyWeight(converttoPiecewiseSpace(ss11)⊗converttoPiecewiseSpace(ss22),O))
end

function GreensFun{HS1<:HierarchicalSpace,HS2<:HierarchicalSpace}(f::Function,ss::AbstractProductSpace{Tuple{HS1,HS2}};method::Symbol=:lowrank,kwds...)
    (ss11,ss22),(ss21,ss12) = partition(ss)
    meth1 = method == :Cholesky || method == :lowrank ? :standard : method
    G11 = GreensFun(f,ss11;method=method,kwds...)
    G22 = GreensFun(f,ss22;method=method,kwds...)
    G21 = GreensFun(LowRankFun(f,ss21;method=meth1,kwds...))
    G12 = method == :Cholesky ? transpose(G21) : GreensFun(LowRankFun(f,ss12;method=meth1,kwds...))
    return HierarchicalMatrix((G11,G22),(G21,G12))
end

blocksize{F<:GreensFun,G<:GreensFun}(H::HierarchicalMatrix{F,G}) = map(length,domain(H).domains)

function domain{F<:GreensFun,G<:GreensFun}(H::HierarchicalMatrix{F,G})
    H11,H22 = diagonaldata(H)
    H21,H12 = offdiagonaldata(H)
    m1,n2 = domain(H12)[1],domain(H12)[2]
    m2,n1 = domain(H21)[1],domain(H21)[2]
    @assert (m1,n1) == (domain(H11)[1],domain(H11)[2])
    @assert (m2,n2) == (domain(H22)[1],domain(H22)[2])
    (m1∪m2)*(n1∪n2)
end

function Base.getindex{G<:GreensFun,L<:LowRankFun,T}(⨍::DefiniteLineIntegral,H::HierarchicalMatrix{G,GreensFun{L,T}})
    H11,H22 = diagonaldata(H)
    wsp = domainspace(⨍)
    if length(domain(H11)[2]) ≥ 2
        ⨍1 = DefiniteLineIntegral(PiecewiseSpace(wsp[1:length(domain(H11)[2])]))
    else
        ⨍1 = DefiniteLineIntegral(wsp[1])
    end
    if length(domain(H22)[2]) ≥ 2
        ⨍2 = DefiniteLineIntegral(PiecewiseSpace(wsp[end-length(domain(H22)[2])+1:end]))
    else
        ⨍2 = DefiniteLineIntegral(wsp[end])
    end
    HierarchicalOperator((qrfact(⨍1[H11]),qrfact(⨍2[H22])),map(LowRankIntegralOperator,offdiagonaldata(H)))
end

Base.qrfact(H::HierarchicalOperator) = H # trivial no-op for now.


text:
# Autogenerated wrapper script for Chuffed_jll for aarch64-linux-musl-cxx03
export fznchuffed

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Chuffed")
JLLWrappers.@declare_executable_product(fznchuffed)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        fznchuffed,
        "bin/fzn-chuffed",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    DebiasedWhittleLikelihood(model::Type{<:TimeSeriesModel}, ts, Δ; lowerΩcutoff, upperΩcutoff, taper)
    DebiasedWhittleLikelihood(model::Type{<:TimeSeriesModel}, timeseries::TimeSeries; lowerΩcutoff, upperΩcutoff, taper)

Generate a function to evaluate the Debiased Whittle likelihood it's gradient and expected Hessian.

Create a callable struct which prealloctes memory appropriately.

    (f::DebiasedWhittleLikelihood)(θ)

Evaluates the Whittle likelihood at θ.

    (f::DebiasedWhittleLikelihood)(F,G,EH,θ)

Evaluates the Whittle likelihood at θ and stores the gradient and expected Hessian in G and EH respectively.
If F, G or EH equal nothing, then the function, gradient or expected Hessian are not evaluated repsectively.

# Aruguments
- `model`: the model for the process. Should be of type TimeSeriesModel, so `OU` and not `OU(1,1)`.
- `ts`: the timeseries in the form of an n by d matrix (where d is the dimension of the time series model).
- `Δ`: the sampling rate of the time series.
- `timeseries`: can be provided instead of `ts` and `Δ`. Must be of type `TimeSeries`.
- `lowerΩcutoff`: the lower bound of the frequency range included in the likelihood.
- `upperΩcutoff`: the upper bound of the frequency range included in the likelihood.
- `taper`: optional taper which should be a vector of length `size(ts,1)`, or `nothing` in which case no taper will be used.

Note that to use the gradient the model must have `grad_add_sdf!` or `grad_acv!` specified.
Similarly, to use the Hessian, the model must have `hess_add_sdf!` or `hess_acv!` specified.

# Examples
```julia-repl
julia> obj = DebiasedWhittleLikelihood(OU, ones(1000), 1)
Debiased Whittle likelihood for the OU model.

julia> obj([1.0, 1.0])
-1982.0676530999626

julia> F, G, EH = 0.0, zeros(2), zeros(2,2)
(0.0, [0.0, 0.0], [0.0 0.0; 0.0 0.0])

julia> obj(F, G, EH, [1.0, 1.0])
-1982.0676530999626

julia> G
2-element Vector{Float64}:
 1998.3136810970122
 -685.7904154568779

julia> H
2×2 Matrix{Float64}:
 -0.00967179   0.0607696
  0.0607696   -0.381827
```
"""
struct DebiasedWhittleLikelihood{T,S<:TimeSeriesModelStorage,M}
    data::DebiasedWhittleData{T}
    memory::S
    function DebiasedWhittleLikelihood(model, timeseries::TimeSeries; lowerΩcutoff = 0, upperΩcutoff = Inf, taper = nothing)
        DebiasedWhittleLikelihood(model, timeseries.ts, timeseries.Δ; lowerΩcutoff = lowerΩcutoff, upperΩcutoff = upperΩcutoff, taper = taper)
    end
    function DebiasedWhittleLikelihood(
        model::Type{<:TimeSeriesModel{D}}, ts, Δ;
        lowerΩcutoff = 0, upperΩcutoff = Inf, taper = nothing) where {D}

        Δ > 0 || throw(ArgumentError("Δ should be a positive."))
        D == size(ts,2) || throw(ArgumentError("timeseries is $(size(ts,2)) dimensional, but model is $D dimensional."))
        if taper !== nothing
            taper .*= inv(sqrt(sum(abs2,taper))) # normalise the taper
        end
        wdata = DebiasedWhittleData(model, ts, Δ, lowerΩcutoff = lowerΩcutoff, upperΩcutoff = upperΩcutoff, taper = taper)
        mem = allocate_memory_EI_FG(model, size(ts,1), Δ, taper = taper)
        new{eltype(wdata.I),typeof(mem),model}(wdata,mem)
    end
end
(f::DebiasedWhittleLikelihood{T,S,M})(θ) where {T,S,M} = debiasedwhittle!(f.memory,M(θ),f.data)
function (f::DebiasedWhittleLikelihood{T,S,M})(F,G,EH,θ) where {T,S,M}
    checksize(G,EH,θ)
    debiasedwhittle_Fisher!(F,G,EH,f.memory,M(θ),f.data)
end
Base.show(io::IO, ::DebiasedWhittleLikelihood{T,S,M}) where {T,S,M} = print(io, "Debiased Whittle likelihood for the $(M) model.")
getmodel(::DebiasedWhittleLikelihood{T,S,M}) where {T,S,M} = M

## Interior functions

"""
    debiasedwhittle!(store::TimeSeriesModelStorage, model::TimeSeriesModel, data::GenWhittleData)

Function to compute the debiased Whittle likelihood using a preallocated store.
"""
function debiasedwhittle!(store::TimeSeriesModelStorage, model::TimeSeriesModel, data::GenWhittleData)
    EI!(store, model)
    return generalwhittle(store, data)
end

"""
    debiasedwhittle_FG!(F, G, store, model::TimeSeriesModel, data::GenWhittleData)

Compute the debiased Whittle likelihood and its gradient using a preallocated store.
"""
function debiasedwhittle_FG!(F, G, store, model::TimeSeriesModel ,data::GenWhittleData)
    if F !== nothing || G !== nothing
        EI!(store, model)
    end
    if G !== nothing
        grad_EI!(store, model)
        grad_generalwhittle!(G, store, data)
    end
    if F !== nothing
        return generalwhittle(store, data)
    end
    return nothing
end

"""
    debiasedwhittle_FGH!(F, G, H, store, model::TimeSeriesModel, data::GenWhittleData)

Compute the debiased Whittle likelihood and its gradient and hessian using a preallocated store.
"""
function debiasedwhittle_FGH!(F, G, H, store, model::TimeSeriesModel, data::GenWhittleData)
    if F !== nothing || G !== nothing || H !== nothing
        EI!(store, model)
    end
    if G !== nothing || H !== nothing
        grad_EI!(store, model)
    end
    if H !== nothing
        hess_EI!(store, model)
        hess_generalwhittle!(H, store, data)
    end
    if G !== nothing
        grad_generalwhittle!(G, store, data)
    end
    if F !== nothing
        return generalwhittle(store, data)
    end
    return nothing
end

"""
    getallderiv(store::AdditiveStorage)
    getallderiv(store::TimeSeriesModelStorageGradient)

Extract all of the derivatives from a generic store.
"""
function getallderiv(store::AdditiveStorage)
    Vcat(getallderiv(store.store1), getallderiv(store.store2))
end
getallderiv(store::TimeSeriesModelStorageGradient) = extract_array(store.gradmemory)

"""
    extract_array(store::Sdf2EIStorage)
    extract_array(store::Acv2EIStorage)
    extract_array(store::Sdf2EIStorageUni)
    extract_array(store::Acv2EIStorageUni)

Extract array from storage.
"""
extract_array(store::Sdf2EIStorage) = extract_array(store.acv2EI)
extract_array(store::Acv2EIStorage) = store.hermitianarray
extract_array(store::Sdf2EIStorageUni) = extract_array(store.acv2EI)
extract_array(store::Acv2EIStorageUni) = store.allocatedarray

"""
    debiasedwhittle_Ehess!(EH, store, data::GenWhittleData)

Compute the expected hessian of the de-biased Whittle likelihood.
"""
function debiasedwhittle_Ehess!(EH, store, data::GenWhittleData)
    EH .= zero(eltype(EH))
    ∇S = getallderiv(store)
    S = extract_array(extract_S(store))
    length(S) == size(∇S,2) || error("S should be same length as second dim of ∇S")
    @inbounds for iω = data.Ω_used_index
        invS = inv(S[iω]) # account for double frequency resolution
        for jpar = 1:size(EH, 1), kpar = 1:jpar
            EH[jpar,kpar] += real(tr(∇S[jpar, iω] * invS*∇S[kpar, iω] * invS))
        end
    end
    for jpar = 1:size(EH, 1)-1, kpar = jpar+1:size(EH, 1) # fill upper triangle from symmetry
        @views EH[jpar, kpar] = EH[kpar, jpar]
    end
    return nothing
end

"""
    debiasedwhittle_Fisher!(F, G, H, store, model::TimeSeriesModel, data::GenWhittleData)

Compute the debiased Whittle likelihood, its gradient and fisher information.
"""
function debiasedwhittle_Fisher!(F, G, H, store, model::TimeSeriesModel, data::GenWhittleData)
    if F !== nothing || G !== nothing || H !== nothing
        EI!(store, model)
    end
    if G !== nothing || H !== nothing
        grad_EI!(store, model)
    end
    if H !== nothing
        debiasedwhittle_Ehess!(H, store, data)
    end
    if G !== nothing
        grad_generalwhittle!(G, store, data)
    end
    if F !== nothing
        return generalwhittle(store, data)
    end
    return nothing
end


text:
function Base.convert(
    ::Type{Distributions.GeneralizedPareto{T}}, d::GeneralizedPareto{S}
) where {T<:Real,S<:Real}
    return Distributions.GeneralizedPareto{T}(T(d.μ), T(d.σ), T(d.k))
end
function Base.convert(
    ::Type{Distributions.GeneralizedPareto}, d::GeneralizedPareto{T}
) where {T<:Real}
    return convert(Distributions.GeneralizedPareto{T}, d)
end


text:
minmaxrescale(x::AbstractArray{T}, minval, maxval) where {T<:Real} = minval .+ (
    maxval - minval) .* (x .-minimum(x)) ./ (maximum(x) - minimum(x))


text:
# Random Variable Application #

"Post-projection application"
function ppapl end

"Apply function to argument"
function apl end

"Apply if randvar otherwise just return value (treat as constant randvar)"
subapl(x, ω) = x
@inline subapl(rv::RandVar, ω) = apl(rv, ω)

"Project `ω` to `x`"
proj(ω::ΩBase, x::RandVar) = proj(ω, x.id, 1)

proj(ω::O, i...) where {I, O <: ΩBase{I}} =
  ΩProj{O, I}(ω, base(I, i...))

proj(tω::TaggedΩ, x::RandVar) = tag(proj(tω.taggedω, x), tω.tags)
@spec _res.tags == tω.tags "tags are preserved in projection"

"Project `ω` to `rv` then apply"
@inline apl(rv::RandVar, ω::ΩBase) =  ppapl(rv, proj(ω, rv))

"Reproject back to parent random variable"
@inline apl(rv::RandVar, πω::ΩProj) = apl(rv, parentω(πω))

"Reify arguments (resolve random variables to values)"
@inline reify(ω, args) = map(x -> subapl(x, ω), args)
@spec all([r isa elemtype(a) for (a, r) in zip(args, _res)])

@inline reify(ω, (x1, x2)::Tuple{T1, T2}) where {T1, T2} =
 subapl(x1, ω), subapl(x2, ω)
@inline reify(ω, (x1, x2, x3)::Tuple{T1, T2, T3}) where {T1, T2, T3} =
 subapl(x1, ω), subapl(x2, ω), subapl(x3, ω)
@inline reify(ω, (x1, x2, x3, x4)::Tuple{T1, T2, T3, T4}) where {T1, T2, T3, T4} =
 subapl(x1, ω), subapl(x2, ω), subapl(x3, ω), subapl(x4, ω)

# TODO use generated function to avoid runtime iteration in reify

# "Reify random variable args.. i.e. map(x -> apl(x, ω), args)"
# @generated function reify(ω, args)
#   if any(isa.(args RandVar))
#     map(t -> t isa RandVar, args)
#     quote
#       (apl(x, ))
#     end
#   else
#     quote
#       args
#     end
#   end
# end

@inline apl(rv::RandVar, tω::TaggedΩ{I, T, ΩT}) where {I, T, ΩT <: ΩProj}  =
  apl(rv, TaggedΩ(parentω(tω.taggedω), tω.tags))

# Generated funtion for typed dispatch on different tags (might be unnecssary)
@generated function apl(rv::RandVar, tω::TaggedΩ{I, Tags{K, V}, ΩT}) where {I, K, V, ΩT <: ΩBase}
  if hastags(tω, :replmap)
    :(Omega.Causal.replaceapl(rv, tω))
  elseif hastags(tω, :cache)
    :(memapl(rv, tω))
  else
    :(ppapl(rv, proj(tω, rv)))
  end
end


text:
function Base.get{S<:State1D,D,I,R}(::Feature_PosFt, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(0.0) # no lateral offset
end
function Base.get{S<:State1D,D,I,R}(::Feature_PosFyaw, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(0.0) # no lane-relative heading
end
function Base.get{S<:State1D,D,I,R}(::Feature_Speed, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(rec[pastframe][vehicle_index].state.v)
end
function Base.get{S<:State1D,D,I,R}(::Feature_VelFs, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(rec[pastframe][vehicle_index].state.v)
end
function Base.get{S<:State1D,D,I,R}(::Feature_VelFt, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(0.0) # no lateral speed
end