text:
const USE_GPU = true
const VISUALIZE = false
const BENCHMARK = true
include("elastic_wave_3D.jl")

# Grid sizes for the scaling experiment
grid_size = 16 * 2 .^ (1:4)

# Array allocation
T_effs = zeros(4)

# Scaling experiment loop
for it = 1:4
    T_effs[it] = elastic_wave_3D(grid_size[it], 8192.0 / 8^it)
end

# Visualization of Results
gr()
max_bandwith = 360.0 # Theoretical Max Memory Bandwidth of NVIDIA GeForce RTX 3060

plot(grid_size, max_bandwith .* ones(4), label="Theoretical Max Bandwidth", lw=3)
plot!(grid_size, T_effs, xlabel="grid_size", ylabel="T_eff GB/s", ylims=(0.0, 400.0), title="GPU Scaling Experiment", label="elastic_wave_3D.jl", lw=3)
png("./docs/img/elastic_wave_3D_scaling_experiment_gpu")


text:
using FileIO
using Test

@testset "Basic checks" begin
  fname = tempname()*".gslib"

  prop1 = rand(10,10,10)
  prop2 = rand(10,10,10)

  save(fname, [prop1,prop2])
  grid = load(fname)
  @test grid[:prop1] == prop1
  @test grid[:prop2] == prop2

  save(fname, grid)
  grid = load(fname)
  @test grid[:prop1] == prop1
  @test grid[:prop2] == prop2

  save(fname, prop1)
  grid = load(fname)
  @test grid[:prop1] == prop1

  rm(fname)
end

ext: jl


text:
using FFTW: r2r, r2r!, R2HC, HC2R
using Random

@doc raw"""
    surrogate(A::AbstractArray{Float64, 1})

Generates a new surrogate time-series for the given one by means of phase-randomization.
"""
function surrogate(A::AbstractArray{Float64, 1})
    N = length(A)
    # compute real -> half-complex DFT
    res = r2r(A, R2HC, 1)
    # shuffle phases in half-complex representation
    ϕ = rand(Float64, N÷2)
    for i in 1:(N÷2)
        z = res[1+i] + 1im*res[end-i];
        z = abs(z)*(exp(2im*π*ϕ[i]))
        res[1+i] = real(z)
        res[end-i] = imag(z)
    end
    # compute inverse half-complex -> real DFT inplace
    r2r!(res, HC2R, 1)
    return res ./ N;
end

@doc raw"""
    surrogate(A::AbstractArray{Float64, 1})

Generates a new surrogate time-series in-place for the given one by means of phase-randomization.
"""
function surrogate!(A::AbstractArray{Float64, 1})
    N = length(A)
    # compute real -> half-complex DFT
    r2r!(A, R2HC, 1)
    # shuffle phases in half-complex representation
    ϕ = rand(Float64, N÷2)
    for i in 1:(N÷2)
        z = A[1+i] + 1im*A[end-i];
        z = abs(z)*(exp(2im*π*ϕ[i]))
        A[1+i] = real(z)
        A[end-i] = imag(z)
    end
    # compute inverse half-complex -> real DFT inplace
    r2r!(A, HC2R, 1)
    A ./= N
end

@doc raw"""
    surrogate(A::AbstractArray{Float64, 1})

Generates a new surrogate time-series for each column of `A` by means of phase-randomization.
"""
function surrogate(A::AbstractArray{Float64, 2})
    N,M = size(A)
    # compute real -> half-complex DFT
    res = r2r(A, R2HC, 1)
    # shuffle phases in half-complex representation
    ϕ = rand(Float64, M)
    for i in 1:(N÷2)
        z = res[1+i,:] + 1im*res[end-i,:];
        z = map(abs,z).*(exp.(2im*π*ϕ))
        res[1+i,:] = map(real,z)
        res[end-i,:] = map(imag,z)
    end
    # compute inverse half-complex -> real DFT inplace
    r2r!(res, HC2R, 1)
    return res ./ N;
end

@doc raw"""
    surrogate(A::AbstractArray{Float64, 1})

Generates a new surrogate time-series in-place for each column of `A` by means of phase-randomization.
"""
function surrogate!(A::AbstractArray{Float64, 2})
    N,M = size(A)
    # compute real -> half-complex DFT inplace
    r2r!(A, R2HC, 1)
    # shuffle phases in half-complex representation
    ϕ = rand(Float64, M)
    for i in 1:(N÷2)
        z = A[1+i,:] + 1im*A[end-i,:];
        z = map(abs,z) .* (exp.(2im*π*ϕ))
        A[1+i,:] = map(real,z)
        A[end-i,:] = map(imag,z)
    end
    # compute inverse half-complex -> real DFT inplace
    r2r!(A, HC2R, 1)
    A ./= N
    return A;
end


text:
using Documenter, PowerModels

makedocs(
    modules = [PowerModels],
    format = :html,
    sitename = "PowerModels",
    authors = "Carleton Coffrin, Russell Bent, and contributors.",
    analytics = "UA-367975-10",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Getting Started" => "quickguide.md",
            "Network Data Format" => "network-data.md",
            "Result Data Format" => "result-data.md"
        ],
        "Library" => [
            "Network Formulations" => "formulations.md",
            "Problem Specifications" => "specifications.md",
            "Modeling Components" => [
                "PowerModel" => "model.md",
                "Objective" => "objective.md",
                "Variables" => "variables.md",
                "Constraints" => "constraints.md"
            ],
            "Relaxation Schemes" => "relaxations.md",
            "File IO" => "parser.md"
        ],
        "Developer" => "developer.md",
        "Experiment Results" => "experiment-results.md"
    ]
)

