text:
# Autogenerated wrapper script for FastJet_jll for armv7l-linux-gnueabihf-cxx11
export libClusteringVetoPlugin, libConstituentSubtractor, libEnergyCorrelator, libFlavorCone, libGenericSubtractor, libJetCleanser, libJetFFMoments, libJetsWithoutJets, libLundPlane, libNsubjettiness, libQCDAwarePlugin, libRecursiveTools, libScJet, libSoftKiller, libSubjetCounting, libValenciaPlugin, libVariableR, libfastjet, libfastjetplugins, libfastjettools, libsiscone, libsiscone_spherical

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `libClusteringVetoPlugin`
const libClusteringVetoPlugin_splitpath = ["lib", "libClusteringVetoPlugin.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libClusteringVetoPlugin_path = ""

# libClusteringVetoPlugin-specific global declaration
# This will be filled out by __init__()
libClusteringVetoPlugin_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libClusteringVetoPlugin = "libClusteringVetoPlugin.so"

# Relative path to `libConstituentSubtractor`
const libConstituentSubtractor_splitpath = ["lib", "libConstituentSubtractor.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libConstituentSubtractor_path = ""

# libConstituentSubtractor-specific global declaration
# This will be filled out by __init__()
libConstituentSubtractor_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libConstituentSubtractor = "libConstituentSubtractor.so"

# Relative path to `libEnergyCorrelator`
const libEnergyCorrelator_splitpath = ["lib", "libEnergyCorrelator.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libEnergyCorrelator_path = ""

# libEnergyCorrelator-specific global declaration
# This will be filled out by __init__()
libEnergyCorrelator_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libEnergyCorrelator = "libEnergyCorrelator.so"

# Relative path to `libFlavorCone`
const libFlavorCone_splitpath = ["lib", "libFlavorCone.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libFlavorCone_path = ""

# libFlavorCone-specific global declaration
# This will be filled out by __init__()
libFlavorCone_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libFlavorCone = "libFlavorCone.so"

# Relative path to `libGenericSubtractor`
const libGenericSubtractor_splitpath = ["lib", "libGenericSubtractor.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libGenericSubtractor_path = ""

# libGenericSubtractor-specific global declaration
# This will be filled out by __init__()
libGenericSubtractor_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libGenericSubtractor = "libGenericSubtractor.so"

# Relative path to `libJetCleanser`
const libJetCleanser_splitpath = ["lib", "libJetCleanser.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libJetCleanser_path = ""

# libJetCleanser-specific global declaration
# This will be filled out by __init__()
libJetCleanser_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libJetCleanser = "libJetCleanser.so"

# Relative path to `libJetFFMoments`
const libJetFFMoments_splitpath = ["lib", "libJetFFMoments.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libJetFFMoments_path = ""

# libJetFFMoments-specific global declaration
# This will be filled out by __init__()
libJetFFMoments_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libJetFFMoments = "libJetFFMoments.so"

# Relative path to `libJetsWithoutJets`
const libJetsWithoutJets_splitpath = ["lib", "libJetsWithoutJets.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libJetsWithoutJets_path = ""

# libJetsWithoutJets-specific global declaration
# This will be filled out by __init__()
libJetsWithoutJets_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libJetsWithoutJets = "libJetsWithoutJets.so"

# Relative path to `libLundPlane`
const libLundPlane_splitpath = ["lib", "libLundPlane.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libLundPlane_path = ""

# libLundPlane-specific global declaration
# This will be filled out by __init__()
libLundPlane_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libLundPlane = "libLundPlane.so"

# Relative path to `libNsubjettiness`
const libNsubjettiness_splitpath = ["lib", "libNsubjettiness.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libNsubjettiness_path = ""

# libNsubjettiness-specific global declaration
# This will be filled out by __init__()
libNsubjettiness_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libNsubjettiness = "libNsubjettiness.so"

# Relative path to `libQCDAwarePlugin`
const libQCDAwarePlugin_splitpath = ["lib", "libQCDAwarePlugin.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libQCDAwarePlugin_path = ""

# libQCDAwarePlugin-specific global declaration
# This will be filled out by __init__()
libQCDAwarePlugin_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libQCDAwarePlugin = "libQCDAwarePlugin.so"

# Relative path to `libRecursiveTools`
const libRecursiveTools_splitpath = ["lib", "libRecursiveTools.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libRecursiveTools_path = ""

# libRecursiveTools-specific global declaration
# This will be filled out by __init__()
libRecursiveTools_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libRecursiveTools = "libRecursiveTools.so"

# Relative path to `libScJet`
const libScJet_splitpath = ["lib", "libScJet.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libScJet_path = ""

# libScJet-specific global declaration
# This will be filled out by __init__()
libScJet_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libScJet = "libScJet.so"

# Relative path to `libSoftKiller`
const libSoftKiller_splitpath = ["lib", "libSoftKiller.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libSoftKiller_path = ""

# libSoftKiller-specific global declaration
# This will be filled out by __init__()
libSoftKiller_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libSoftKiller = "libSoftKiller.so"

# Relative path to `libSubjetCounting`
const libSubjetCounting_splitpath = ["lib", "libSubjetCounting.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libSubjetCounting_path = ""

# libSubjetCounting-specific global declaration
# This will be filled out by __init__()
libSubjetCounting_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libSubjetCounting = "libSubjetCounting.so"

# Relative path to `libValenciaPlugin`
const libValenciaPlugin_splitpath = ["lib", "libValenciaPlugin.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libValenciaPlugin_path = ""

# libValenciaPlugin-specific global declaration
# This will be filled out by __init__()
libValenciaPlugin_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libValenciaPlugin = "libValenciaPlugin.so"

# Relative path to `libVariableR`
const libVariableR_splitpath = ["lib", "libVariableR.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libVariableR_path = ""

# libVariableR-specific global declaration
# This will be filled out by __init__()
libVariableR_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libVariableR = "libVariableR.so"

# Relative path to `libfastjet`
const libfastjet_splitpath = ["lib", "libfastjet.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libfastjet_path = ""

# libfastjet-specific global declaration
# This will be filled out by __init__()
libfastjet_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libfastjet = "libfastjet.so.0"

# Relative path to `libfastjetplugins`
const libfastjetplugins_splitpath = ["lib", "libfastjetplugins.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libfastjetplugins_path = ""

# libfastjetplugins-specific global declaration
# This will be filled out by __init__()
libfastjetplugins_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libfastjetplugins = "libfastjetplugins.so.0"

# Relative path to `libfastjettools`
const libfastjettools_splitpath = ["lib", "libfastjettools.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libfastjettools_path = ""

# libfastjettools-specific global declaration
# This will be filled out by __init__()
libfastjettools_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libfastjettools = "libfastjettools.so.0"

# Relative path to `libsiscone`
const libsiscone_splitpath = ["lib", "libsiscone.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libsiscone_path = ""

# libsiscone-specific global declaration
# This will be filled out by __init__()
libsiscone_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libsiscone = "libsiscone.so.0"

# Relative path to `libsiscone_spherical`
const libsiscone_spherical_splitpath = ["lib", "libsiscone_spherical.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libsiscone_spherical_path = ""

# libsiscone_spherical-specific global declaration
# This will be filled out by __init__()
libsiscone_spherical_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libsiscone_spherical = "libsiscone_spherical.so.0"

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    global libClusteringVetoPlugin_path = normpath(joinpath(artifact_dir, libClusteringVetoPlugin_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libClusteringVetoPlugin_handle = dlopen(libClusteringVetoPlugin_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libClusteringVetoPlugin_path))

    global libConstituentSubtractor_path = normpath(joinpath(artifact_dir, libConstituentSubtractor_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libConstituentSubtractor_handle = dlopen(libConstituentSubtractor_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libConstituentSubtractor_path))

    global libEnergyCorrelator_path = normpath(joinpath(artifact_dir, libEnergyCorrelator_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libEnergyCorrelator_handle = dlopen(libEnergyCorrelator_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libEnergyCorrelator_path))

    global libFlavorCone_path = normpath(joinpath(artifact_dir, libFlavorCone_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libFlavorCone_handle = dlopen(libFlavorCone_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libFlavorCone_path))

    global libGenericSubtractor_path = normpath(joinpath(artifact_dir, libGenericSubtractor_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libGenericSubtractor_handle = dlopen(libGenericSubtractor_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libGenericSubtractor_path))

    global libJetCleanser_path = normpath(joinpath(artifact_dir, libJetCleanser_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libJetCleanser_handle = dlopen(libJetCleanser_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libJetCleanser_path))

    global libJetFFMoments_path = normpath(joinpath(artifact_dir, libJetFFMoments_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libJetFFMoments_handle = dlopen(libJetFFMoments_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libJetFFMoments_path))

    global libJetsWithoutJets_path = normpath(joinpath(artifact_dir, libJetsWithoutJets_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libJetsWithoutJets_handle = dlopen(libJetsWithoutJets_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libJetsWithoutJets_path))

    global libLundPlane_path = normpath(joinpath(artifact_dir, libLundPlane_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libLundPlane_handle = dlopen(libLundPlane_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libLundPlane_path))

    global libNsubjettiness_path = normpath(joinpath(artifact_dir, libNsubjettiness_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libNsubjettiness_handle = dlopen(libNsubjettiness_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libNsubjettiness_path))

    global libQCDAwarePlugin_path = normpath(joinpath(artifact_dir, libQCDAwarePlugin_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libQCDAwarePlugin_handle = dlopen(libQCDAwarePlugin_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libQCDAwarePlugin_path))

    global libRecursiveTools_path = normpath(joinpath(artifact_dir, libRecursiveTools_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libRecursiveTools_handle = dlopen(libRecursiveTools_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libRecursiveTools_path))

    global libScJet_path = normpath(joinpath(artifact_dir, libScJet_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libScJet_handle = dlopen(libScJet_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libScJet_path))

    global libSoftKiller_path = normpath(joinpath(artifact_dir, libSoftKiller_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libSoftKiller_handle = dlopen(libSoftKiller_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libSoftKiller_path))

    global libSubjetCounting_path = normpath(joinpath(artifact_dir, libSubjetCounting_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libSubjetCounting_handle = dlopen(libSubjetCounting_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libSubjetCounting_path))

    global libValenciaPlugin_path = normpath(joinpath(artifact_dir, libValenciaPlugin_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libValenciaPlugin_handle = dlopen(libValenciaPlugin_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libValenciaPlugin_path))

    global libVariableR_path = normpath(joinpath(artifact_dir, libVariableR_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libVariableR_handle = dlopen(libVariableR_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libVariableR_path))

    global libfastjet_path = normpath(joinpath(artifact_dir, libfastjet_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libfastjet_handle = dlopen(libfastjet_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libfastjet_path))

    global libfastjetplugins_path = normpath(joinpath(artifact_dir, libfastjetplugins_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libfastjetplugins_handle = dlopen(libfastjetplugins_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libfastjetplugins_path))

    global libfastjettools_path = normpath(joinpath(artifact_dir, libfastjettools_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libfastjettools_handle = dlopen(libfastjettools_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libfastjettools_path))

    global libsiscone_path = normpath(joinpath(artifact_dir, libsiscone_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libsiscone_handle = dlopen(libsiscone_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libsiscone_path))

    global libsiscone_spherical_path = normpath(joinpath(artifact_dir, libsiscone_spherical_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libsiscone_spherical_handle = dlopen(libsiscone_spherical_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libsiscone_spherical_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
using UnPack, OrderedCollections, ForwardDiff, JLD2, NLsolve, SparseArrays

mutable struct NKCapital{T <: Real}
    β::T
    γ::T
    φ::T
    ν::T
    χ::T
    δ::T
    α::T
    ϵ::T
    θ::T
    π_ss::T
    ϕ_r::T
    ϕ_π::T
    ϕ_y::T
    ρ_β::T
    ρ_l::T
    ρ_a::T
    ρ_r::T
    σ_β::T
    σ_l::T
    σ_a::T
    σ_r::T
    N_approx::Int
    S::OrderedDict{Symbol, Int}
    J::OrderedDict{Symbol, Int}
    E::OrderedDict{Symbol, Int}
    SH::OrderedDict{Symbol, Int}
end

function NKCapital(; β::T = .99, γ::T = 3.8, φ::T = 1., ν::T = 1., χ::T = 4.,
                   δ::T = 0.025, α::T = 0.33, ϵ::T = 10., θ::T = 0.7,
                   π_ss::T = 0., ϕ_r::T = 0.5,
                   ϕ_π::T = 1.3, ϕ_y::T = 0.25, ρ_β::T = 0.1,
                   ρ_l::T = 0.1, ρ_a::T = 0.9, ρ_r::T = 0.,
                   σ_β::T = 0.01, σ_l::T = 0.01, σ_a::T = 0.01, σ_r::T = 0.01,
                   N_approx::Int = 1) where {T <: Real}

    @assert N_approx > 0 "N_approx must be at least 1."

    ## Create Indexing dictionaries.

    # Note that for the exogenous shock
    # state variables, instead of e.g. η_L and η_A, I use η_l and η_a
    # since the uppercase variable will not appear in the jumps/states.
    S_init  = [:k₋₁, :v₋₁, :r₋₁, :output₋₁, :η_β, :η_l, :η_a, :η_r] # State Variables
    J_init  = [:output, :c, :l, :w, :r, :π, :q, :x, :rk, :mc,
               :s₁, :s₂, :v, :ω] # Jump variables
    E_init  = [:wage, :euler, :tobin, :cap_ret,
               :eq_mc, :kl_ratio, :eq_s₁, :eq_s₂,
               :phillips_curve, :price_dispersion,
               :mp, :output_market_clear, :production] # Equations
    SH_init = [:ε_β, :ε_l, :ε_a, :ε_r] # Exogenous shocks

    # Add approximations for forward-difference equations
    push!(E_init, :eq_omega)
    for var in [:q, :s₁, :s₂]
        inds = (var == :q) ? (1:N_approx) : (0:(N_approx - 1))
        push!(J_init, [Symbol(:d, var, "$(i)") for i in inds]...)
        push!(J_init, [Symbol(:p, var, "$(i)") for i in 1:N_approx]...)
        push!(E_init, [Symbol(:eq_d, var, "$(i)") for i in inds]...)
        push!(E_init, [Symbol(:eq_p, var, "$(i)") for i in 1:N_approx]...)
    end

    S  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(S_init))
    J  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(J_init))
    E  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(E_init))
    SH = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(SH_init))

    return NKCapital{T}(β, γ, φ, ν, χ, δ, α, ϵ, θ, π_ss, ϕ_r, ϕ_π, ϕ_y,
                        ρ_β, ρ_l, ρ_a, ρ_r, σ_β, σ_l, σ_a, σ_r,
                        N_approx, S, J, E, SH)
end

function nk_capital(m::NKCapital{T}; sparse_arrays::Bool = false,
                    sparse_jacobian::Vector{Symbol} = Symbol[]) where {T <: Real}

    # Get parameters
    @unpack β, γ, φ, ν, χ, δ, α, ϵ, θ, π_ss, ϕ_r, ϕ_π, ϕ_y = m
    @unpack ρ_β, ρ_l, ρ_a, ρ_r, σ_β, σ_l, σ_a, σ_r = m
    r_ss = π_ss - log(β)
    X̄    = δ * χ / (χ + 1.)

    # Unpack indexing dictionaries
    @unpack N_approx, S, J, E, SH = m
    @unpack k₋₁, v₋₁, r₋₁, output₋₁, η_β, η_l, η_a, η_r = S
    @unpack output, c, l, w, r, π, q, x, rk, ω, mc, s₁, s₂, v = J
    @unpack wage, euler, tobin, cap_ret, eq_mc, kl_ratio, eq_s₁, eq_s₂ = E
    @unpack phillips_curve, price_dispersion, mp = E
    @unpack output_market_clear, production, eq_omega = E
    @unpack ε_β, ε_l, ε_a, ε_r = SH

    Nz = length(S)
    Ny = length(J)
    Nε = length(SH)

    ## Define nonlinear equations

    # Some helper functions
    _Φ(Xin, Kin)  = X̄ ^ (1. / χ) / (1. - 1. / χ) * (Xin / Kin) ^ (1. - 1. / χ) - X̄ / (χ * (χ - 1.))
    _Φ′(Xin, Kin) = X̄ ^ (1. / χ) * (Xin / Kin) ^ (- 1. / χ)
    Φ(z, y)  = _Φ(exp(y[x]), exp(z[k₋₁]))
    Φ′(z, y) = _Φ′(exp(y[x]), exp(z[k₋₁]))
    m_ξ(z, y) = log(β) - z[η_β] + γ * y[c]
    function m_fwd!(i, Γ₅, Γ₆)
        Γ₅[i, η_β] = 1.
        Γ₆[i, c]   = -γ
    end
    pstar(y) = log(ϵ / (ϵ - 1.)) + y[s₁] - y[s₂]

    function μ(F, z, y)
        F[k₋₁]      = log(1 + X̄ ^ (1. / χ) / (1. - 1. / χ) *
            (exp(y[x] - z[k₋₁])) ^ (1. - 1. / χ) -
            X̄ / (1. - 1. / χ)) + z[k₋₁]
        F[v₋₁]      = y[v]
        F[r₋₁]      = y[r]
        F[output₋₁] = y[output]
        F[η_β]      = ρ_β * z[η_β]
        F[η_l]      = ρ_l * z[η_l]
        F[η_a]      = ρ_a * z[η_a]
        F[η_r]      = ρ_r * z[η_r]
    end

    function ξ(F, z, y)
        F_type = eltype(F)

        ## Pre-evaluate (just once) some terms
        Φv     = Φ(z, y)
        Φ′v    = Φ′(z, y)
        pstarv = pstar(y)
        m_ξv   = m_ξ(z, y)

        ## Non-forward-difference equations
        F[wage]                = log(φ) + z[η_l] + ν * y[l] - (-γ * y[c] + y[w])
        F[euler]               = y[r] + m_ξv
        F[tobin]               = y[q] + log(Φ′v)
        F[eq_mc]               = (1. - α) * y[w] + α * y[rk] - z[η_a] -
            (1. - α) * log(1. - α) - α * log(α) - y[mc]
        F[kl_ratio]            = z[k₋₁] - y[l] - log(α / (1. - α)) - (y[w] - y[rk])
        F[phillips_curve]      = (1. - ϵ) * y[π] - log((1. - θ) * exp((1. - ϵ) * (pstarv + y[π])) + θ)
        F[price_dispersion]    = y[v] - ϵ * y[π] - log((1. - θ) * exp(-ϵ * (pstarv + y[π])) + θ * exp(z[v₋₁]))
        F[mp]                  = (1. - ϕ_r) * r_ss + ϕ_r * z[r₋₁] + (1. - ϕ_r) .*
            (ϕ_π * (y[π] - π_ss) + ϕ_y * (y[output] - z[output₋₁])) + z[η_r] - y[r]
        F[output_market_clear] = y[output] - log(exp(y[c]) + exp(y[x]))
        F[production]          = z[η_a] + α * z[k₋₁] + (1. - α) * y[l] - y[v] - y[output]

        ## Forward-difference equations separately handled b/c recursions
        F[eq_omega] = log(1. - δ + Φv - Φ′v * exp(y[x] - z[k₋₁])) - y[ω]
        F[cap_ret]  = y[q] - log(sum([exp(y[J[Symbol("dq$(i)")]]) for i in 1:N_approx]) +
                                exp(y[J[Symbol("pq$(N_approx)")]]))
        F[eq_s₁]    = y[s₁] - log(sum([exp(y[J[Symbol("ds₁$(i)")]]) for i in 0:(N_approx - 1)]) +
                               exp(y[J[Symbol("ps₁$(N_approx)")]]))
        F[eq_s₂]    = y[s₂] - log(sum([exp(y[J[Symbol("ds₂$(i)")]]) for i in 0:(N_approx - 1)]) +
                               exp(y[J[Symbol("ps₂$(N_approx)")]]))

        # Set initial boundary conditions
        F[E[:eq_dq1]]  = -y[J[:dq1]] + m_ξv
        F[E[:eq_pq1]]  = -y[J[:pq1]] + m_ξv
        F[E[:eq_ds₁0]] = y[J[:ds₁0]] - y[mc] - y[output]
        F[E[:eq_ps₁1]] = log(θ) - y[J[:ps₁1]] + m_ξv
        F[E[:eq_ds₂0]] = y[J[:ds₂0]] - y[output]
        F[E[:eq_ps₂1]] = log(θ) - y[J[:ps₂1]] + m_ξv

        # Recursions for forward-difference equations
        for i in 2:N_approx
            F[E[Symbol("eq_dq$(i)")]]    = -y[J[Symbol("dq$(i)")]] + m_ξv
            F[E[Symbol("eq_pq$(i)")]]    = -y[J[Symbol("pq$(i)")]] + m_ξv
            F[E[Symbol("eq_ds₁$(i-1)")]] = log(θ) - y[J[Symbol("ds₁$(i-1)")]] + m_ξv
            F[E[Symbol("eq_ps₁$(i)")]]   = log(θ) - y[J[Symbol("ps₁$(i)")]]   + m_ξv
            F[E[Symbol("eq_ds₂$(i-1)")]] = log(θ) - y[J[Symbol("ds₂$(i-1)")]] + m_ξv
            F[E[Symbol("eq_ps₂$(i)")]]   = log(θ) - y[J[Symbol("ps₂$(i)")]]   + m_ξv
        end
    end

    # The cache is initialized as zeros so we only need to fill non-zero elements
    Λ = zeros(T, Nz, Ny)

    # The cache is initialized as zeros so we only need to fill non-zero elements
    function Σ(F, z)
        F[η_β, ε_β] = σ_β
        F[η_l, ε_l] = σ_l
        F[η_a, ε_a] = σ_a
        F[η_r, ε_r] = σ_r
    end

    function ccgf(F, α, z)
        # F .= .5 * diag(α * α') # slower but this is the underlying math
        sum!(F, α.^2) # faster implementation
        F .*= .5
    end

    ## Forward-looking variables
    Γ₅ = zeros(T, Ny, Nz)
    Γ₆ = zeros(T, Ny, Ny)

    # Equations w/out SDF terms and are not forward-difference equations
    Γ₆[euler, π] = -one(T)

    # Equations with SDF terms but are not forward-difference equations
    m_fwd!(euler, Γ₅, Γ₆)

    # Forward difference equations: boundary conditions
    m_fwd!(E[:eq_dq1], Γ₅, Γ₆)
    Γ₆[E[:eq_dq1], rk] = one(T)

    m_fwd!(E[:eq_pq1], Γ₅, Γ₆)
    Γ₆[E[:eq_pq1], q] = one(T)
    Γ₆[E[:eq_pq1], ω] = one(T)

    m_fwd!(E[:eq_ps₁1], Γ₅, Γ₆)
    Γ₆[E[:eq_ps₁1], π]  = convert(T, ϵ)
    Γ₆[E[:eq_ps₁1], s₁] = one(T)

    m_fwd!(E[:eq_ps₂1], Γ₅, Γ₆)
    Γ₆[E[:eq_ps₂1], π]  = convert(T, ϵ) - one(T)
    Γ₆[E[:eq_ps₂1], s₂] = one(T)

    # Forward difference equations: recursions
    for i in 2:N_approx
        m_fwd!(E[Symbol("eq_dq$(i)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_dq$(i)")], ω] = one(T)
        Γ₆[E[Symbol("eq_dq$(i)")], J[Symbol("dq$(i-1)")]] = one(T)

        m_fwd!(E[Symbol("eq_pq$(i)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_pq$(i)")], ω] = one(T)
        Γ₆[E[Symbol("eq_pq$(i)")], J[Symbol("pq$(i-1)")]] = one(T)

        m_fwd!(E[Symbol("eq_ds₁$(i-1)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_ds₁$(i-1)")], π] = convert(T, ϵ)
        Γ₆[E[Symbol("eq_ds₁$(i-1)")], J[Symbol("ds₁$(i-2)")]] = one(T)

        m_fwd!(E[Symbol("eq_ps₁$(i)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_ps₁$(i)")], π] = convert(T, ϵ)
        Γ₆[E[Symbol("eq_ps₁$(i)")], J[Symbol("ps₁$(i-1)")]] = one(T)

        m_fwd!(E[Symbol("eq_ds₂$(i-1)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_ds₂$(i-1)")], π] = convert(T, ϵ) - one(T)
        Γ₆[E[Symbol("eq_ds₂$(i-1)")], J[Symbol("ds₂$(i-2)")]] = one(T)

        m_fwd!(E[Symbol("eq_ps₂$(i)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_ps₂$(i)")], π] = convert(T, ϵ) - one(T)
        Γ₆[E[Symbol("eq_ps₂$(i)")], J[Symbol("ps₂$(i-1)")]] = one(T)
    end

    if sparse_arrays
        Γ₅ = sparse(Γ₅)
        Γ₆ = sparse(Γ₆)
    end

    ## Mapping from states to jump variables
    Ψ = zeros(T, Ny, Nz)

    ## Deterministic steady state as initial guess
    z = Vector{T}(undef, Nz)
    y = Vector{T}(undef, Ny)

    # AR(1) start at 0
    η_β0 = 0.
    η_l0 = 0.
    η_a0 = 0.
    η_r0 = 0.

    # Variables known outright
    M0 = β
    Q0 = 1.
    RK0 = 1. / β + X̄ - 1.

    # Guesses
    L0 = .5548429
    V0 = 1. # true if π_ss = 0, otherwise this is only a reasonable guess

    # Implied values given guesses
    C0_fnct = Cin -> Cin[1] + X̄ * (α / (1. - α) * φ * L0 ^ ν / Cin[1] ^ (-γ) / RK0 * L0) -
        (α / (1. - α) * φ * L0 ^ ν / Cin[1] ^ (-γ) / RK0) ^ α * L0 / V0
    C0_guess = NaN
    for theguess in .5:.5:10.
        try
            C0_fnct([theguess])
            C0_guess = theguess
        catch e
        end
    end
    C0 = nlsolve(C0_fnct, [C0_guess]).zero[1]
    W0 = φ * L0 ^ ν / C0 ^ (-γ)
    MC0 = (1. / (1. - α)) ^ (1. - α) * (1. / α) ^ α * W0 ^ (1. - α) * RK0 ^ α
    K0  = α / (1. - α) * W0 / RK0 * L0
    X0  = X̄ * K0
    Y0  = K0 ^ α * L0 ^ (1. - α) / V0
    S₁0  = MC0 * Y0 / (1. - θ * exp(π_ss) ^ ϵ)
    S₂0  = Y0 / (1. - θ * exp(π_ss) ^ (ϵ - 1.))
    Π0  = exp(π_ss)
    R0  = exp(r_ss)
    Ω0  = 1. - δ + _Φ(X0, K0) - _Φ′(X0, K0) * X0 / K0
    z .= [convert(T, x) for x in log.([K0, V0, R0, Y0, exp.([η_β0, η_l0, η_a0, η_r0])...])]
    y[1:14] = [convert(T, x) for x in log.([Y0, C0, L0, W0, R0, Π0, Q0, X0, RK0, Ω0, MC0, S₁0, S₂0, V0])]

    y[J[:dq1]] = convert(T, log(M0 * RK0))
    y[J[:pq1]] = convert(T, log(Ω0 * M0 * Q0))
    y[J[:ds₁0]] = convert(T, log(MC0 * Y0))
    y[J[:ps₁1]] = convert(T, log(θ * M0 * Π0^ϵ * S₁0))
    y[J[:ds₂0]] = convert(T, log(Y0))
    y[J[:ps₂1]] = convert(T, log(θ * M0 * Π0^(ϵ - 1.) * S₂0))

    for i in 2:N_approx
        y[J[Symbol("dq$(i)")]] = convert(T, log(M0) + log(Ω0) + y[J[Symbol("dq$(i-1)")]])
        y[J[Symbol("pq$(i)")]] = convert(T, log(M0) + log(Ω0) + y[J[Symbol("pq$(i-1)")]])
        y[J[Symbol("ds₁$(i-1)")]] = convert(T, log(θ) + log(M0) + ϵ * π_ss + y[J[Symbol("ds₁$(i-2)")]])
        y[J[Symbol("ps₁$(i)")]] = convert(T, log(θ) + log(M0) + ϵ * π_ss + y[J[Symbol("ps₁$(i-1)")]])
        y[J[Symbol("ds₂$(i-1)")]] = convert(T, log(θ) + log(M0) + (ϵ - 1.) * π_ss + y[J[Symbol("ds₂$(i-2)")]])
        y[J[Symbol("ps₂$(i)")]] = convert(T, log(θ) + log(M0) + (ϵ - 1.) * π_ss + y[J[Symbol("ps₂$(i-1)")]])
    end

    if sparse_arrays
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, vec(z), vec(y), Ψ, Nε; sparse_jacobian = sparse_jacobian,
                                         Σ_cache_init = dims -> spzeros(dims...),
                                         Λ_cache_init = dims -> spzeros(dims...))
    else
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, vec(z), vec(y), Ψ, Nε; sparse_jacobian = sparse_jacobian)
    end
end

nk_cₜ(m, zₜ) = exp(m.y[2] + (m.Ψ * (zₜ - m.z))[2])
nk_qₜ(m, zₜ) = exp(m.y[7] + (m.Ψ * (zₜ - m.z))[7])
nk_dqₜ(m, zₜ, i, J) = exp(m.y[J[Symbol("dq$(i)")]] + (m.Ψ * (zₜ - m.z))[J[Symbol("dq$(i)")]])
nk_pqₜ(m, zₜ, i, J) = exp(m.y[J[Symbol("pq$(i)")]] + (m.Ψ * (zₜ - m.z))[J[Symbol("pq$(i)")]])
nk_s₁ₜ(m, zₜ) = exp(m.y[12] + (m.Ψ * (zₜ - m.z))[12])
nk_s₂ₜ(m, zₜ) = exp(m.y[13] + (m.Ψ * (zₜ - m.z))[13])
nk_ds₁ₜ(m, zₜ, i, J) = exp(m.y[J[Symbol("ds₁$(i)")]] + (m.Ψ * (zₜ - m.z))[J[Symbol("ds₁$(i)")]])
nk_ps₁ₜ(m, zₜ, i, J) = exp(m.y[J[Symbol("ps₁$(i)")]] + (m.Ψ * (zₜ - m.z))[J[Symbol("ps₁$(i)")]])
nk_ds₂ₜ(m, zₜ, i, J) = exp(m.y[J[Symbol("ds₂$(i)")]] + (m.Ψ * (zₜ - m.z))[J[Symbol("ds₂$(i)")]])
nk_ps₂ₜ(m, zₜ, i, J) = exp(m.y[J[Symbol("ps₂$(i)")]] + (m.Ψ * (zₜ - m.z))[J[Symbol("ps₂$(i)")]])

# Evaluates Euler equation errors in log terms
function nk_log_euler(m, zₜ, εₜ₊₁, Cₜ; β::T = .99, γ::T = 3.8,
                      J::AbstractDict = NKCapital().J, S::AbstractDict = NKCapital().S) where {T <: Real}
    yₜ = m.y + m.Ψ * (zₜ - m.z)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)
    return log(β) - γ * (yₜ₊₁[J[:c]] - log(Cₜ)) +
        zₜ₊₁[S[:η_β]] - zₜ[S[:η_β]] + yₜ[J[:r]] - yₜ₊₁[J[:π]]
end

function nk_log_dq(m, zₜ, εₜ₊₁, DQₜ; β::T = .99, γ::T = 3.8,
                   i::Int = 1, J::AbstractDict = NKCapital().J, S::AbstractDict = NKCapital().S) where {T <: Real}
    yₜ = m.y + m.Ψ * (zₜ - m.z)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)
    mₜ₊₁ = log(β) - γ * (yₜ₊₁[J[:c]] - yₜ[J[:c]]) +
        zₜ₊₁[S[:η_β]] - zₜ[S[:η_β]]
    if i == 1
        return mₜ₊₁ + yₜ₊₁[J[:rk]] - log(DQₜ)
    else
        return yₜ₊₁[J[:ω]] + mₜ₊₁ + yₜ₊₁[J[Symbol("dq$(i-1)")]] - log(DQₜ)
    end
end

function nk_log_pq(m, zₜ, εₜ₊₁, PQₜ; β::T = .99, γ::T = 3.8,
                   i::Int = 1, J::AbstractDict = NKCapital().J, S::AbstractDict = NKCapital().S) where {T <: Real}
    yₜ = m.y + m.Ψ * (zₜ - m.z)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)
    mₜ₊₁ = log(β) - γ * (yₜ₊₁[J[:c]] - yₜ[J[:c]]) +
        zₜ₊₁[S[:η_β]] - zₜ[S[:η_β]]
    if i == 1
        return yₜ₊₁[J[:ω]] + mₜ₊₁ + yₜ₊₁[J[:q]] - log(PQₜ)
    else
        return yₜ₊₁[J[:ω]] + mₜ₊₁ + yₜ₊₁[J[Symbol("pq$(i-1)")]] - log(PQₜ)
    end
end

function nk_log_ds₁(m, zₜ, εₜ₊₁, DS₁ₜ; β::T = .99, γ::T = 3.8,
                    θ::T = 0.7, ϵ::T = 10., i::Int = 0,
                    J::AbstractDict = NKCapital().J, S::AbstractDict = NKCapital().S) where {T <: Real}
    yₜ = m.y + m.Ψ * (zₜ - m.z)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)
    mₜ₊₁ = log(β) - γ * (yₜ₊₁[J[:c]] - yₜ[J[:c]]) +
        zₜ₊₁[S[:η_β]] - zₜ[S[:η_β]]
    if i == 0
        return yₜ[J[:mc]] + yₜ[J[:output]]
    else
        return log(θ) + mₜ₊₁ + ϵ * yₜ₊₁[J[:π]] + yₜ₊₁[J[Symbol("ds₁$(i-1)")]] - log(DS₁ₜ)
    end
end

function nk_log_ps₁(m, zₜ, εₜ₊₁, PS₁ₜ; β::T = .99, γ::T = 3.8,
                    θ::T = 0.7, ϵ::T = 10., i::Int = 0,
                    J::AbstractDict = NKCapital().J, S::AbstractDict = NKCapital().S) where {T <: Real}
    yₜ = m.y + m.Ψ * (zₜ - m.z)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)
    mₜ₊₁ = log(β) - γ * (yₜ₊₁[J[:c]] - yₜ[J[:c]]) +
        zₜ₊₁[S[:η_β]] - zₜ[S[:η_β]]
    if i == 1
        return log(θ) + mₜ₊₁ + ϵ * yₜ₊₁[J[:π]] + yₜ₊₁[J[:s₂]] - log(PS₁ₜ)
    else
        return log(θ) + mₜ₊₁ + ϵ * yₜ₊₁[J[:π]] + yₜ₊₁[J[Symbol("ps₁$(i-1)")]] - log(PS₁ₜ)
    end
end

function nk_log_ds₂(m, zₜ, εₜ₊₁, DS₂ₜ; β::T = .99, γ::T = 3.8,
                    θ::T = 0.7, ϵ::T = 10., i::Int = 0,
                    J::AbstractDict = NKCapital().J, S::AbstractDict = NKCapital().S) where {T <: Real}
    yₜ = m.y + m.Ψ * (zₜ - m.z)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)
    mₜ₊₁ = log(β) - γ * (yₜ₊₁[J[:c]] - yₜ[J[:c]]) +
        zₜ₊₁[S[:η_β]] - zₜ[S[:η_β]]
    if i == 0
        return yₜ[J[:output]]
    else
        return log(θ) + mₜ₊₁ + (ϵ - 1.) * yₜ₊₁[J[:π]] + yₜ₊₁[J[Symbol("ds₂$(i-1)")]] - log(DS₂ₜ)
    end
end

function nk_log_ps₂(m, zₜ, εₜ₊₁, PS₂ₜ; β::T = .99, γ::T = 3.8,
                    θ::T = 0.7, ϵ::T = 10., i::Int = 0,
                    J::AbstractDict = NKCapital().J, S::AbstractDict = NKCapital().S) where {T <: Real}
    yₜ = m.y + m.Ψ * (zₜ - m.z)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)
    mₜ₊₁ = log(β) - γ * (yₜ₊₁[J[:c]] - yₜ[J[:c]]) +
        zₜ₊₁[S[:η_β]] - zₜ[S[:η_β]]
    if i == 1
        return log(θ) + mₜ₊₁ + (ϵ - 1.) * yₜ₊₁[J[:π]] + yₜ₊₁[J[:s₂]] - log(PS₂ₜ)
    else
        return log(θ) + mₜ₊₁ + (ϵ - 1.) * yₜ₊₁[J[:π]] + yₜ₊₁[J[Symbol("ps₂$(i-1)")]] - log(PS₂ₜ)
    end
end

# Evaluates n-period ahead Euler equation errors in log terms

# Calculate Euler equation via quadrature
std_norm_mean = zeros(4)
std_norm_sig  = ones(4)
if !isdefined(Main, :n_GH)
    n_GH = 5
end
nk_𝔼_quadrature(f::Function) = gausshermite_expectation(f, std_norm_mean, std_norm_sig, n_GH)

ext: jl


text:
# #define Abs(x)    ((x) < 0 ? -(x) : (x))
# #define Max(a, b) ((a) > (b) ? (a) : (b))

# double RelDif(double a, double b)
# {
# 	double c = Abs(a);
# 	double d = Abs(b);

# 	d = Max(c, d);

# 	return d == 0.0 ? 0.0 : Abs(a - b) / d;
# }
# Typical usage is
# 	if(RelDif(a, b) <= TOLERANCE) ...

using Base.Test

P = rand(Float64, 1000, 1000, 3);
P = P ./ sum(P, 1);

function test1{T<:FloatingPoint}(A::AbstractArray{T})
  r, c = size(A)[1:2]
  return (r == c) && all(abs(sum(A, 1) .- one(T)) .< 10eps(T))
end
@test test1(P)

function test2{T<:FloatingPoint}(A::AbstractArray{T})
  r, c = size(A)[1:2]
  ϵ = convert(T, 10)*eps(T)
  return (r == c) && all([abs(x - one(T)) < ϵ for x = sum(A, 1)])
end
@test test2(P)

function test2b{T<:FloatingPoint}(A::AbstractArray{T})
  r, c = size(A)[1:2]
  X = sum(A, 1)
  return (r == c) && all([abs(x - one(T)) < 5eps(max(x, one(T))) for x = X])
end
@test test2b(P)

function test3{T<:FloatingPoint}(A::AbstractArray{T})
  r, c = size(A)[1:2]
  ϵ = convert(T, 10)*eps(T)
  return (r == c) && all(Bool[abs(x - one(T)) < ϵ for x = sum(A, 1)])
end
@test test3(P)

function test4{T<:FloatingPoint}(A::AbstractArray{T})
  r, c = size(A)[1:2]
  ϵ = convert(T, 10)*eps(T)
  for x = sum(A, 1)
    abs(x - one(T)) < ϵ || return false
  end
  true
end
@test test4(P)

function test5{T<:FloatingPoint}(A::AbstractArray{T})
  r, c, d = size(A)
  r == c || return false
  ϵ = convert(T, 10)*eps(T)
  for k = 1:d
    for j = 1:c
      S = zero(T)
      @simd for i = 1:r
        @inbounds S += A[i, j, k]
      end
      abs(S - one(T)) < ϵ || return false
    end
  end
  true
end
@test test5(P)

function approx_equal_one{T}(x::T, max_rel_diff::T)
  oneT = one(T)
  diff = abs(x - oneT)
  x = abs(x)
  largest = (oneT > x) ? oneT : x
  diff <= largest * max_rel_diff || return false
  return true
end

function test6{T<:FloatingPoint}(A::AbstractArray{T})
  r, c, d = size(A)
  r == c || return false
  for k = 1:d
    for j = 1:c
      S = zero(T)
      @simd for i = 1:r
        @inbounds S += A[i, j, k]
      end
      approx_equal_one(S, 3eps(T)) || return false
    end
  end
  true
end
@test test6(P)

function test7{T<:FloatingPoint}(A::AbstractArray{T})
  r, c, d = size(A)
  r == c || return false
  for k = 1:d
    for j = 1:c
      S = zero(T)
      @simd for i = 1:r
        @inbounds S += A[i, j, k]
      end
      x = S < one(T) ? one(T) : S
      abs(S - one(T)) < 4eps(x) || return false
    end
  end
  true
end
@test test7(P)

# -------
# Timings
# -------

print("Test 1: ")
@time for i = 1:1000
  test1(P)
end

print("Test 2: ")
@time for i = 1:1000
  test2(P)
end

print("Test 2b: ")
@time for i = 1:1000
  test2b(P)
end

print("Test 3: ")
@time for i = 1:1000
  test3(P)
end

print("Test 4: ")
@time for i = 1:1000
  test4(P)
end

print("Test 5: ")
@time for i = 1:1000
  test5(P)
end

print("Test 6: ")
@time for i = 1:1000
  test6(P)
end

print("Test 7: ")
@time for i = 1:1000
  test6(P)
end


text:
type Reporter
    ctx::Context
    sock::Socket
    port::Int
    key::String
end

function Reporter(port::Int, key::String)
    ctx = Context()
    sock = Socket(ctx, REQ)
    connect(sock, "tcp://localhost:$port")
    return Reporter(ctx, sock, port, key)
end

Base.show(io::IO, r::Reporter) =
    print(io, "Reporter($(r.sock),$(r.port))")

@doc """Send reporting message""" ->
function report(r::Reporter, msg::String)
    send(r.sock, "r:$(r.key) $msg")
    resp = bytestring(recv(r.sock))
    if resp != "ok"
        error("Reporting failed with message: $resp")
    end
end


text:
module ReliabilityDiagrams

using Makie: Makie, @recipe
using ConsistencyResampling: ConsistencyResampling
using DataStructures: DataStructures
using RecipesBase: RecipesBase
using Statistics: Statistics
using StatsBase: StatsBase
using StructArrays: StructArrays

using Random: Random

export ConsistencyBars
export EqualSize, EqualMass

include("binning.jl")
include("consistency_bars.jl")
include("utils.jl")

include("recipes/makie.jl")
include("recipes/plots.jl")

include("conversions.jl")

end


text:
"""
Supertype of domains that are defined by an indicator function.

An indicator function is a function `f : S -> [0,1]` that indicates membership
of `x` to a domain `D` with `D ⊂ S`. The indicator function corresponds exactly
to the `in` function of a domain: `f(x) = x ∈ D`.

Concrete subtypes of `AbstractIndicatorFunction` store a representation of this
indicator function and implement `in` using that representation, rather than
implementing `in` directly.
"""
abstract type AbstractIndicatorFunction{T} <: Domain{T} end

"The indicator function of a domain is the function `f(x) = x ∈ D`."
indicatorfunction(d::Domain) = x -> x ∈ d

indomain(x, d::AbstractIndicatorFunction) = _indomain(x, d, indicatorfunction(d))
_indomain(x, d::AbstractIndicatorFunction, f) = f(x)

show(io::IO, d::AbstractIndicatorFunction) =
    print(io, "indicator domain defined by function f = $(indicatorfunction(d))")

"An `IndicatorFunction` is a domain that implements `f(x) = x ∈ D` by storing `f`."
struct IndicatorFunction{T,F} <: AbstractIndicatorFunction{T}
    f       ::  F
end

IndicatorFunction(f) = IndicatorFunction{Float64}(f)
IndicatorFunction{T}(f::F) where {T,F} = IndicatorFunction{T,F}(f)

indicatorfunction(d::IndicatorFunction) = d.f

similardomain(d::IndicatorFunction, ::Type{T}) where {T} = IndicatorFunction{T}(d.f)

convert(::Type{IndicatorFunction}, d::AbstractIndicatorFunction) = d
convert(::Type{IndicatorFunction}, d::Domain{T}) where {T} =
    IndicatorFunction{T}(indicatorfunction(d))

==(d1::IndicatorFunction, d2::IndicatorFunction) = indicatorfunction(d1)==indicatorfunction(d2)

intersectdomain1(d1::IndicatorFunction, d2) = BoundedIndicatorFunction(d1.f, d2)
intersectdomain2(d1, d2::IndicatorFunction) = BoundedIndicatorFunction(d2.f, d1)

"An indicator function with a known bounding domain."
struct BoundedIndicatorFunction{F,D,T} <: AbstractIndicatorFunction{T}
    f       ::  F
    domain  ::  D
end

BoundedIndicatorFunction(f::F, domain::D) where {F,T,D<:Domain{T}} =
    BoundedIndicatorFunction{F,D,T}(f, domain)

indicatorfunction(d::BoundedIndicatorFunction) = d.f

boundingdomain(d::BoundedIndicatorFunction) = d.domain

indomain(x, d::BoundedIndicatorFunction) = in(x, boundingdomain(d)) && d.f(x)

==(d1::BoundedIndicatorFunction, d2::BoundedIndicatorFunction) =
    indicatorfunction(d1)==indicatorfunction(d2) && boundingdomain(d1)==boundingdomain(d2)

similardomain(d::BoundedIndicatorFunction, ::Type{T}) where {T} =
    BoundedIndicatorFunction(d.f, convert(Domain{T}, d.domain))

Domain(gen::Base.Generator) = generator_domain(gen)

generator_domain(gen::Base.Generator{<:Domain}) = BoundedIndicatorFunction(gen.f, gen.iter)
generator_domain(gen::Base.Generator{<:Base.Iterators.ProductIterator}) =
    productgenerator_domain(gen, gen.iter.iterators)

function productgenerator_domain(gen, domains::Tuple{Vararg{Domain,N} where N})
    domain = TupleProductDomain(gen.iter.iterators)
    BoundedIndicatorFunction(gen.f, domain)
end

boundingbox(d::BoundedIndicatorFunction) = boundingbox(boundingdomain(d))

function show(io::IO, d::BoundedIndicatorFunction)
    print(io, "indicator function bounded by: ")
    show(io, boundingdomain(d))
end


text:
# proportion.jl
# =============
#
# Types and methods for computing p-distances.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE.md

"""
Proportion{T} is a distance which is the count of the mutations of type T that
exist between the two biological sequences, divided by the number of valid sites
examined (sites which don't have gap or ambiguous symbols).
In other words this so called p-distance is simply the proportion of sites
between each pair of sequences, that are mutated (again where T determines
what kind of mutation).
"""
struct Proportion{T} end

@inline function distance(::Type{Proportion{T}}, x::Sequence, y::Sequence) where T
    cs = count(T, x, y)
    return _pcorrection(Proportion{T}, cs, x, y)
end

@inline function pdistance(::Type{T}, x, y) where T
    return distance(Proportion{T}, x, y)
end

@inline function _pcorrection(k::Int, x::Sequence, y::Sequence)
    return k / min(length(x), length(y))
end

@inline function _pcorrection(cs::Tuple{Int,Int}, x::Sequence, y::Sequence)
    return cs[1] / cs[2]
end


text:
using FunctionalCollections
using FactCheck

facts("Persistent Queues") do

    context("equality") do
        @fact PersistentQueue(1:100) --> PersistentQueue(1:100)
        @fact PersistentQueue(1:100) --> not(PersistentQueue(1:99))
    end

    context("isempty") do
        @fact PersistentQueue{Int}() --> isempty
        @fact PersistentQueue([1])   --> not(isempty)
    end

    context("peek") do
        @fact peek(PersistentQueue(1:100)) --> 100
        @fact try peek(PersistentQueue{Int}()); false catch e true end --> true
    end

    context("pop") do
        @fact pop(PersistentQueue(1:100)) --> PersistentQueue(1:99)
        @fact try pop(PersistentQueue{Int}()); false catch e true end --> true
    end

    context("enq") do
        q = PersistentQueue{Int}()
        @fact peek(enq(q, 1)) --> 1
        @fact peek(pop(pop(enq(enq(enq(q, 1), 2), 3)))) --> 3
    end

    context("iteration") do
        arr2 = Int[]
        for i in PersistentQueue(1:1000)
            push!(arr2, i)
        end
        @fact arr2 --> collect(1000:-1:1)
    end

    context("hash") do
        @fact hash(PersistentQueue(1:1000)) --> hash(PersistentQueue(1:1000))
    end

    context("length") do
        @fact length(PersistentQueue([1, 2, 3])) --> 3
        @fact length(PersistentQueue{Int}()) --> 0

        @fact length(pop(PersistentQueue([1, 2, 3]))) --> 2
    end

end


text:
include(joinpath(@__DIR__, "testcommon.jl"))

@testset "$(basename(@__FILE__)[1:end-3])" begin

CxxWrap.CxxWrapCore.libexcept

@test ccall((:internalthrow, CxxWrap.CxxWrapCore.libexcept), Cint, (Cint,), -1) == 1
@test ccall((:internalthrow, CxxWrap.CxxWrapCore.libexcept), Cint, (Cint,), -2) == 2
@test ccall((:internalthrow, CxxWrap.CxxWrapCore.libexcept), Cint, (Cint,), -1) == 1

# This crashes when linking CxxWrap compiled with VS for 32-bit targets due to a mingw32 bug, see:
# https://ghc.haskell.org/trac/ghc/ticket/10435
# Only happens with SEHOP enabled, as is the case on appveyor
try
  ccall((:internalthrow, CxxWrap.CxxWrapCore.libexcept), Cint, (Cint,), 1) == 0
  @test false
catch
  println("exception 1")
  @test true
end

try
  ccall((:internalthrow, CxxWrap.CxxWrapCore.libexcept), Cint, (Cint,), 2) == 0
  @test false
catch
  println("exception 2")
  @test true
end

try
  ccall((:internalthrow, CxxWrap.CxxWrapCore.libexcept), Cint, (Cint,), 3) == 0
  @test false
catch
  println("exception 3")
  @test true
end

end


text:
using DrWatson
@quickactivate "INN_Velocity-Migration"

using JLD, PyPlot, LinearAlgebra, Statistics, JUDI.TimeModeling
include(srcdir("plot_rtm.jl"))

####################################################################################################
# This data is available in the SLIM dropbox folder at
# https://www.dropbox.com/sh/koxgvxuzg838s9x/AADivZzhLafitt4Q06XRWhYGa?dl=0
D = load(datadir("overthrust_images_train.jld"))
num = length(D["m"])

# Crop images and models to n2 x n1
n1 = 60
n2 = 200
m_all = zeros(Float32, n2, n1, num*4)
m0_all = zeros(Float32, n2, n1, num*4)

for i = 1:num

    m_all[:,:,4*i-3] = D["m"][i][1:n2, 6:n1+5] # keep 5 water layers out of 10 water layers
    m_all[:,:,4*i-2] = D["m"][i][n2+1:n2*2, 6:n1+5]
    m_all[:,:,4*i-1] = D["m"][i][1:n2, n1+1:n1*2]
    m_all[:,:,4*i] = D["m"][i][n2+1:n2*2, n1+1:n1*2]

    m0_all[:,:,4*i-3] = D["m0"][i][1:n2, 6:n1+5] # keep 5 water layers out of 10 water layers
    m0_all[:,:,4*i-2] = D["m0"][i][n2+1:n2*2, 6:n1+5]
    m0_all[:,:,4*i-1] = D["m0"][i][1:n2, n1+1:n1*2]
    m0_all[:,:,4*i] = D["m0"][i][n2+1:n2*2, n1+1:n1*2]

end

m_all[:,1:5,:] .= m_all[1,1,1]
m0_all[:,1:5,:] .= m0_all[1,1,1]

# add 4 water layers, i.e., totally 9 water layers
n2, n1, num = size(m_all)
m_temp = m_all[1,1,1] .* ones(n2, 4, num)
m_all = cat(m_temp, m_all; dims = 2)
m0_all = cat(m_temp, m0_all; dims = 2)

save(datadir("overthrust_8k_models_200x64.jld"), "m_all", m_all, "m0_all", m0_all)

# seperate the models w.r.t velocity
idx = sort([ collect(range(1,num, step=4)); collect(range(2,num, step=4))])
m_all1 = m_all[:,:,idx]
m_all2 = m_all[:,:,setdiff(1:end, idx)]
m0_all1 = m0_all[:,:,idx]
m0_all2 = m0_all[:,:,setdiff(1:end, idx)]

m_mean1 = mean(mean(m_all1, dims=3), dims=1)
m_mean2 = mean(mean(m_all2, dims=3), dims=1)
v_mean1 = sqrt.(1f0 ./ m_mean1) # velocity(km/s): 1.8, 3.5->4.3
v_mean2 = sqrt.(1f0 ./ m_mean2) # velocity(km/s): 1.8, 4.3->5.5

save(datadir("overthrust_4k_models_200x64.jld"), "m_all1", m_all1, "m0_all1", m0_all1, "m_all2", m_all2, "m0_all2", m0_all2)

####################################################################################################
m_all, m0_all = load(datadir("overthrust_4k_models_200x64.jld"), "m_all1", "m0_all1")

figfolder = "figs/rtm_64x200_lin_vel_lsrtmobj"
mkpath(datadir(figfolder))

## Computing RTM images
i = 1

m = m_all[:,:,i]
n = size(m)
# m0 = m0_all[:,:,i] # migration velocity
# m0 = ones(Float32, n) .* m[1,1] # constant velocity

# build a migration velocity m0 whose value increases linearly with depth
idx_wb = 9 # index of the water bottom
v0 = ones(1, n[2])
v0[idx_wb+1:end] = collect(range(3.5, 4.3, length = n[2]-idx_wb))
v0 = (1f0 ./ v0).^2
v0[1:idx_wb] .= m[1,1]
m0 = repeat(v0, n[1], 1)

dm = m - m0

figure(figsize=[8,8])
ax1 = subplot(3, 1, 1); imshow(adjoint(m), cmap="jet"); plt.colorbar(); title(L"True velocity ($s^2/km^2$)")
ax2 = subplot(3, 1, 2); imshow(adjoint(m0), cmap="jet"); plt.colorbar(); title(L"Migration velocity ($s^2/km^2$)")
ax3 = subplot(3, 1, 3); imshow(adjoint(dm), cmap="jet"); plt.colorbar(); title(L"Perturbation ($s^2/km^2$)")
savefig(datadir(figfolder, "vel_example.png"))

# Computing RTM images
n = size(m)
dm = vec(m - m0)

# Setup info and model structure
d = (25., 25.)
o = (0., 0.)
nsrc = 21	# number of sources

## Set up receiver geometry in meters
nxrec = 201
xrec = range(0f0, stop=5000f0, length=nxrec)
yrec = 0f0
zrec = range(0f0, stop=0f0, length=nxrec)

# receiver sampling and recording time
timeR = 2000f0   # receiver recording time [ms]
dtR = 1f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(range(0f0, stop=5000f0, length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(0f0, stop=0f0, length=nsrc))

# source sampling and number of time steps
timeS = timeR  # ms
dtS = dtR   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.006f0     # KHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet);

# Set up model and info structure for linear operators
model = Model(n, d, o, m; nb=200)
model0 = Model(n, d, o, m0; nb=200)
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

opt = Options(space_order=16, isic=true) # Use linearized inverse scattering imaging condition for the Jacobian

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr * F0 * adjoint(Ps), q) # Born forward modeling operator

# Forward modeling and RTM
d_obs = Pr * F * adjoint(Ps) * q  # Practical observation data generated on m
# d_syn = Pr*F0*adjoint(Ps)*q # Synthetic observation data generated on m0
# fval, rtm = fwi_objective(model0, q, d_obs; options=opt)

# or using lsrtm_objective() like
# x = zeros(Float32, model0.n)
# fval, rtm = lsrtm_objective(model0, q, d_obs, x; nlind=true, options=opt) # Mr is a preconditioner like model topmute

# or compute RTM directly ()
# d_syn = Pr * F0 * adjoint(Ps) * q # Synthetic observation data generated on m0
# rtm = adjoint(J) * (d_obs - d_syn)

# When the migarion velocity is close to the true velocity
m1 = ones(Float32, n) .* m[1,1] # constant velocity
m1[:, idx_wb+1:end] .= m[:, idx_wb+1]
model1 = Model(n, d, o, m1; nb=200)
F1 = judiModeling(info, model1; options=opt)
d_syn = Pr*F1*adjoint(Ps)*q # forward modeling on the constant velocity to mute the direct wave
rtm = adjoint(J) * (d_obs - d_syn)

# rtm = rtm_isic(m, m0)
rtm1 = adjoint(reshape(rtm, n))
figure(figsize=[8,2]); imshow(rtm1, cmap="gray", vmin=-2e-2, vmax=2e-2); plt.colorbar(); title("RTM_isic result")
savefig(datadir(figfolder, "rtm.png"))

plot_rtm(reshape(rtm, n), d)
savefig(datadir(figfolder, "rtm1.png"))

figure(figsize=[12,8])
ax1 = subplot(1,3,1); imshow(d_obs.data[1], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_obs from src1")
ax2 = subplot(1,3,2); imshow(d_obs.data[11], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_obs from src11")
ax3 = subplot(1,3,3); imshow(d_obs.data[21], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_obs from src21")
savefig(datadir(figfolder, "d_obs.png"))

# figure(figsize=[12,8])
# ax1 = subplot(1,3,1); imshow(d_syn.data[1], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_syn from src1")
# ax2 = subplot(1,3,2); imshow(d_syn.data[11], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_syn from src11")
# ax3 = subplot(1,3,3); imshow(d_syn.data[21], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_syn from src21")
# savefig(datadir(figfolder, "d_syn.png"))

# figure(figsize=[12,8])
# ax1 = subplot(1,3,1); imshow(d_obs.data[1]-d_syn.data[1], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_diff from src1")
# ax2 = subplot(1,3,2); imshow(d_obs.data[11]-d_syn.data[11], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_diff from src11")
# ax3 = subplot(1,3,3); imshow(d_obs.data[21]-d_syn.data[21], vmin=-1, vmax=1, cmap="gray", aspect="auto"); title("d_diff from src21")
# savefig(datadir(figfolder, "d_diff.png"))

# Right-hand preconditioners (model topmute)
# assign a mask where depth 0 to idx_wb-2 should be exactly 0,
# idx_wb-2 to idx_wb go from 0 to 1 smoothly,
# after idx_wb is exactly 1
d = (25., 25.)
o = (0., 0.)
model0 = Model(n, d, o, m0)
idx_wb = find_water_bottom(reshape(dm, n)) # find the index of the water bottom (you could also manually add this)
Tm = judiTopmute(model0.n, idx_wb, 2)  # Mute water column
S = judiDepthScaling(model0) # design a depth scaling operator
Mr = S * Tm
rtm = Mr' * rtm # mute water layers
rtm2 = adjoint(reshape(rtm, n))
figure(figsize=[8,2]); imshow(rtm2, cmap="gray", vmin=-2e-2, vmax=2e-2); plt.colorbar(); title("RTM_isic result with topmute")
savefig(datadir(figfolder, "rtm_TopMute.png"))

plot_rtm(reshape(rtm, n), d)
savefig(datadir(figfolder, "rtm2.png"))

figure(figsize=[8,2]); imshow(rtm1 - rtm2, cmap="gray", vmin=-2e-2, vmax=2e-2); plt.colorbar(); title("Difference of RTM_isic results")
savefig(datadir(figfolder, "rtm_difference.png"))


text:
type ParameterizedFunction{isinplace,F,P} <: ConstructedParameterizedFunction{isinplace}
  f::F
  params::P
end

Base.@pure function ParameterizedFunction(f,p)
  isinplace = numargs(f)>=4
  ParameterizedFunction{isinplace,typeof(f),typeof(p)}(f,p)
end

(pf::ParameterizedFunction{true,F,P}){F,P}(t,u,du) = pf.f(t,u,pf.params,du)
(pf::ParameterizedFunction{true,F,P}){F,P}(t,u,params,du) = pf.f(t,u,params,du)
(pf::ParameterizedFunction{false,F,P}){F,P}(t,u) = pf.f(t,u,pf.params)
(pf::ParameterizedFunction{false,F,P}){F,P}(t,u,params) = pf.f(t,u,params)

type DAEParameterizedFunction{isinplace,F,P} <: ConstructedParameterizedFunction{isinplace}
  f::F
  params::P
end

Base.@pure function DAEParameterizedFunction(f,p)
  isinplace = numargs(f)>=5
  DAEParameterizedFunction{isinplace,typeof(f),typeof(p)}(f,p)
end

(pf::DAEParameterizedFunction{true,F,P}){F,P}(t,u,du,out) = pf.f(t,u,pf.params,du,out)
(pf::DAEParameterizedFunction{true,F,P}){F,P}(t,u,params,du,out) = pf.f(t,u,params,du,out)
(pf::DAEParameterizedFunction{false,F,P}){F,P}(t,u,du) = pf.f(t,u,pf.params,du)
(pf::DAEParameterizedFunction{false,F,P}){F,P}(t,u,params,du) = pf.f(t,u,params,du)

type DDEParameterizedFunction{isinplace,F,P} <: ConstructedParameterizedFunction{isinplace}
  f::F
  params::P
end

Base.@pure function DDEParameterizedFunction(f,p)
  isinplace = numargs(f)>=5
  DDEParameterizedFunction{isinplace,typeof(f),typeof(p)}(f,p)
end

(pf::DDEParameterizedFunction{true,F,P}){F,P}(t,u,h,du) = pf.f(t,u,h,pf.params,du)
(pf::DDEParameterizedFunction{true,F,P}){F,P}(t,u,h,params,du) = pf.f(t,u,h,params,du)
(pf::DDEParameterizedFunction{false,F,P}){F,P}(t,u,h) = pf.f(t,u,h,pf.params)
(pf::DDEParameterizedFunction{false,F,P}){F,P}(t,u,h,params) = pf.f(t,u,h,params)

type ProbParameterizedFunction{F,P} <: ConstructedParameterizedFunction{false}
  f::F
  params::P
end
(pf::ProbParameterizedFunction{F,P}){F,P}(prob,i) = pf.f(prob,i,params)

type OutputParameterizedFunction{F,P} <: ConstructedParameterizedFunction{false}
  f::F
  params::P
end
(pf::OutputParameterizedFunction{F,P}){F,P}(sol) = pf.f(sol,params)

### Interface

DiffEqBase.param_values(pf::ConstructedParameterizedFunction) = pf.params
DiffEqBase.num_params(pf::ConstructedParameterizedFunction) = length(pf.params)


text:
function tst_comp_refs()
    ## Create plink files of the 3 groups of data
    smpl = "dat/v4r"
    wdir = "$smpl/tmp"

    @info "Merging v4 genotypes into plink format"
    #-- Note: I substitued original ID names with integers, as they are not important here.
    dir = "$smpl/18k_v4"
    lmp = "$smpl/map/V4_used/18k_v4.map"
    tgt = "$wdir/v4"
    make_bed(dir, lmp, tgt)

    @info "Merging v2 genotypes into plink format"
    dir = "$smpl/18k_v2"
    lmp = "$smpl/map/V4_used/18k_v2.map"
    tgt = "$wdir/v2"
    make_bed(dir, lmp, tgt)

    @info "Merging 8k genotypes into plink format"
    dir = "$smpl/8k"
    lmp = "$smpl/map/V4_used/8k_rudi.map"
    tgt = "$wdir/8k"
    make_bed(dir, lmp, tgt)

    ## QC of these 3 groups of plink files.
    #--- attention: if the smpl is too small, then many loci may be removed
    #---   change QC criteria, or use larger smpl
    QC("$wdir/v4", "$wdir/v4-qc")
    QC("$wdir/v2", "$wdir/v2-qc")
    QC("$wdir/8k", "$wdir/8k-qc")

    ## Make references
    ### v4
    make_ref("$wdir/v4-qc", "$wdir/v4-ref")
    ### v2
    make_ref("$wdir/v2-qc", "$wdir/v2-ref")

    ## Do the comparison and report results
    comp_refs("$wdir/v4-ref", "$wdir/v2-ref", "$wdir/8k-qc")
end


text:
module tdpe

include("common.jl")

#include individual problems
include("001.jl")
include("002.jl")
include("003.jl")
include("004.jl")
include("005.jl")
include("006.jl")
include("007.jl")
include("008.jl")
include("009.jl")
include("010.jl")
#end includes

function print_solutions()
    #print individual solutions
    solution001()
    solution002()
    solution003()
    solution004()
    solution005()
    solution006()
    solution007()
    solution008()
    solution009()
    solution010()
    #end prints
end

end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
# Structured optimizer interface
# ========================== #
"""
    supports_structure(optimizer::StochasticProgramOptimizerType, structure::AbstractStochasticStructure)

Return a `Bool` indicating whether `optimizer` supports the stochastic `structure`. That is, `load_structure!(optimizer, structure)` will not throw `UnsupportedStructure`
"""
function supports_structure(optimizer::StochasticProgramOptimizerType, structure::AbstractStochasticStructure)
    return false
end
"""
    check_loadable(optimizer::AbstractStructuredOptimizer, structure::AbstractStochasticStructure)

Throws an `UnloadableStructure` exception if `structure` is not loadable by `optimizer`.

    See also: [`load_structure!`](@ref)
"""
function check_loadable end
"""
    load_structure!(optimizer::AbstractStructuredOptimizer, structure::AbstractStochasticStructure, x₀::AbstractVector)

Instantiate the `optimizer` with the stochastic program represented in memory by the given `structure` and initial decision `x₀`.

See also: [`optimize!`](@ref)
"""
function load_structure! end
"""
    restore_structure!(optimizer::AbstractStructuredOptimizer)

Restore the stochastic program to the state it was in before a call to `optimize!`

See also: [`load_structure!`](@ref)
"""
function restore_structure! end
"""
    optimize!(optimizer::AbstractStructuredOptimizer)

Start the solution procedure for `optimizer` after a call to [`load_structure!`](@ref).

See also: [`load_structure!`](@ref)
"""
function optimize! end
"""
    optimizer_name(optimizer::AbstractStructuredOptimizer)

Optionally, return a string identifier of `AbstractStructuredOptimizer`.
"""
function optimizer_name(::AbstractStructuredOptimizer)
    return "SolverName() attribute not implemented by the optimizer."
end
"""
    master_optimizer(optimizer::AbstractStructuredOptimizer)

Return a MOI optimizer constructor
"""
function master_optimizer end
"""
   subproblem_optimizer(optimizer::AbstractStructuredOptimizer)

Return a MOI optimizer constructor for solving subproblems
"""
function subproblem_optimizer end
# Sample-based solver interface
# ========================== #
"""
    load_model!(optimizer::AbstractSampledOptimizer, model::StochasticModel, x₀::AbstractVector)

Instantiate the `optimizer` with the stochastic model and initial decision `x₀`.

See also: [`optimize!`](@ref)
"""
function load_model! end
"""
    optimizer_name(optimizer::AbstractSampledOptimizer)

Optionally, return a string identifier of `AbstractSampledOptimizer`.
"""
function optimizer_name(::AbstractSampledOptimizer)
    return "SolverName() attribute not implemented by the optimizer."
end
"""
    optimal_instance(optimizer::AbstractSampledOptimizer)

Return a stochastic programming instance of the stochastic model after a call to [`optimize!`](@ref).
"""
function optimal_instance end


text:
using Test

# Helper function for numerical tests
function check_numerical(
  chain,
  symbols::Vector{Symbol},
  exact_vals::Vector;
  eps=0.2,
)
	for (sym, val) in zip(symbols, exact_vals)
		E = mean(chain[sym])
		print("  $sym = $E ≈ $val (eps = $eps) ?")
		cmp = abs.(sum(mean(chain[sym]) - val)) <= eps
		if cmp
			printstyled("./\n", color = :green)
			printstyled("    $sym = $E, diff = $(abs.(E - val))\n", color = :green)
		else
		printstyled(" X\n", color = :red)
		printstyled("    $sym = $E, diff = $(abs.(E - val))\n", color = :red)
	end
  end
end

"""
  getteststoskip(filepath)

Returns a list of files (tests) to skip within a test group.
"""
function getteststoskip(filepath)
  if isfile(filepath)
	lines = readlines(filepath)
	lines = filter(line -> endswith(line, ".jl"), lines)
	lines = filter(line -> !startswith(line, "#"), lines)
	lines = filter(line -> length(split(line)) == 1, lines)
	lines = map(line -> strip(line), lines)
	return Set{String}(lines)
  else
	return Set{String}()
  end
end

function insdelim(c, deli=",")
	return reduce((e, res) -> append!(e, [res, deli]), c; init = [])[1:end-1]
end

"""
  runtests(; test_folders = ["all"], exclude = [])

Run specified Turing tests. By default this function runs all Turing tests.
The user can specify a list of test sets that she wants to run if necessary.
Specific files can be excluded by adding paths to the `exclude` argument.

Example:
```julia
runtests(test_folders = ["util.jl", "trace.jl"])
```

Excluding a specific file:
```julia
runtests(exclude = [joinpath("hmc.jl", "matrix_support.jl"])
```

Running a specific test:
```julia
runtests(specific_tests = [joinpath("hmc.jl", "matrix_support.jl"])
```

Note that calling `runtests` with a nonzero length `specific_tests` argument
overwrites all other function arguments.
"""
function runtests(; test_folders = ["all"], exclude = [], specific_tests = [])
	# test groups
	CORE_TESTS = ["ad.jl", "compiler.jl", "container.jl", "varinfo.jl",
		# "io.jl",
		"util.jl"]
	SAMPLER_TESTS = ["resample.jl", "adapt.jl", "vectorisation.jl", "gibbs.jl", "nuts.jl",
		"hmcda.jl", "hmc_core.jl", "hmc.jl", "sghmc.jl", "sgld.jl", "is.jl",
		"mh.jl",
		# "pmmh.jl", "ipmcmc.jl", "pgibbs.jl", "smc.jl"
		]
	TRACE_TESTS = ["trace.jl"]
	ALL = union(CORE_TESTS, SAMPLER_TESTS, TRACE_TESTS)

	# test groups that should be executed
	TEST_GROUPS = "all" ∈ test_folders ? ALL : test_folders

	# If we have specific tests, only run those.
	has_specific_tests = length(specific_tests) > 0
	if has_specific_tests
		TEST_GROUPS = []
		for i in specific_tests
			TEST_GROUPS = TEST_GROUPS ∪ [dirname(i)]
		end
	end

	# Run tests
	tr = @testset "Turing tests" begin
	for test_group in TEST_GROUPS
		teststoskip = getteststoskip(joinpath(test_group, "skip_tests"))
			@testset "$(test_group)" begin
				for test in filter(f -> endswith(f, ".jl")
					&& !(f ∈ teststoskip), readdir(test_group))

					@testset "$(test)" begin
						test_path = joinpath(test_group, test)

						# If we've got specific tests, examine those.
						if has_specific_tests
							if test_path in specific_tests
								include(joinpath(test_group, test))
							end
						# Otherwise look at every file in the folder.
						else
							if test_path in exclude
								@info "Skipping $test_path."
							else
								include(joinpath(test_group, test))
							end
						end

					end
				end
			end
		end
	end

	println()
	Test.print_test_results(tr)
end

using Pkg;
"""
	isinstalled(x::String)
Check if a package is installed.
"""
isinstalled(x::AbstractString) = x ∈ keys(Pkg.installed())

# NOTE: Remove the code below when DynamicHMC is registered.
using Pkg;
isinstalled("TransformVariables") || pkg"add https://github.com/tpapp/TransformVariables.jl#master";
isinstalled("LogDensityProblems") || pkg"add https://github.com/tpapp/LogDensityProblems.jl#master";
isinstalled("DynamicHMC") || pkg"add https://github.com/tpapp/DynamicHMC.jl#master";


text:
const _logScales = [:ln, :log2, :log10]
const _logScaleBases = Dict(:ln => ℯ, :log2 => 2.0, :log10 => 10.0)

# NOTE: This file was moved from Gadfly.jl, and the original author is Daniel Jones (@dcjones)

# Find the smallest order of magnitude that is larger than xspan This is a
# little opaque because I want to avoid assuming the log function is defined
# over typeof(xspan)
function bounding_order_of_magnitude(xspan::DT, base) where DT
    one_dt = convert(DT, one(DT))

    a = 1
    step = 1
    while xspan < base^a * one_dt
        a -= step
    end

    b = 1
    step = 1
    while xspan > base^b * one_dt
        b += step
    end

    while a + 1 < b
        c = div(a + b, 2)
        if xspan < base^c * one_dt
            b = c
        else
            a = c
        end
    end

    return b
end

const float_digit_range = floor(Int,log10(floatmin())):ceil(Int,log10(floatmax()))
postdecimal_digits(x) = first(i for i in float_digit_range if x==floor(x; digits=i))

# Empty catchall
optimize_ticks() = Any[]

"""
    optimize_ticks(xmin, xmax; extend_ticks::Bool=false,
                           Q=[(1.0,1.0), (5.0, 0.9), (2.0, 0.7), (2.5, 0.5), (3.0, 0.2)],
                           k_min::Int=2, k_max::Int=10, k_ideal::Int=5,
                           granularity_weight::Float64=1/4, simplicity_weight::Float64=1/6,
                           coverage_weight::Float64=1/3, niceness_weight::Float64=1/4,
                           strict_span=true, span_buffer = nothing
        )

Find some reasonable values for tick marks.

This is basically Wilkinson's ad-hoc scoring method that tries to balance
tight fit around the data, optimal number of ticks, and simple numbers.

## Arguments:

*    `xmax`:

    The maximum value occurring in the data.

*    `xmin`:

    The minimum value occurring in the data.

*    `extend_ticks`:

    Determines whether to extend tick computation. Defaults to false.

*    `strict_span`:

    True if no ticks should be outside [x_min, x_max]. Defaults to true.

*    `Q`:

    A distribution of nice numbers from which labellings are sampled. Stored in the form (number, score).

*    `k_min`:

    The minimum number of ticks.

*    `k_max`:

    The maximum number of ticks.

*    `k_ideal`:

    The ideal number of ticks.

*    `granularity_weight`:

    Encourages returning roughly the number of labels requested.

*    `simplicity_weight`:

    Encourages nicer labeling sequences by preferring step sizes that appear earlier in Q.
    Also rewards labelings that include 0 as a way to ground the sequence.

*    `coverage_weight`:

    Encourages labelings that do not extend far beyond the range of the data, penalizing unnecessary whitespace.

*    `niceness_weight`:

    Encourages labellings to produce nice ranges.

## Returns:
  `(ticklocations::Vector{Float64}, x_min, x_max)`

## Mathematical details

Wilkinson’s optimization function is defined as the sum of three
components. If the user requests m labels and a possible labeling has
k labels, then the components are `simplicity`, `coverage` and `granularity`.

These components are defined as follows:
```math
\\begin{aligned}
  &\\text{simplicity} = 1 - \\frac{i}{|Q|} + \\frac{v}{|Q|}\\\\
  &\\text{coverage}   = \\frac{x_{max} - x_{min}}{\\mathrm{label}_{max} - \\mathrm{label}_{min}}\\\\
  &\\text{granularity}= 1 - \\frac{\\left|k - m\\right|}{m}
\\end{aligned}
```

and the variables here are:

*  `q`: element of `Q`.
*  `i`: index of `q` in `Q`.
*  `v`: 1 if label range includes 0, 0 otherwise.
"""
function optimize_ticks(x_min::T, x_max::T; extend_ticks::Bool=false,
                           Q=[(1.0,1.0), (5.0, 0.9), (2.0, 0.7), (2.5, 0.5), (3.0, 0.2)],
                           k_min::Int=2, k_max::Int=10, k_ideal::Int=5,
                           granularity_weight::Float64=1/4, simplicity_weight::Float64=1/6,
                           coverage_weight::Float64=1/3, niceness_weight::Float64=1/4,
                           strict_span=true, span_buffer=nothing, scale=nothing) where T

    Qv = [(Float64(q[1]), Float64(q[2])) for q in Q]
    optimize_ticks_typed(x_min, x_max, extend_ticks, Qv, k_min, k_max, k_ideal,
                         granularity_weight, simplicity_weight,
                         coverage_weight, niceness_weight, strict_span, span_buffer, scale)
end

function optimize_ticks_typed(x_min::T, x_max::T, extend_ticks,
                           Q::Vector{Tuple{Float64,Float64}}, k_min, k_max, k_ideal,
                           granularity_weight::Float64, simplicity_weight::Float64,
                           coverage_weight::Float64, niceness_weight::Float64,
                           strict_span, span_buffer, scale) where T
    one_t = convert(T, one(T))
    if x_max - x_min < eps() * one_t
        R = typeof(1. * one_t)
        return R[x_min], x_min - one_t, x_min + one_t
    end

    n = length(Q)
    is_log_scale = scale ∈ _logScales
    base = get(_logScaleBases, scale, 10.)

    # generalizing "order of magnitude"
    xspan = x_max - x_min
    z = bounding_order_of_magnitude(xspan, base)

    # find required significant digits for ticks with q*base^z spacing,
    # for q values specified in Q
    x_digits = bounding_order_of_magnitude(max(abs(x_min), abs(x_max)), base)
    q_extra_digits = maximum(postdecimal_digits(q[1]) for q in Q)
    sigdigits(z) = max(1, x_digits - z + q_extra_digits)

    ib = Int(base)
    round_base = (
        isinteger(base) ?
        v -> round(v, sigdigits=sigdigits(z), base=ib) :
        v -> round(v, sigdigits=sigdigits(z))
    )

    high_score = -Inf
    S_best = Array{typeof(1. * one_t)}(undef, 1)
    viewmin_best, viewmax_best = x_min, x_max

    # we preallocate arrays that hold all required S arrays for every given
    # the k parameter, so we don't have to create them again and again, which
    # saves many allocations
    prealloc_Ss = if extend_ticks
        [Array{typeof(1. * one_t)}(undef, Int(3k)) for k in k_min:2k_max]
    else
        [Array{typeof(1. * one_t)}(undef, k) for k in k_min:2k_max]
    end

    while 2k_max * base^(z + 1) * one_t > xspan
        for (ik, k) in enumerate(k_min:2k_max)
            for (q, qscore) in Q
                tickspan = q * base^z * one_t
                span = (k - 1) * tickspan
                span < xspan && continue

                stp = q*base^z
                stp < eps() && continue

                if is_log_scale && !isinteger(stp)
                    continue  # prefer integer exponents for log scales
                end
                r = (x_max - span) / (stp * one_t)
                isfinite(r) || continue
                r = ceil(Int64, r)

                while r*stp * one_t <= x_min
                    # Filter or expand ticks
                    if extend_ticks
                        S = prealloc_Ss[ik]
                        for i in 0:(3k - 1)
                            S[i+1] = (r + i - k) * tickspan
                        end
                        # round only those values that end up as viewmin and viewmax
                        # to save computation time
                        S[k + 1] = round_base(S[k + 1])
                        S[2 * k] = round_base(S[2 * k])
                        viewmin, viewmax = S[k + 1], S[2 * k]
                    else
                        S = prealloc_Ss[ik]
                        for i in 0:(k - 1)
                            S[i+1] = (r + i) * tickspan
                        end
                        # round only those values that end up as viewmin and viewmax
                        # to save computation time
                        S[1] = round_base(S[1])
                        S[k] = round_base(S[k])
                        viewmin, viewmax = S[1], S[k]
                    end
                    if strict_span
                        viewmin = max(viewmin, x_min)
                        viewmax = min(viewmax, x_max)
                        buf = something(span_buffer, 0) * (viewmax - viewmin)

                        # filter the S array while reusing its own memory to do so
                        # this works because S is sorted, and we will only overwrite
                        # values that are not needed anymore going forward in the loop

                        # we do this because it saves allocations and leaves S type stable
                        counter = 0
                        @inbounds for i in 1:length(S)
                            if (viewmin - buf) <= S[i] <= (viewmax + buf)
                                counter += 1
                                S[counter] = S[i]
                            end
                        end
                        S = view(S, 1:counter)
                    end

                    # evaluate quality of ticks

                    has_zero = r <= 0 && abs(r) < k

                    # simplicity
                    s = has_zero ? 1 : 0

                    # granularity
                    g = 0 < length(S) < 2k_ideal ? 1 - abs(length(S) - k_ideal) / k_ideal : 0

                    # coverage
                    effective_span = (length(S) - 1) * tickspan
                    c = 1.5xspan / effective_span

                    score = granularity_weight * g +
                            simplicity_weight * s +
                            coverage_weight * c +
                            niceness_weight * qscore

                    # strict limits on coverage
                    if strict_span && span > xspan
                        score -= 10000
                    end
                    if span >= 2xspan
                        score -= 1000
                    end

                    if score > high_score && (k_min <= length(S) <= k_max)
                        if strict_span
                            # make S a copy because it is a view and
                            # could otherwise be mutated in the next runs
                            S = collect(S)
                        end
                        S_best, viewmin_best, viewmax_best = S, viewmin, viewmax
                        high_score = score
                    end
                    r += 1
                end
            end
        end
        z -= 1
    end

    if isinf(high_score)
        if strict_span
            @warn("No strict ticks found")
            return optimize_ticks_typed(x_min, x_max, extend_ticks, Q, k_min, k_max, k_ideal,
                                       granularity_weight, simplicity_weight,
                                       coverage_weight, niceness_weight,
                                       false, span_buffer, scale)
        else
            R = typeof(1. * one_t)
            return R[x_min], x_min - one_t, x_min + one_t
        end
    end

    return S_best, viewmin_best, viewmax_best
end

function optimize_ticks(x_min::Date, x_max::Date; extend_ticks::Bool=false,
                        k_min=nothing, k_max=nothing, scale=:auto,
                        granularity_weight=nothing, simplicity_weight=nothing,
                        coverage_weight=nothing, niceness_weight=nothing,
                        strict_span=true, span_buffer = nothing)
    return optimize_ticks(convert(DateTime, x_min), convert(DateTime, x_max),
                          extend_ticks=extend_ticks, scale=scale)
end

function optimize_ticks(x_min::DateTime, x_max::DateTime; extend_ticks::Bool=false,
                        k_min=nothing, k_max=nothing, scale=:auto,
                        granularity_weight=nothing, simplicity_weight=nothing,
                        coverage_weight=nothing, niceness_weight=nothing,
                        strict_span=true, span_buffer = nothing)
    if x_min == x_max
        x_max += Second(1)
    end

    if year(x_max) - year(x_min) <= 1 && scale != :year
        if year(x_max) == year(x_min) && month(x_max) - month(x_min) <= 1 && scale != :month
            ticks = DateTime[]

            scales = [
                Day(1), Hour(1), Minute(1), Second(1), Millisecond(100),
                Millisecond(10), Millisecond(1)
            ]

            # ticks on week boundries
            if x_min + Day(7) < x_max || scale == :week
                push!(ticks, x_min)
                while true
                    next_month = Date(year(ticks[end]), month(ticks[end])) + Month(1)
                    while ticks[end] + Week(1) < next_month - Day(2)
                        push!(ticks, ticks[end] + Week(1))
                    end
                    push!(ticks, next_month)
                    if next_month >= x_max
                        break
                    end
                end
            else
                scale = nothing
                if scale != :auto
                    # TODO: manually setting scale with :day, :minute, etc
                end

                if scale === nothing
                    for proposed_scale in [Day(1), Hour(1), Minute(1),
                                           Second(1), Millisecond(100),
                                           Millisecond(10), Millisecond(1)]
                        if x_min + proposed_scale < x_max
                            scale = proposed_scale
                            break
                        end
                    end
                end

                if scale === nothing
                    scale = Millisecond(1)
                end

                # round x_min down
                if scale === Day(1)
                    first_tick = DateTime(year(x_min), month(x_min), day(x_min))
                elseif scale === Hour(1)
                    first_tick = DateTime(year(x_min), month(x_min), day(x_min),
                                          hour(x_min))
                elseif scale === Minute(1)
                    first_tick = DateTime(year(x_min), month(x_min), day(x_min),
                                          hour(x_min), minute(x_min))
                elseif scale === Second(1)
                    first_tick = DateTime(year(x_min), month(x_min), day(x_min),
                                          hour(x_min), minute(x_min), second(x_min))
                elseif scale === Millisecond(100)
                    first_tick = DateTime(year(x_min), month(x_min), day(x_min),
                                          hour(x_min), minute(x_min),
                                          second(x_min), millisecond(x_min) % 100)
                elseif scale === Millisecond(10)
                    first_tick = DateTime(year(x_min), month(x_min), day(x_min),
                                          hour(x_min), minute(x_min),
                                          second(x_min), millisecond(x_min) % 10)
                else
                    first_tick = x_min
                end
                push!(ticks, first_tick)

                while ticks[end] < x_max
                    push!(ticks, ticks[end] + scale)
                end
            end

            viewmin, viewmax = ticks[1], ticks[end]
            return ticks, viewmin, viewmax
        else
            ticks = DateTime[]
            push!(ticks, Date(year(x_min), month(x_min)))
            while ticks[end] < x_max
                push!(ticks, ticks[end] + Month(1))
            end
            viewmin, viewmax = ticks[1], ticks[end]

            return ticks, x_min, x_max
        end
    else
        ticks, viewmin, viewmax =
            optimize_ticks(year(x_min), year(x_max + Year(1) - Day(1)), extend_ticks=extend_ticks)

        return DateTime[DateTime(round(y)) for y in ticks],
                        DateTime(round(viewmin)), DateTime(round(viewmax))
    end
end

# Generate ticks suitable for multiple scales.
function multilevel_ticks(viewmin::T, viewmax::T;
                          scales=[0.5, 5.0, 10.0]) where T

    ticks = Dict()
    for scale in scales
        ticks[scale] = optimize_ticks(viewmin, viewmax,
                                      k_min=max(1, round(Int, 2*scale)),
                                      k_max=max(3, round(Int, 10*scale)),
                                      k_ideal=max(2, round(Int, 15*scale)))[1]
    end

    return ticks
end

function multilevel_ticks(viewmin::Date, viewmax::Date;
                          scales=[:year, :month, :day])
    return multilevel_ticks(convert(DateTime, viewmin),
                            convert(DateTime, viewmax),
                            scales=scales)
end

function multilevel_ticks(viewmin::DateTime, viewmax::DateTime;
                          scales=[:year, :month, :day])
    # TODO: This needs to be improved for DateTime
    span = convert(Float64, Dates.toms(viewmax - viewmin))
    ticks = Dict()
    for scale in scales
        if scale == :year
            s = span / Dates.toms(Day(360))
        elseif scale == :month
            s = span / Dates.toms(Day(90))
        else
            s = span / Dates.toms(Day(1))
        end

        ticks[s/20] = optimize_ticks(viewmin, viewmax, scale=scale)[1]
    end

    return ticks
end

# Choose "round" (full seconds/minutes/hours/days/months/years) DateTime ticks
# between x_min and x_max:
function optimize_datetime_ticks(a_min, a_max; k_min = 2, k_max = 4)
    x_min = DateTime(Dates.UTM(Int64(round(a_min))))
    x_max = DateTime(Dates.UTM(Int64(round(a_max))))

    Δt = x_max - x_min
    if Δt > Dates.Day(365 * k_min)
        P = Dates.Year
        steplength = Δt / (k_max * Dates.Millisecond(Dates.Day(365)))
    elseif Δt > Dates.Day(30 * k_min)
        P = Dates.Month
        steplength = Δt / (k_max * Dates.Millisecond(Dates.Day(30)))
    elseif Δt > Dates.Day(k_min)
        P = Dates.Day
        steplength = Δt / (k_max * Dates.Millisecond(Dates.Day(1)))
    elseif Δt > Dates.Hour(k_min)
        P = Dates.Hour
        steplength = Δt / (k_max * Dates.Millisecond(Dates.Hour(1)))
    elseif Δt > Dates.Minute(k_min)
        P = Dates.Minute
        steplength = Δt / (k_max * Dates.Millisecond(Dates.Minute(1)))
    elseif Δt > Dates.Second(k_min)
        P = Dates.Second
        steplength = Δt / (k_max * Dates.Millisecond(Dates.Second(1)))
    else
        P = Dates.Millisecond
        steplength = Δt / (k_max * Dates.Millisecond(1))
    end
    steplength = P(max(1, Int(round(steplength))))

    period_hierarchy = [Dates.Millisecond, Dates.Second, Dates.Minute,
        Dates.Hour, Dates.Day, Dates.Month, Dates.Year]
    i = findfirst(period_hierarchy .== P)
    showtype = i >= 5 ? Date : DateTime
    start = DateTime((PH(x_min) for PH in period_hierarchy[end:-1:i])...) + P(1)
    ticks = collect(start:steplength:x_max)
    labels = string.(showtype.(ticks))

    return Dates.value.(ticks), labels
end


text:
import Test
import StaticArrays:
    SVector,
    MVector
using StateSpaceReconstruction: Simplices

# Constructors
@testset "Immutable simplices" begin
    @test Simplex([rand(3) for i = 1:4]) isa Simplex
    @test Simplex(rand(3, 4)) isa Simplex
    @test Simplex(rand(4, 3)) isa Simplex
    @test Simplex(rand(4, 5)) isa Simplex
    @test Simplex(rand(5, 4)) isa Simplex
    @test SSimplex([SVector{3}(rand(3)) for i = 1:4]) isa SSimplex

    @test_throws DomainError Simplex(rand(5, 3))

    s = Simplex([rand(3) for i = 1:4])

    @test npoints(s) == 4
    @test dimension(s) == 3

    @test volume(s) isa Float64
    @test volume(s) >= 0

    @test orientation(s) isa Float64

    s1 = Simplex(rand(3, 4))
    s2 = Simplex(rand(4, 3))

    # Intersection volume
    @test intersect(s1, s2) >= 0

    # Intersecting vertices
    @test intersect_vertices(s1, s2) isa Array{Float64, 2}
end

# Constructors
@testset "Mutable simplices" begin
    @test MutableSimplex([rand(3) for i = 1:4]) isa MutableSimplex
    @test MutableSimplex([rand(3) for i = 1:4]) isa MutableSimplex
    @test MutableSSimplex([MVector{3}(rand(3)) for i = 1:4]) isa MutableSSimplex

    @test MutableSimplex(rand(3, 4)) isa MutableSimplex
    @test MutableSimplex(rand(4, 3)) isa MutableSimplex
    @test MutableSimplex(rand(4, 5)) isa MutableSimplex
    @test MutableSimplex(rand(5, 4)) isa MutableSimplex
    @test_throws DomainError MutableSimplex(rand(5, 3))

    s = MutableSimplex([rand(3) for i = 1:4])

    @test volume(s) isa Float64
    @test volume(s) >= 0

    @test orientation(s) isa Float64

    s1 = MutableSimplex(rand(3, 4))
    s2 = MutableSimplex(rand(4, 3))

    # Intersection volume
    @test intersect(s1, s2) >= 0

    # Intersecting vertices
    @test intersect_vertices(s1, s2) isa Array{Float64, 2}
end

@testset "Volume intersections" begin
    # Generate the vertices of two intersecting simplices, then construct simplices of
    # different types from those vertices. All types should give the same intersecting
    # volume and intersecting vertices.
    vertices_s1 = [rand(3) for i = 1:4]

    s1 = Simplex(vertices_s1)
    s1_mutable = MutableSimplex(vertices_s1)
    #ss1 = SSimplex([SVector{3, Float64}(vertices_s1[i]) for i = 1:4])
    ss1_mutable = MutableSSimplex([MVector{3, Float64}(vertices_s1[i]) for i = 1:4])

    s2 = generate_intersecting_simplex(s1)
    vertices_s2 = s2.vertices

    s2_mutable = MutableSimplex(vertices_s2)
    ss2 = SSimplex([SVector{3, Float64}(vertices_s2[i]) for i = 1:4])
    ss2_mutable = MutableSSimplex([MVector{3, Float64}(vertices_s2[i]) for i = 1:4])

    # Individual intersections should all be nonzero
    @test s1 ∩ s2 > 1e-12
    @test s1 ∩ s2_mutable > 1e-12
    @test s1 ∩ ss2 > 1e-12
    @test s1 ∩ ss2_mutable > 1e-12

    # Using the different types should give the same answer. Just try a few different
    # combinations.
    @test s1 ∩ s2 == s1 ∩ s2
    @test s1 ∩ s2 == s1 ∩ ss2
    @test s1 ∩ s2 == s1 ∩ s2_mutable
    @test s1 ∩ s2 == s1 ∩ ss2_mutable

end

@testset "Does a point lie inside simplex?" begin
    s = Simplex(rand(3, 4))
    n = 10000
    contained = Vector{Bool}(undef, n)

    # Interior points
    for i = 1:n
        p = generate_interior_points(s, 1)[1]
        contained[i] = point_contained(s, p)
    end

    @test all(x -> (x .== true), contained)

    # Exterior points
    for i = 1:n
        p = generate_exterior_points(s, 1)[1]
        contained[i] = point_contained(s, p)
    end

    @test all(x -> (x .== false), contained)
end

@testset "Refining simplices" begin

    pts = [rand(3) for i = 1:20]
    inv_pts = invariantize(pts);
    triang = triangulate(inv_pts[1:end - 1])
    triang_im = triangulate(inv_pts[2:end])

    simplices = [Simplex(inv_pts[triang.simplexindices[i]]) for i in 1:length(triang.simplexindices)];
    simplices_im = [Simplex(inv_pts[triang_im.simplexindices[i]]) for i in 1:length(triang_im.simplexindices)];

    @test refine(simplices[1]) isa Vector{Simplex{3,Float64}}
    @test refine(simplices[1], simplices[2]) isa Tuple{Vector{Simplex{3,Float64}}, Vector{Simplex{3,Float64}}}

    k = 2
    dim = 3
    rules = simplex_splitting_rules(k, dim)
    @test refine(simplices[1], k = 2, splitting_rules = rules) isa Vector{Simplex{3,Float64}}
    @test refine(simplices[1], simplices[2], k = 2, splitting_rules = rules) isa Tuple{Vector{Simplex{3,Float64}}, Vector{Simplex{3,Float64}}}
end

@testset "Subsampling" begin
    dim = 3
    inv_pts = invariantize([rand(dim) for i = 1:20])
    triang = triangulate(inv_pts[1:end-1])
    triang_im = triangulate(inv_pts[2:end])
    simplices = [Simplex(inv_pts[triang.simplexindices[i]]) for i in 1:length(triang.simplexindices)];
    simplices_im = [Simplex(inv_pts[triang_im.simplexindices[i]]) for i in 1:length(triang_im.simplexindices)];

    coeffs = subsample_coeffs(dim, 50, false)

    # On single simplices.
    @test subsample(simplices[1], coeffs) isa Vector{SVector{3, Float64}}
    @test subsample(simplices[1]) isa Vector{SVector{3, Float64}}

    # Pairwise
    ss_simplices_coeffs_precomputed = subsample(simplices[1], simplices_im[1], coeffs)
    ss_simplices_coeffs_computed_onthefly = subsample(simplices[1], simplices_im[1])
    @test ss_simplices_coeffs_precomputed isa Tuple{Vector{SVector{3, Float64}}, Vector{SVector{3, Float64}}}
    @test ss_simplices_coeffs_computed_onthefly isa Tuple{Vector{SVector{3, Float64}}, Vector{SVector{3, Float64}}}
end


text:
using UUIDs
import Sockets
is_ch_type(::Val{:Nullable})  = true
can_be_nullable(::Val{:Nullable}) = false

missings_vector_type(::Type{Vector{T}}) where {T} =
                                                 Vector{Union{T, Missing}}
missings_vector_type(::Type{CategoricalVector{T}}) where {T} =
                                                CategoricalVector{Union{T, Missing}}

result_type(::Val{:Nullable}, nested)  = missings_vector_type(result_type(nested))

convert_to_missings(data::Vector{T}) where {T} =
                                 convert(Vector{Union{T, Missing}}, data)

convert_to_missings(data::CategoricalVector{T}) where {T} =
                            convert(CategoricalVector{Union{T, Missing}}, data)

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt,
                                            ::Val{:Nullable}, nested::TypeAst)

    missing_map = chread(sock, Vector{UInt8}, num_rows)
    unmissing = read_col_data(sock, num_rows, nested)
    result = convert_to_missings(unmissing)
    for i in 1:length(missing_map)
        (missing_map[i] == 0x1) && (result[i] = missing)
    end
    return result
end

missing_replacement(::Type{T}) where {T <: Number} = zero(T)
missing_replacement(::Type{UUID}) = UUID(0)
missing_replacement(::Type{Date}) = Date(1970)
missing_replacement(::Type{DateTime}) = unix2datetime(0)
missing_replacement(::Type{String}) = ""
missing_replacement(::Type{Sockets.IPv4}) = Sockets.IPv4(0)
missing_replacement(::Type{Sockets.IPv6}) = Sockets.IPv6(0)

uint8_ismissing(v)::UInt8 = ismissing(v) ? 1 : 0

function write_col_data(sock::ClickHouseSock,
                                data::AbstractVector{Union{Missing, T}},
                                ::Val{:Nullable}, nested::TypeAst) where {T}
    !can_be_nullable(nested.name) &&
            error("$(nested.name) cannot be inside Nullable")
    missing_map = uint8_ismissing.(data)
    chwrite(sock, missing_map)
    unmissing = if !any(x -> x > 0, missing_map)
        convert(Vector{T}, data)
    else
        replacement = missing_replacement(T)
        [ismissing(v) ? replacement : v for v in data]
    end

    write_col_data(sock, unmissing, nested)
end

function write_col_data(sock::ClickHouseSock,
                                data::AbstractVector{T},
                                ::Val{:Nullable}, nested::TypeAst) where {T}
    !can_be_nullable(nested.name) &&
            error("$(nested.name) cannot be inside Nullable")

    missing_map = fill(Int8(0), length(data))
    chwrite(sock, missing_map)
    write_col_data(sock, data, nested)
end

function write_col_data(sock::ClickHouseSock,
                                data::AbstractCategoricalVector{Union{Missing, T}},
                                ::Val{:Nullable}, nested::TypeAst) where {T}
    !can_be_nullable(nested.name) &&
            error("$(nested.name) cannot be inside Nullable")
    missing_map = uint8_ismissing.(data)
    chwrite(sock, missing_map)
    unmissing = if !any(x -> x > 0, missing_map)
        convert(CategoricalVector{T}, data)
    else
        tmp = deepcopy(data)
        #replace missing (it's always 0 in refs of CategorialVector)
        #with something valid
        replace!(tmp.refs, 0=>1)
        convert(CategoricalVector{T}, tmp)
    end

    write_col_data(sock, unmissing, nested)
end

#For Nothing columns,
#(i.e. column of nulls in CHtypes of CH, of missing column in Julia types)
function write_col_data(sock::ClickHouseSock,
                                data::AbstractVector{Missing},
                                ::Val{:Nullable}, nested::TypeAst) where {T}
    nested.name != :Nothing &&
        error("Vector{Missing} can be writen only to Nullable(Nothing) column")
    missing_map = Vector{UInt8}(undef, length(data))
    fill!(missing_map, 0x1)
    chwrite(sock, missing_map)
    write_col_data(sock, data, nested)
end


text:
export Block
import Base: getindex, setindex!, length, size, firstindex, lastindex, push!, pop!, pushfirst!, popfirst!, append!, prepend!

struct Block
    layers
    name                        # only used for display
    Block(ls...; name=nothing) = new(Any[ls...], name)
end

function (s::Block)(x)
    for l in s.layers
        x = l(x)
    end
    return x
end

# AbstractArray interface
function getindex(s::Block, i)
    if length(i) == 1
        getindex(s.layers, i)
    else
        name=(s.name === nothing ? nothing : "$(s.name)[$i]")
        Block(getindex(s.layers, i)...; name)
    end
end

setindex!(s::Block, v, i) = setindex!(s.layers, v, i)
length(s::Block) = length(s.layers)
size(s::Block) = size(x.layers)
firstindex(s::Block) = 1
lastindex(s::Block) = length(s)

# Dequeue interface
push!(s::Block, x...) = (push!(s.layers, x...); s)
pop!(s::Block) = pop!(s.layers)
pushfirst!(s::Block, x...) = (pushfirst!(s.layers, x...); s)
popfirst!(s::Block) = popfirst!(s.layers)
append!(s::Block, t::Block) = (append!(s.layers, t.layers); s)
append!(s::Block, t) = (append!(s.layers, t); s)
prepend!(s::Block, t::Block) = (prepend!(s.layers, t.layers); s)
prepend!(s::Block, t) = (prepend!(s.layers, t); s)

# TODO: iterator interface


text:
# This file was generated, do not modify it. # hide
N=10; @show collect(range(-2,2,length=N))