deploydocs(
    deps = nothing,
    make = nothing,
    target = "build",
    repo = "github.com/lanl-ansi/PowerModels.jl.git",
    julia = "0.6"
)


text:
#!/usr/bin/env julia

#Start Test Script
using Latexify
using LaTeXStrings
using ParameterizedFunctions
using Missings
using Test

# Run tests

@testset "macro test" begin include("macros.jl") end
@testset "recipe test" begin include("recipe_test.jl") end
@testset "latexify tests" begin include("latexify_test.jl") end
@testset "latexraw tests" begin include("latexraw_test.jl") end
@testset "latexalign tests" begin include("latexalign_test.jl") end
@testset "latexarray tests" begin include("latexarray_test.jl") end
@testset "latexequation tests" begin include("latexequation_test.jl") end
@testset "latexinline tests" begin include("latexinline_test.jl") end
@testset "latextabular tests" begin include("latextabular_test.jl") end
@testset "mdtable tests" begin include("mdtable_test.jl") end
@testset "chemical_arrows test" begin include("chemical_arrows_test.jl") end
@testset "DataFrame Plugin" begin include("plugins/DataFrames.jl") end
@testset "unocode2latex" begin include("unicode2latex.jl") end
@testset "cdot test" begin include("cdot_test.jl") end
@testset "numberformatters" begin include("numberformatters_test.jl") end


text:
Base.@kwdef struct DGPMBeta <: AbstractModel
    # Data
    y0::Vector{Float64}
    X0::Matrix{Float64}
    y1::Vector{Float64}
    X1::Matrix{Float64}
    # Hyperparameters
    a_v0::Float64 = 1.0
    b_v0::Float64 = 1.0
    # Parameters
    μ::Vector{Float64} = rand(1)
    v::Vector{Float64} = rand(InverseGamma(a_v0, b_v0), 1)
    μ_shadow::Vector{Float64} = randexp(1) * 10.0
    v_shadow::Vector{Float64} = randexp(1) * 10.0
    # Transformed parameters
    sumlogy1::Vector{Float64} = zeros(1)
    sumlogy2::Vector{Float64} = zeros(1)
    # Skeleton
    skl::Skeleton = Skeleton(; y0, y1, X0, X1)
end

function skeleton(m::DGPMBeta)
    return m.skl
end

function kernel_pdf(m::DGPMBeta, y0::Float64, j::Int)
    (; μ, v) = m
    kernel = Beta(μ[j] * v[j], (1 - μ[j]) * v[j])
    return pdf(kernel, y0)
end

function update_atoms!(m::DGPMBeta)
    (; a_v0, b_v0, μ, v, μ_shadow, v_shadow, skl) = m
    n = cluster_sizes(skl)
    update_suffstats!(m)
    while length(μ) < rmax(skl)
        push!(μ, rand())
        push!(v, rand(InverseGamma(a_v0, b_v0)))
        push!(μ_shadow, 10 * randexp())
        push!(v_shadow, 10 * randexp())
    end
    s_μ = Walker2020Sampler(10.0, 0.0, 1.0);
    s_v = Walker2020Sampler(10.0, 0.0, Inf);
    for j in 1:rmax(skl)
        if n[j] > 0
            μ[j], μ_shadow[j] =
                rand(Random.GLOBAL_RNG, s_μ, x -> logpμ(m, x, j), μ[j], μ_shadow[j])
            v[j], v_shadow[j] =
                rand(Random.GLOBAL_RNG, s_v, x -> logpv(m, x, j), v[j], v_shadow[j])
        else
            μ[j] = rand()
            v[j] = rand(InverseGamma(a_v0, b_v0))
            μ_shadow[j] = 10 * randexp()
            v_shadow[j] = 10 * randexp()
        end
    end
end

function update_suffstats!(m::DGPMBeta)
    (; y0, sumlogy1, sumlogy2, skl) = m
    d = cluster_labels(skl)
    while length(sumlogy1) < rmax(skl)
        push!(sumlogy1, 0.0)
        push!(sumlogy2, 0.0)
    end
    for j in 1:rmax(skl)
        sumlogy1[j] = 0.0
        sumlogy2[j] = 0.0
    end
    for i in 1:length(y0)
        sumlogy1[d[i]] += log(y0[i])
        sumlogy2[d[i]] += log(1.0 - y0[i])
    end
    return nothing
end

function logpμ(m::DGPMBeta, μ0::Float64, j::Int)
    (; sumlogy1, sumlogy2, v, skl) = m
    n = cluster_sizes(skl)
    μv = μ0 * v[j]
    return (
        (μv - 1) * sumlogy1[j] +
        (v[j] - μv - 1) * sumlogy2[j] -
        n[j] * logbeta(μv, v[j] - μv)
    )
end

function logpv(m::DGPMBeta, v0::Float64, j::Int)
    (; sumlogy1, sumlogy2, a_v0, b_v0, μ, skl) = m
    n = cluster_sizes(skl)
    μv = μ[j] * v0
    return (
        (μv - 1) * sumlogy1[j] +
        (v0 - μv - 1) * sumlogy2[j] -
        n[j] * logbeta(μv, v0 - μv) -
        (a_v0 + 1) * log(v0) -
        b_v0 / v0
    )
end


text:
# ALL TESTS THAT CONCERN MODFICIATION

# include
using LatticePhysics
using PyPlot
using Base.Test

# begin the Modification testblock
bandstructure_testset = @testset "Bandstructures" begin

################################################################################
#
#   METHODS FOR CONSTRUCTION OF BANDSTRUCTURES (ALONG PATHS)
#
################################################################################

################################################################################
#
#   1) TYPE BANDSTRUCTURE
#       - type definition
#       - TODO printInfo function
#
################################################################################

################################################################################
#
#   2) CALCULATION OF BAND STRUCTURES OF UNTICELL OBJECTS
#
################################################################################

# begin the testset
@testset "Calculation of bandstructures" begin

    # begin the testset
    @testset "2D (honeycomb)" begin

        # get a unitcell
        unitcell = getUnitcellHoneycomb()
        # get a path
        path = getDefaultPathTriangular(resolution=30)

        # test the calculations (bare)
        @test typeof(getBandstructure(unitcell, path)) == Bandstructure

        # test the calculations (resolution set)
        @test typeof(getBandstructure(unitcell, path, resolution=60)) == Bandstructure

        # test the calculations (hermitian)
        @test typeof(getBandstructure(unitcell, path, enforce_hermitian=true)) == Bandstructure

    # end the testset here
    end;

    # begin the testset
    @testset "2D (square)" begin

        # get a unitcell
        unitcell = getUnitcellSquare()
        # get a path
        path = getDefaultPathSquare(resolution=30)

        # test the calculations (bare)
        @test typeof(getBandstructure(unitcell, path)) == Bandstructure

        # test the calculations (resolution set)
        @test typeof(getBandstructure(unitcell, path, resolution=60)) == Bandstructure

        # test the calculations (hermitian)
        @test typeof(getBandstructure(unitcell, path, enforce_hermitian=true)) == Bandstructure

    # end the testset here
    end;

    # begin the testset
    @testset "3D (diamond)" begin

        # get a unitcell
        unitcell = getUnitcellDiamond()
        # get a path
        path = getDefaultPathFCC(resolution=60)

        # test the calculations (bare)
        @test typeof(getBandstructure(unitcell, path)) == Bandstructure

        # test the calculations (resolution set)
        @test typeof(getBandstructure(unitcell, path, resolution=120)) == Bandstructure

        # test the calculations (hermitian)
        @test typeof(getBandstructure(unitcell, path, enforce_hermitian=true)) == Bandstructure

    # end the testset here
    end;

# end the testset here
end;

################################################################################
#
#   3) PLOTTING OF BAND STRUCTURES
#       - plotting of Bandstructure objects
#       - plotting of bandstructures of unitcells along paths
#
################################################################################

# begin the testset
@testset "Plotting of bandstructures" begin

    # begin the testset
    @testset "2D (honeycomb)" begin

        # get a unitcell
        unitcell = getUnitcellHoneycomb()
        # get a path
        path = getDefaultPathTriangular(resolution=30)

        # get the bandstructure (bare)
        bandstructure = getBandstructure(unitcell, path)

        # test the plotting
        @test_nowarn close(plotBandstructure(bandstructure, showPlot=false))
        @test_nowarn close(plotBandstructure(bandstructure, limits_energy=[-1, 2.0], showPlot=false))
        @test_nowarn close(plotBandstructure(bandstructure, plot_title="myplot", showPlot=false))
        @test_nowarn close(plotBandstructure(bandstructure, plot_color="r", showPlot=false))
        @test_nowarn close(plotBandstructure(bandstructure, figsize=(10,8), showPlot=false))
        @test_nowarn close(plotBandstructure(bandstructure, save_filename="$(FOLDER_SPECTRA)testplot.pdf", showPlot=false))

        # test calculation and plotting
        @test_nowarn close(plotBandstructure(unitcell, path, showPlot=false))
        @test_nowarn close(plotBandstructure(unitcell, path, limits_energy=[-1, 2.0], showPlot=false))
        @test_nowarn close(plotBandstructure(unitcell, path, plot_title="myplot", showPlot=false))
        @test_nowarn close(plotBandstructure(unitcell, path, plot_color="r", showPlot=false))
        @test_nowarn close(plotBandstructure(unitcell, path, figsize=(10,8), showPlot=false))
        @test_nowarn close(plotBandstructure(unitcell, path, save_filename="$(FOLDER_SPECTRA)testplot.pdf", showPlot=false))

    # end the testset here
    end;

    # begin the testset
    @testset "2D (square)" begin

        # get a unitcell
        unitcell = getUnitcellSquare()
        # get a path
        path = getDefaultPathSquare(resolution=30)

        # get the bandstructure (bare)
        bandstructure = getBandstructure(unitcell, path)

        # test the plotting
        @test_nowarn close(plotBandstructure(bandstructure, showPlot=false))

    # end the testset here
    end;

    # begin the testset
    @testset "3D (diamond)" begin

        # get a unitcell
        unitcell = getUnitcellDiamond()
        # get a path
        path = getDefaultPathFCC(resolution=60)

        # get the bandstructure (bare)
        bandstructure = getBandstructure(unitcell, path)

        # test the plotting
        @test_nowarn close(plotBandstructure(bandstructure, showPlot=false))

    # end the testset here
    end;

# end the testset here
end;

# end the testset here
end;


text:
using Agents
using AgentsPlots
using Random

mutable struct SugarSeeker <: AbstractAgent
	id::Int
	pos::Tuple{Int, Int}
	vision::Int
	metabolic_rate::Int
	age::Int
	max_age::Int
	wealth::Int
end

function distances(pos, sugar_peaks, max_sugar)
	all_dists = Array{Int, 1}(undef, length(sugar_peaks))
	for (ind, peak) in enumerate(sugar_peaks)
		d = round(Int, sqrt(sum((pos .- peak).^2)))
		all_dists[ind] = d
	end
	return minimum(all_dists)
end

function sugar_caps(dims, sugar_peaks, max_sugar, dia=4)
	sugar_capacities = zeros(Int, dims)
	for i in 1:dims[1], j in 1:dims[2]
		sugar_capacities[i, j] = distances((i,j), sugar_peaks, max_sugar)
	end
	for i in 1:dims[1]
		for j in 1:dims[2]
			sugar_capacities[i, j] = max(0, max_sugar - (sugar_capacities[i, j] ÷ dia))
		end
	end
	return sugar_capacities
end

"Start a sugarscape simulation"
function sugarscape(;dims=(50,50), sugar_peaks=((10, 40), (40, 10)), growth_rate=1, N=250, w0_dist=(5, 25), metabolic_rate_dist=(1,4), vision_dist=(1,6), 	max_age_dist=(60, 100), max_sugar = 4)
	sugar_capacities = sugar_caps(dims, sugar_peaks, max_sugar, 6)
	sugar_values = deepcopy(sugar_capacities)
	space = GridSpace(dims, periodic=true, moore=true)
	properties = Dict(
		:growth_rate => growth_rate,
		:N => N,
		:w0_dist => w0_dist,
		:metabolic_rate_dist => metabolic_rate_dist,
		:vision_dist => vision_dist,
		:max_age_dist => max_age_dist,
		:sugar_values => sugar_values,
		:sugar_capacities => sugar_capacities
		)
	model = AgentBasedModel(SugarSeeker, space, scheduler=random_activation,
		properties=properties)
	for ag in 1:N
		add_agent_single!(model,
			rand(vision_dist[1]:vision_dist[2]),
			rand(metabolic_rate_dist[1]:metabolic_rate_dist[2]),
			0,
			rand(max_age_dist[1]:max_age_dist[2]),
			rand(w0_dist[1]:w0_dist[2]),
		)
	end
	return model, sugarscape_agent_step!, sugarscape_env!
end

function sugarscape_env!(model)
	## At each cell, sugar grows back at a rate of $\alpha$ units per time-step up to the cell's capacity c.
	togrow = findall(x -> model.sugar_values[x] < model.sugar_capacities[x], 1:prod(model.space.dimensions))
	model.sugar_values[togrow] .+= model.growth_rate
end

function movement!(agent, model)
	posvertex = coord2vertex(agent.pos, model)
	newsite = posvertex
	## find all unoccupied cells within vision
	neighbors = node_neighbors(coord2vertex(agent.pos, model), model, agent.vision)
	empty_nodes = [i for i in neighbors if isempty(i, model)]
	if length(empty_nodes) > 0
		## identify the one(s) with greatest amount of sugar
		maxsugar = maximum(model.sugar_values[empty_nodes])
		if maxsugar > 0
		sugary_sites_inds = findall(x -> x == maxsugar, model.sugar_values[empty_nodes])
		sugary_sites = empty_nodes[sugary_sites_inds]
		## select the nearest one (randomly if more than one)
		for dia in 1:agent.vision
			nn = node_neighbors(posvertex, model, dia)
			suitable = intersect(nn, sugary_sites)
			if length(suitable) > 0
				newsite = rand(suitable)
				break
			end
		end
		## move there and collect all the sugar in it
		newsite != posvertex && move_agent!(agent, newsite, model)
		end
	end
	## update wealth (collected - consumed)
	agent.wealth += (model.sugar_values[newsite] - agent.metabolic_rate)
	model.sugar_values[newsite] = 0
	## age
	agent.age += 1
end

function replacement!(agent, model)
	## If the agent's sugar wealth become zero or less, it dies
	if agent.wealth <= 0 || agent.age >= agent.max_age
		kill_agent!(agent, model)
		## Whenever an agent dies, a young one is added to a random pos.
		## New agent has random attributes
		add_agent_single!(model,
			rand(model.vision_dist[1]:model.vision_dist[2]),
			rand(model.metabolic_rate_dist[1]:model.metabolic_rate_dist[2]),
			0,
			rand(model.max_age_dist[1]:model.max_age_dist[2]),
			rand(model.w0_dist[1]:model.w0_dist[2]),
		)
	end
end

function sugarscape_agent_step!(agent, model)
	movement!(agent, model)
	replacement!(agent, model)
end


text:
# Command Forward
"""
	forward!(session::Session)

This command causes the browser to traverse one step forwards in the joint session history of the current top-level browsing context.
"""
function forward!(session::Session)
    @unpack addr, id = session
    response =
        HTTP.post("$addr/session/$id/forward", [("Content-Type" => "application/json")])
    @assert response.status == 200
    nothing
end


text:
using LinearAlgebra
using Random
using NNlib: softplus
using Roots # for inverse

################################################################################
#                            Planar and Radial Flows                           #
#             Ref: Variational Inference with Normalizing Flows,               #
#               D. Rezende, S. Mohamed(2015) arXiv:1505.05770                  #
################################################################################

###############
# RadialLayer #
###############

mutable struct RadialLayer{T1 <: Real, T2 <: AbstractVector{<:Real}} <: Bijector{1}
    α_::T1
    β::T1
    z_0::T2
end

function RadialLayer(dims::Int, container=Array)
    α_ = randn()
    β = randn()
    z_0 = container(randn(dims))
    return RadialLayer(α_, β, z_0)
end

h(α, r) = 1 ./ (α .+ r)     # for radial flow from eq(14)
#dh(α, r) = .- (1 ./ (α .+ r)) .^ 2   # for radial flow; derivative of h()

# An internal version of transform that returns intermediate variables
function _transform(flow::RadialLayer, z::AbstractVecOrMat)
    return _radial_transform(flow.α_, flow.β, flow.z_0, z)
end
function _radial_transform(α_, β, z_0, z)
    α = softplus(α_)            # from A.2
    β_hat = -α + softplus(β)    # from A.2
    if z isa AbstractVector
        r = norm(z .- z_0)
    else
        r = vec(sqrt.(sum(abs2, z .- z_0; dims = 1)))
    end
    transformed = z .+ β_hat ./ (α .+ r') .* (z .- z_0)   # from eq(14)
    return (transformed = transformed, α = α, β_hat = β_hat, r = r)
end

(b::RadialLayer)(z::AbstractMatrix{<:Real}) = _transform(b, z).transformed
(b::RadialLayer)(z::AbstractVector{<:Real}) = vec(_transform(b, z).transformed)

function forward(flow::RadialLayer, z::AbstractVecOrMat)
    transformed, α, β_hat, r = _transform(flow, z)
    # Compute log_det_jacobian
    d = size(flow.z_0, 1)
    h_ = h(α, r)
    if transformed isa AbstractVector
        T = eltype(transformed)
    else
        T = typeof(vec(transformed))
    end
    log_det_jacobian::T = @. (
        (d - 1) * log(1 + β_hat * h_)
        + log(1 +  β_hat * h_ + β_hat * (- h_ ^ 2) * r)
    )   # from eq(14)
    return (rv = transformed, logabsdetjac = log_det_jacobian)
end

function (ib::Inverse{<:RadialLayer})(y::AbstractVector{<:Real})
    flow = ib.orig
    T = promote_type(eltype(flow.α_), eltype(flow.β), eltype(flow.z_0), eltype(y))
    TV = vectorof(T)
    α = softplus(flow.α_)            # from A.2
    β_hat = - α + softplus(flow.β)   # from A.2
    # Define the objective functional
    f(y) = r -> norm(y .- flow.z_0) - r * (1 + β_hat / (α + r))   # from eq(26)
    # Run solver
    rs::T = find_zero(f(y), zero(T), Order16())
    return (flow.z_0 .+ (y .- flow.z_0) ./ (1 .+ β_hat ./ (α .+ rs)))::TV
end
function (ib::Inverse{<:RadialLayer})(y::AbstractMatrix{<:Real})
    flow = ib.orig
    T = promote_type(eltype(flow.α_), eltype(flow.β), eltype(flow.z_0), eltype(y))
    TM = matrixof(T)
    α = softplus(flow.α_)            # from A.2
    β_hat = - α + softplus(flow.β)   # from A.2
    # Define the objective functional
    f(y) = r -> norm(y .- flow.z_0) - r * (1 + β_hat / (α + r))   # from eq(26)
    # Run solver
    rs = mapvcat(eachcol(y)) do c
        find_zero(f(c), zero(T), Order16())
    end
    return (flow.z_0 .+ (y .- flow.z_0) ./ (1 .+ β_hat ./ (α .+ rs')))::TM
end

logabsdetjac(flow::RadialLayer, x::AbstractVecOrMat) = forward(flow, x).logabsdetjac
isclosedform(b::Inverse{<:RadialLayer}) = false


text:
using Test, GpABC

@testset "Matern Kernel Tests" begin
    mi1 = ExponentialIsoKernel()
    mi3 = MaternIsoKernel(3)
    mi5 = MaternIsoKernel(5)

    ma1 = ExponentialArdKernel()
    ma3 = MaternArdKernel(3)
    ma5 = MaternArdKernel(5)

    R57 = [
        0.791406184803418	0.290820793849798	0.273243301209641	0.976630407840307	0.919637321922759	0.373113582084474	0.167134671084180;
        0.163747991045820	0.825082233559844	0.400163430441247	0.227264572040767	0.344469555814038	0.919974093753832	0.0324169118108307;
        0.112485507089871	0.157099440301349	0.546775426038902	0.635662289166666	0.499824566533707	0.00459372066175273	0.831036804630904;
        0.543917659457735	0.370293144774760	0.654585435825673	0.576094024151607	0.551489732019366	0.206663162078086	0.235575431601878;
        0.624102506145883	0.516884229386637	0.617419922080552	0.412263916066039	0.887926049478318	0.387765046389664	0.316278932763393;
    ]
    R55 = [
        0.700022954354073	0.281477719210776	0.145081796750930	0.243184048330031	0.327940463706823;
        0.285793828911052	0.796026775518445	0.389757992881238	0.194300282837032	0.695373520578897;
        0.336882508051061	0.822016138500542	0.662586354959098	0.604853641989107	0.142760369527508;
        0.984312814187404	0.625239778601211	0.0888677975576687	0.775160910922643	0.373530648822678;
        0.345468047701506	0.404958950544223	0.601073551737272	0.129429876864622	0.929614762742824;
    ]
    x1 = [
        0.289686565422307	0.400690696388805;
        0.762105052714285	0.408534206640224;
        0.699620968987989	0.234411002245923;
        0.529625781127078	0.745714747228527;
        0.420503402790289	0.227756669850907;
    ]
    x2 = [
        0.906350184188387	0.112500416632975;
        0.0275456542277456	0.794502350406024;
        0.959810605719597	0.201003274042627;
        0.536850329272845	0.262026834423666;
        0.843484088705612	0.885006158559120;
        0.437078875227402	0.298506903135454;
        0.984655946951665	0.00552099757153587;
    ]

    @testset "Check hyperparameters" begin
        @test 2 == get_hyperparameters_size(mi1, zeros(0, 0))
        @test 2 == get_hyperparameters_size(mi3, zeros(0, 0))
        @test 5 == get_hyperparameters_size(ma5, zeros(0, 4))
        @test 4 == get_hyperparameters_size(ma3, zeros(0, 3))
    end

    log_theta = log.([4, 2, 3])
    @testset "Matern ARD ν = 1 (Exponential ARD)" begin
        M = [
            11.5841611485303	13.2909793837280	11.3700092357433	14.0222586904381	11.6123040446907	14.7523433872054	11.0339312696684
            14.1591849555517	10.8419960392437	14.1814529341518	14.1516238340854	13.5805061233645	13.5446082674416	13.4388537336879
            14.3180788074261	10.8935494907303	14.0415131311797	14.7417626826706	12.7318338836143	14.0074644436926	13.6117512667987
            12.0575422606651	12.4413109254846	12.0746218530243	13.6169992534182	13.5845752777948	13.6877679310735	11.4383386327386
            12.5114734394976	12.1827554326725	12.2165024727785	15.0790185190317	11.7997688428430	15.6050010869727	11.9525830328299
        ]
        @test M ≈ covariance(ma1, log_theta, x1, x2)

        M = [
            16	12.6336641564393	12.9391758368811	13.5501472509310	14.6641249897463;
            12.6336641564393	16	14.9793620729469	13.6112655419223	13.3493771316819;
            12.9391758368811	14.9793620729469	16	13.2253470207306	13.9156250513432;
            13.5501472509310	13.6112655419223	13.2253470207306	16	13.3500485055730;
            14.6641249897463	13.3493771316819	13.9156250513432	13.3500485055730	16;
        ]
        @test M ≈ covariance_training(ma1, log_theta, x1)
        @test M ≈ covariance_training(ma1, log_theta, x1) # test hitting the cache
        dM = [343.921478384293, 10.3213449835906, 6.38933255398455]
        @test dM ≈ covariance_grad(ma1, log_theta, x1, R55)
        @test dM ≈ covariance_grad(ma1, log_theta, x1, R55)
    end

    @testset "Matern ARD ν = 3" begin
        M = [
            14.2606234566307	15.3314244048386	14.0932626383531	15.6407043245208	14.2820869992591	15.8558852837985	13.8164059103774;
            15.6882654783931	13.6506275988849	15.6956880091364	15.6857252510439	15.4649217521938	15.4491509635950	15.4014022972021;
            15.7393100068678	13.6956968927578	15.6475915262437	15.8534468545973	15.0333663530316	15.6353683076306	15.4784676598925;
            14.6050405864475	14.8579266649518	14.6168010730980	15.4807263285850	15.4666954161076	15.5107194381009	14.1474345463901;
            14.9015687011262	14.6901700399937	14.7126814838511	15.9212064638772	14.4218893190551	15.9854291910391	14.5317433595824;
        ]
        @test M ≈ covariance(ma3, log_theta, x1, x2)

        M = [
            16	14.9756414903414	15.1500190215362	15.4515988541525	15.8349379559396;
            14.9756414903414	16	15.9033310682782	15.4782583694961	15.3595059433015;
            15.1500190215362	15.9033310682782	16	15.2991709179202	15.6013855729400;
            15.4515988541525	15.4782583694961	15.2991709179202	16	15.3598254022276;
            15.8349379559396	15.3595059433015	15.6013855729400	15.3598254022276	16;
        ]
        @test M ≈ covariance_training(ma3, log_theta, x1)
        @test M ≈ covariance_training(ma3, log_theta, x1) # test hitting the cache
        dM = [372.084465282425, 4.78654375158504, 2.63809480215839]
        @test dM ≈ covariance_grad(ma3, log_theta, x1, R55)
        @test dM ≈ covariance_grad(ma3, log_theta, x1, R55)
    end

    @testset "Matern ARD ν = 5" begin
        M = [
            14.7343138231780	15.5570668558234	14.5972409963104	15.7722064253454	14.7517474101827	15.9127738901835	14.3662583343162;
            15.8040342966367	14.2255566733735	15.8089663802905	15.8023441542700	15.6515275009666	15.6404870077301	15.6068623435658;
            15.8377502311411	14.2639800921978	15.7768388221507	15.9112366875537	15.3387491805615	15.7686120399016	15.6609839364099;
            15.0098474591547	15.2060276946967	15.0190909685638	15.6625583124446	15.6527671021026	15.6833984211323	14.6418265876983;
            15.2393134357353	15.0764988027662	15.0940223135409	15.9533299666004	14.8644632222004	15.9916747989222	14.9519841042042;
        ]
        @test M ≈ covariance(ma5, log_theta, x1, x2)

        M = [
            16	15.2954029325564	15.4253279991210	15.6422028359281	15.8995197460673;
            15.2954029325564	16	15.9423599184652	15.6608380182498	15.5771213867627;
            15.4253279991210	15.9423599184652	16	15.5339169075397	15.7456140607006;
            15.6422028359281	15.6608380182498	15.5339169075397	16	15.5773489802123;
            15.8995197460673	15.5771213867627	15.7456140607006	15.5773489802123	16;
        ]
        @test M ≈ covariance_training(ma5, log_theta, x1)
        @test M ≈ covariance_training(ma5, log_theta, x1) # test hitting the cache
        dM = [374.899359986320, 3.42341729639196, 1.84924721006274]
        @test dM ≈ covariance_grad(ma5, log_theta, x1, R55)
        @test dM ≈ covariance_grad(ma5, log_theta, x1, R55)
    end

    log_theta = log.([3, 2])
    @testset "Matern ISO ν = 1 (Exponential ISO)" begin
        M = [
            6.40374962810124	7.10418580254636	6.34459259548824	7.81092023069155	6.22999187130739	8.22805896019984	6.03449962280165;
            7.63369048671289	5.94369110438302	7.79833100550643	7.86852381596106	7.06772605737364	7.58105534779673	7.14944854706554;
            7.98229092078633	5.81123101055548	7.89367475245528	8.28690232217822	6.44992796519164	7.86244889530660	7.49655498249347;
            6.22654782625103	6.99364663956730	6.36094189979183	7.06641758711903	7.58016585194165	7.16267020929581	5.82866729459857;
            7.01155324353779	6.37510026018038	6.87051722595412	8.47041933396979	6.08865655531065	8.67887226766782	6.64623836344120;
        ]
        @test M ≈ covariance(mi1, log_theta, x1, x2)

        M = [
            9.00000000000000	7.10630756931683	7.21409934498584	7.29433612106788	8.07525942876702;
            7.10630756931683	9.00000000000000	8.20486339969566	7.33344827035029	7.41853111640394;
            7.21409934498584	8.20486339969566	9.00000000000000	6.87446096867600	7.82736665827119;
            7.29433612106788	7.33344827035029	6.87446096867600	9.00000000000000	6.90717354045213;
            8.07525942876702	7.41853111640394	7.82736665827119	6.90717354045213	9.00000000000000;
        ]
        @test M ≈ covariance_training(mi1, log_theta, x1)
        @test M ≈ covariance_training(mi1, log_theta, x1) # test hitting the cache
        dM = [1.897023548035826e+02, 10.937241608256016]
        @test dM ≈ covariance_grad(mi1, log_theta, x1, R55)
        @test dM ≈ covariance_grad(mi1, log_theta, x1, R55)
    end

    @testset "Matern ISO ν = 3" begin
        M = [
            7.93393113569130	8.42245612528808	7.88637218712751	8.76945360658600	7.79153343442712	8.90205231203270	7.62163123534655;
            8.69668280017932	7.53928310045817	8.76460460450300	8.79101280527447	8.40050220472439	8.67321166928089	8.44916746012979;
            8.83055114163477	7.41533561669713	8.80010200913708	8.91632373211294	7.97038942178510	8.78878784805798	8.63376437028541;
            7.78862840945221	8.35469878668682	7.89961169913462	8.39970708049495	8.67280774994901	8.45685621289903	7.43190859244578;
            8.36591711705520	7.91101807687226	8.27504244294758	8.95369379680283	7.66971390180964	8.98290915814688	8.11878436808226;
        ]
        @test M ≈ covariance(mi3, log_theta, x1, x2)

        M = [
            9.00000000000000	8.42372172166323	8.48627246543244	8.53059901714994	8.85984522808927;
            8.42372172166323	9.00000000000000	8.89612286961489	8.55151073279344	8.59541763287696;
            8.48627246543244	8.89612286961489	9.00000000000000	8.27766161438813	8.77571419045471;
            8.53059901714994	8.55151073279344	8.27766161438813	9.00000000000000	8.29921452470236;
            8.85984522808927	8.59541763287696	8.77571419045471	8.29921452470236	9.00000000000000;
        ]
        @test M ≈ covariance_training(mi3, log_theta, x1)
        @test M ≈ covariance_training(mi3, log_theta, x1) # test hitting the cache

        dM = [2.075347189222474e+02, 5.640324374231643]
        @test dM ≈ covariance_grad(mi3, log_theta, x1, R55)
        @test dM ≈ covariance_grad(mi3, log_theta, x1, R55)
    end

    @testset "Matern ISO ν = 5" begin
        M = [
            8.21628350836750	8.60265235780503	8.17694984988844	8.85260451476692	8.09770244385133	8.94023400295971	7.95315483120493;
            8.80235090146831	7.88196229213989	8.84929839376336	8.86722556621832	8.58605965183515	8.78586698304239	8.62272975106313;
            8.89368929059150	7.77347792304969	8.87335031095893	8.94928811634894	8.24624704244772	8.86572266891085	8.75788234151801;
            8.09525833173716	8.55118379701449	8.18792735222480	8.58545717782324	8.78558219613035	8.62848598818880	7.78807358815230;
            8.55975743019400	8.19736780846164	8.48973373024520	8.97253479008701	7.99438732934831	8.99011500600917	8.36644127154376;
        ]
        @test M ≈ covariance(mi5, log_theta, x1, x2)

        M = [
            9.00000000000000	8.60360640239698	8.65041199909894	8.68315433586136	8.91297887699073;
            8.60360640239698	9.00000000000000	8.93644682441130	8.69847279148807	8.73035734382252;
            8.65041199909894	8.93644682441130	9.00000000000000	8.49176986934424	8.85686361374794;
            8.68315433586136	8.69847279148807	8.49176986934424	9.00000000000000	8.50848531932549;
            8.91297887699073	8.73035734382252	8.85686361374794	8.50848531932549	9.00000000000000;
        ]
        @test M ≈ covariance_training(mi5, log_theta, x1)
        @test M ≈ covariance_training(mi5, log_theta, x1) # test hitting the cache

        dM = [2.096069712802733e+02; 4.150924998668864]
        @test dM ≈ covariance_grad(mi5, log_theta, x1, R55)
        @test dM ≈ covariance_grad(mi5, log_theta, x1, R55)

    end

end


text:
#this method converts to vector of tuples from DataFrame without removing classification column
function ripserer_convert(data_frame::DataFrame)

    df_x = unique(data_frame)
    df_y = transpose(Array(df_x))
    df_z = [tuple(df_y[:,c]...) for c in 1:size(df_y, 2)]

    return df_z
end

#same as above except it removes the class col by entering its corresponding index
function ripserer_convert(data_frame::DataFrame, column_index::Int64)

    df_w = select!(data_frame, Not(column_index))

    df_z = ripserer_convert(df_w)

    return df_z
end

#this one doesn't work
function ripserer_convert(data_frame::DataFrame, column_index::AbstractString)

    df_w = select!(data_frame, Not(column_index))

    df_z = ripserer_convert(df_w)

    return df_z
end

#take in file handle
#Should I also make a scenario where the user enters a string giving the class col's name?

#converts from .csv file to vector of tuples and doesn't remove class col
#assumes the data has headers.
function ripserer_convert(csv_file::AbstractString; header_status::Bool = true)

    df = CSV.File(csv_file, delim = ",", header = header_status) |> DataFrame

    df_z = ripserer_convert(df)

    return df_z
end

#same as above but removes class col
function ripserer_convert(csv_file::AbstractString, column_index::Int64; header_status::Bool = true)

    df = CSV.File(csv_file, delim = ",", header = header_status) |> DataFrame

    df_w = select!(df, Not(column_index))
    df_z = ripserer_convert(df_w)

    return df_z
end


text:
# script file to reproduce benchmarks in the analysis at
# TODO http ref
using SearchBenchmarks
using Test

path = "" # working directory including path delimiter at end
source = "0x7f,123"

benchmark(200, 1000000, source, [naivesearch,bloom_v0], path*"(1) RND julia impl much better than naive search",false)
benchmark(200, 1000000, source, [bloom_v0,bloom_v1,bloom_v2,bloom_v3], path*"(2) RND julia bloom impl optimizations",false)
benchmark(200, 1000000, source, [bloom_v0,bloom_v3,bloom_v4], path*"(3) RND do bloom test first",false)
benchmark(2000, 1000000, source, [naivesearch,bloom_v0,bloom_v3,bloom_v4], path*"(4) RND degeneration on large patterns",false)
benchmark(2000, 1000000, source, [naivesearch,bloom_v0,bloom_v3,bloom_best], path*"(5) RND limited bits in bloom filter",false)
benchmark(200, 1000000, source, [bloom_v0,bloom_v3,bloom_best,bmsearch], path*"(6) RND boyer moore alternative",false)
benchmark(2000, 1000000, source, [naivesearch,bloom_v0,bloom_v3,bloom_best,bmsearch], path*"(7) RND boyer moore on large patterns",false)


text:
# ---
# notebook: true
# ---