text:
@testset "BLAS" begin
    @testset "dot" begin
        @testset "all entries" begin
            n = 10
            test_frule(BLAS.dot, randn(n), randn(n))
            test_rrule(BLAS.dot, randn(n), randn(n))
        end

        @testset "over strides" begin
            n = 10
            incx = 2
            incy = 3
            test_rrule(
                BLAS.dot,
                n ⊢ nothing,
                randn(n * incx),
                incx ⊢ nothing,
                randn(n * incy),
                incy ⊢ nothing,
            )
        end
    end

    @testset "nrm2" begin
        @testset "all entries" begin
            @testset "$T" for T in (Float64, ComplexF64)
                n = 10
                test_frule(BLAS.nrm2, randn(T, n))
                test_rrule(BLAS.nrm2, randn(T, n))
            end
        end

        @testset "over strides" begin
            dims = (3, 2, 1)
            incx = 2
            @testset "Array{$T,$N}" for N in eachindex(dims), T in (Float64, ComplexF64)
                s = (dims[1] * incx, dims[2:N]...)
                n = div(prod(s), incx)
                test_rrule(
                    BLAS.nrm2, n ⊢ nothing, randn(T, s), incx ⊢ nothing; atol=0, rtol=1e-5,
                )
            end
        end
    end

    @testset "asum" begin
        @testset "all entries" begin
            @testset "$T" for T in (Float64, ComplexF64)
                n = 6
                test_frule(BLAS.asum, randn(T, n))
                test_rrule(BLAS.asum, randn(T, n))
            end
        end

        @testset "over strides" begin
            dims = (2, 2, 1)
            incx = 2
            @testset "Array{$T,$N}" for N in eachindex(dims), T in (Float64, ComplexF64)
                s = (dims[1] * incx, dims[2:N]...)
                n = div(prod(s), incx)
                test_rrule( BLAS.asum, n ⊢ nothing, randn(T, s), incx ⊢ nothing)
            end
        end
    end

    @testset "gemm" begin
        for m in 3:5, n in 3:5, p in 3:5, tA in ('N', 'C', 'T'), tB in ('N', 'C', 'T'), T in (Float64, ComplexF64)
            A = randn(T, tA === 'N' ? (m, n) : (n, m))
            B = randn(T, tB === 'N' ? (n, p) : (p, n))
            test_rrule(gemm, tA ⊢ nothing, tB ⊢ nothing, A, B; check_inferred=false)
            test_rrule(  # 5 arg version with scaling scalar
                gemm, tA ⊢ nothing, tB ⊢ nothing, randn(T), A, B; check_inferred=false,
            )
        end
    end

    @testset "gemv" begin
        for n in 3:5, m in 3:5, t in ('N', 'C', 'T'), T in (Float64, ComplexF64)
            x = randn(T, t === 'N' ? n : m)
            test_rrule(gemv, t ⊢ nothing, randn(T), randn(T, m, n), x; check_inferred=false)
        end
    end
end


text:
module GenBank

using TranscodingStreams
using BioSequences
using BioGenerics
using DataFrames
using CodecZlib

import ..GenomicAnnotations: Record, Gene, AbstractGene, GeneDataView, Locus
import ..GenomicAnnotations: sequence, iscomplement, iscomplete, addgene!, pushproperty!, feature, index, locus, oneline
export sequence, iscomplement, iscomplete, feature, index, locus

include("reader.jl")
include("writer.jl")

end #module

ext: jl


text:
using DrWatson
@quickactivate "StatReth"

# %%
using DataFrames
using CSV
using Turing
using Plots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %% 5.28, 5.29
d = DataFrame!(CSV.File(datadir("exp_raw/milk.csv"), missingstring = "NA"))
# We get the "confusing error message" even earlier because zscore doesn't know how to work
# with missing data, so we'll drop them right out the gate.
dropmissing!(d, [:kcal_per_g, :neocortex_perc, :mass])

d.K = zscore(d.kcal_per_g)
d.N = zscore(d.neocortex_perc)
d.M = zscore(log.(d.mass))

# %% 5.30
# This is why this *would* work on the first try, so we are skipping this

# %% 5.31, 5.32
d.neocortex_perc
# Already did this so it's not going to do anything additional here.
dropmissing!(d, [:kcal_per_g, :neocortex_perc, :mass])

# %% 5.33
@model function milk1(N, K)
    a ~ Normal(0, 1)
    bN ~ Normal(0, 1)
    σ ~ Exponential(1)
    μ = lin(a, bN, N)
    K ~ MvNormal(μ, σ)
end

m5_5_draft = milk1(d.N, d.K)
q5_5_draft = quap(m5_5_draft)

# %% 5.34
prior = sample(m5_5_draft, Prior(), 50) |> DataFrame
xseq = -2:0.1:2
μ = lin(prior.a', prior.bN', xseq)

plot(legend = false, ylims = extrema(xseq))
for c in eachcol(μ)
    plot!(xseq, c, color = :black, alpha = 0.3)
end
plot!()

# %% 5.35
@model function milk2(N, K)
    a ~ Normal(0, 0.2)
    bN ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, bN, N)
    K ~ MvNormal(μ, σ)
end

m5_5 = milk2(d.N, d.K)
q5_5 = quap(m5_5)

# %%
prior = sample(m5_5, Prior(), 50) |> DataFrame
xseq = -2:0.1:2
μ = lin(prior.a', prior.bN', xseq)

plot(legend = false, ylims = extrema(xseq))
for c in eachcol(μ)
    plot!(xseq, c, color = :black, alpha = 0.3)
end
plot!()

# %% 5.36
DataFrame(rand(q5_5.distr, 1000)', q5_5.params) |> precis

# %% 5.37
xseq = range(minimum(d.N) - 0.15, maximum(d.N) + 0.15, length = 30)
post = DataFrame(rand(q5_5.distr, 1000)', q5_5.params)
μ = lin(post.a', post.bN', xseq) |> meanlowerupper

scatter(d.N, d.K, legend = false)
plot!(xseq, μ.mean, ribbon = (μ.mean .- μ.lower, μ.upper .- μ.mean))

# %% 5.38
@model function milk3(M, K)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, bM, M)
    K ~ MvNormal(μ, σ)
end

m5_6 = milk3(d.M, d.K)
q5_6 = quap(m5_6)
post = DataFrame(rand(q5_6.distr, 1000)', q5_6.params)
precis(post)

# %% 5.39
@model function milk4(M, N, K)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    bN ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, bM, M, bN, N)
    K ~ MvNormal(μ, σ)
end

q5_7 = quap(milk4(d.M, d.N, d.K))
post = DataFrame(rand(q5_7.distr, 1000)', q5_7.params)
precis(post)

# %% 5.40
# TODO coefplot

# %% 5.41
xseq = range(minimum(d.M) - 0.15, maximum(d.M) + 0.15, length = 30)
mu = lin(post.a', post.bM', xseq, post.bN', zeros(length(xseq))) |> meanlowerupper
plot(xseq, mu.mean, ribbon = (mu.upper .- mu.mean, mu.mean .- mu.lower))

# %% 5.42
# M -> K <- N
# M -> N
n = 100
M = randn(n)
N = rand.(Normal.(M))
K = rand.(Normal.(N .- M))
d_sim = DataFrame((; K, N, M))

# %% 5.43
# M -> K <- N
# N -> M
n = 100
N = randn(n)
M = rand.(Normal.(N))
K = rand.(Normal.(N .- M))
d_sim2 = DataFrame((; K, N, M))

# M -> K <- N
# M -> U <- N
n = 100
U = randn(n)
N = rand.(Normal.(U))
M = rand.(Normal.(U))
K = rand.(Normal.(N .- M))
d_sim2 = DataFrame((; K, N, M))

# %% 5.44
# TODO working with DAGs

max_stars_repo_path: scripts/5_2 Masked relationship.jl

max_issues_repo_path: scripts/5_2 Masked relationship.jl

max_forks_repo_path: scripts/5_2 Masked relationship.jl


text:
"""
TODO
"""
struct AnalyticalVelocityAndStress <: InitializationStrategy end
function initial_condition(
    ::AnalyticalVelocityAndStress,
    q::Quadrature,
    problem::FluidFlowProblem,
    x::T,
    y::T,
) where {T <: Real}
    u = lattice_velocity(q, problem, x, y)
    ρ = 1.0
    f = equilibrium(q, one(T), u, one(T))
    σ = problem.u_max^2 * deviatoric_tensor(q, problem, x, y, 0.0)
    ∇u = problem.u_max^2 * velocity_gradient(problem, x, y, 0.0)

    cs = q.speed_of_sound_squared
    τ = cs * lattice_viscosity(problem)
    τ_eff = τ + 0.5
    for f_idx in 1:length(f)
        w_i = q.weights[f_idx]
        H_2 = hermite(Val{2}, q.abscissae[:, f_idx], q)

        # TODO: possibly replace with pressure
        temperature = 1.0
        f_neq = w_i * (cs * τ_eff * ρ * temperature) / (2) * dot(H_2, ∇u + ∇u')
        f[f_idx] += -f_neq
    end
    return f
end


text:
using Router, Channels

channel("/tweets/subscribe", "tweets#TweetsChannel.subscribe")

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
using NonlinearSchrodinger
using Documenter

makedocs(;
    modules=[NonlinearSchrodinger],
    authors="Omar Ashour",
    repo="https://github.com/oashour/NonlinearSchrodinger.jl/blob/{commit}{path}#L{line}",
    sitename="NonlinearSchrodinger.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="http://oashour.github.io",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Theory" => [
            "Simulations" => "man/simulations.md",
            "Darboux Transformations" =>  "man/dt.md"],
        "Examples" => [
            "Simulation Examples" => "man/sim_examples.md",
            "DT Examples" =>  "man/dt_examples.md"],
        "Visualization" => "man/visualization.md",
        "Library" => [
            "Public API" => "public.md",
            "Internals" =>  "private.md",
            "Integrators" => "man/integrators.md"]
    ],
)

deploydocs(;
    repo="github.com/oashour/NonlinearSchrodinger.jl",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev"
    #versions = ["stable" => "v^", "v#.#", devurl => devurl],
)


text:
include("test_params.jl")
using SymbolicRegression, SymbolicUtils
using Test
using SymbolicRegression: string_tree
using Random

x1 = 0.1f0;
x2 = 0.1f0;
x3 = 0.1f0;
x4 = 0.1f0;
x5 = 0.1f0;
for i in 0:5
    batching = i in [0, 1]
    weighted = i in [0, 2]

    numprocs = 4
    progress = false
    warmupMaxsizeBy = 0.0f0
    optimizer_algorithm = "NelderMead"
    multi = false
    probPickFirst = 1.0
    multithreading = false
    crossoverProbability = 0.0f0
    skip_mutation_failures = false
    useFrequency = false
    useFrequencyInTournament = false
    print("Testing with batching=$(batching) and weighted=$(weighted), ")
    if i == 0
        println("with serial & progress bar & warmup & BFGS")
        numprocs = 0 #Try serial computation here.
        progress = true #Also try the progress bar.
        warmupMaxsizeBy = 0.5f0 #Smaller maxsize at first, build up slowly
        optimizer_algorithm = "BFGS"
        probPickFirst = 0.8
    elseif i == 1
        println("with multi-output and useFrequency.")
        multi = true
        useFrequency = true
    elseif i == 3
        println("with multi-threading and crossover and useFrequencyInTournament")
        multithreading = true
        numprocs = 0
        crossoverProbability = 0.02f0
        useFrequencyInTournament = true
    elseif i == 4
        println("with crossover and skip mutation failures and both frequencies options")
        crossoverProbability = 0.02f0
        skip_mutation_failures = true
        useFrequency = true
        useFrequencyInTournament = true
    elseif i == 5
        println("with default hyperparameters")
    end
    if i == 5
        options = SymbolicRegression.Options(;
            unary_operators=(cos,),
            batching=batching,
            parsimony=0.0f0, # Required for scoring
        )
    else
        options = SymbolicRegression.Options(;
            default_params...,
            binary_operators=(+, *),
            unary_operators=(cos,),
            npopulations=4,
            batching=batching,
            crossoverProbability=crossoverProbability,
            skip_mutation_failures=skip_mutation_failures,
            seed=0,
            progress=progress,
            warmupMaxsizeBy=warmupMaxsizeBy,
            optimizer_algorithm=optimizer_algorithm,
            probPickFirst=probPickFirst,
            parsimony=0.0f0,
            useFrequency=useFrequency,
            useFrequencyInTournament=useFrequencyInTournament,
        )
    end

    X = randn(MersenneTwister(0), Float32, 5, 100)
    if weighted
        mask = rand(100) .> 0.5
        weights = map(x -> convert(Float32, x), mask)
        # Completely different function superimposed - need
        # to use correct weights to figure it out!
        y = (2 .* cos.(X[4, :])) .* weights .+ (1 .- weights) .* (5 .* X[2, :])
        hallOfFame = EquationSearch(
            X,
            y;
            weights=weights,
            niterations=2,
            options=options,
            numprocs=numprocs,
            multithreading=multithreading,
        )
        dominating = [calculate_pareto_frontier(X, y, hallOfFame, options; weights=weights)]
    else
        y = 2 * cos.(X[4, :])
        if multi
            # Copy the same output twice; make sure we can find it twice
            y = repeat(y, 1, 2)
            y = transpose(y)
        end
        hallOfFame = EquationSearch(
            X, y; niterations=2, options=options, multithreading=multithreading
        )
        if multi
            dominating = [
                calculate_pareto_frontier(X, y[j, :], hallOfFame[j], options) for j in 1:2
            ]
        else
            dominating = [calculate_pareto_frontier(X, y, hallOfFame, options)]
        end
    end

    # Always assume multi
    for dom in dominating
        best = dom[end]
        eqn = node_to_symbolic(best.tree, options)

        local x4 = SymbolicUtils.Sym{Real}(Symbol("x4"))
        true_eqn = 2 * cos(x4)
        residual = simplify(eqn - true_eqn) + x4 * 1e-10

        # Test the score
        @test best.loss < maximum_residual / 10
        # Test the actual equation found:
        # eval evaluates inside global
        @test abs(eval(Meta.parse(string(residual)))) < maximum_residual
    end

    println("Passed.")
end # for i=1...

println("Testing fast-cycle and custom variable names, with mutations")

options = SymbolicRegression.Options(;
    default_params...,
    binary_operators=(+, *),
    unary_operators=(cos,),
    npopulations=4,
    constraints=((*) => (-1, 10), cos => (5)),
    fast_cycle=true,
    skip_mutation_failures=true,
    stateReturn=true,
)
X = randn(MersenneTwister(0), Float32, 5, 100)
y = 2 * cos.(X[4, :])
varMap = ["t1", "t2", "t3", "t4", "t5"]
state, hallOfFame = EquationSearch(X, y; varMap=varMap, niterations=2, options=options)
dominating = calculate_pareto_frontier(X, y, hallOfFame, options)

best = dominating[end]

eqn = node_to_symbolic(best.tree, options; varMap=varMap)

t4 = SymbolicUtils.Sym{Real}(Symbol("t4"))
true_eqn = 2 * cos(t4)
residual = simplify(eqn - true_eqn) + t4 * 1e-10

# Test the score
@test best.loss < maximum_residual / 10
# Test the actual equation found:
t1 = 0.1f0;
t2 = 0.1f0;
t3 = 0.1f0;
t4 = 0.1f0;
t5 = 0.1f0;
residual_value = abs(eval(Meta.parse(string(residual))))
@test residual_value < maximum_residual

# Do search again, but with saved state:
# We do 0 iterations to make sure the state is used.
println("Passed.")
println("Testing whether state saving works.")
state, hallOfFame = EquationSearch(
    X, y; varMap=varMap, niterations=0, options=options, saved_state=(state, hallOfFame)
)

dominating = calculate_pareto_frontier(X, y, hallOfFame, options)
best = dominating[end]
print_tree(best.tree, options)
eqn = node_to_symbolic(best.tree, options; varMap=varMap)
residual = simplify(eqn - true_eqn) + t4 * 1e-10
@test best.loss < maximum_residual / 10

println("Passed.")

println("Testing whether we can stop based on clock time.")
options = Options(; default_params..., timeout_in_seconds=1)
start_time = time()
EquationSearch(X, y; niterations=10000000, options=options, multithreading=true)
end_time = time()
@test end_time - start_time < 100
println("Passed.")

println("Running README example.")
include("../example.jl")
println("Passed.")

println("Testing whether the recorder works.")
include("test_recorder.jl")
println("Passed.")

println("Testing whether deterministic mode works.")
include("test_deterministic.jl")
println("Passed.")


text:
# ---
# title: 150. Evaluate Reverse Polish Notation
# id: problem150
# author: Indigo
# date: 2021-01-19
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/evaluate-reverse-polish-notation/description/>
# hidden: true
# ---
#
# Evaluate the value of an arithmetic expression in [Reverse Polish
# Notation](http://en.wikipedia.org/wiki/Reverse_Polish_notation).
#
# Valid operators are `+`, `-`, `*`, `/`. Each operand may be an integer or
# another expression.
#
# **Note:**
#
#   * Division between two integers should truncate toward zero.
#   * The given RPN expression is always valid. That means the expression would always evaluate to a result and there won't be any divide by zero operation.
#
# **Example 1:**
#
#
#
#     Input: ["2", "1", "+", "3", "*"]
#     Output: 9
#     Explanation: ((2 + 1) * 3) = 9
#
#
# **Example 2:**
#
#
#
#     Input: ["4", "13", "5", "/", "+"]
#     Output: 6
#     Explanation: (4 + (13 / 5)) = 6
#
#
# **Example 3:**
#
#
#
#     Input: ["10", "6", "9", "3", "+", "-11", "*", "/", "*", "17", "+", "5", "+"]
#     Output: 22
#     Explanation:
#       ((10 * (6 / ((9 + 3) * -11))) + 17) + 5
#     = ((10 * (6 / (12 * -11))) + 17) + 5
#     = ((10 * (6 / -132)) + 17) + 5
#     = ((10 * 0) + 17) + 5
#     = (0 + 17) + 5
#     = 17 + 5
#     = 22
#
#
#
## @lc code=start
using LeetCode
function eval_rpn(tokens::Vector{String})
    stk = Stack{Int}()
    for token in tokens
        if token == "+"
            num2 = pop!(stk)
            num1 = pop!(stk)
            push!(stk, num1 + num2)
        elseif token == "-"
            num2 = pop!(stk)
            num1 = pop!(stk)
            push!(stk, num1 - num2)
        elseif token == "*"
            num2 = pop!(stk)
            num1 = pop!(stk)
            push!(stk, num1 * num2)
        elseif token == "/"
            num2 = pop!(stk)
            num1 = pop!(stk)
            push!(stk, num1 ÷ num2)
        else
            push!(stk, parse(Int, token))
        end
    end
    first(stk)
end

## @lc code=end


text:
export COFFHeader

import Base: show

@io struct COFFHeader
    Machine::UInt16
    NumberOfSections::UInt16
    TimeDateStamp::UInt32
    PointerToSymbolTable::UInt32
    NumberOfSymbols::UInt32
    SizeOfOptionalHeader::UInt16
    Characteristics::UInt16
end

function coff_header_is64bit(h::COFFHeader)
    wide_machines = [
        IMAGE_FILE_MACHINE_AMD64,
        IMAGE_FILE_MACHINE_ARM64,
        IMAGE_FILE_MACHINE_IA64,
    ]
    return h.Machine in wide_machines
end

function coff_header_endianness(h::COFFHeader)
    # We don't ever process :BigEndian COFF files; but if we did, we'd probably
    # design a check against `h.Machine` here.
    return :LittleEndian
end

isexecutable(h::COFFHeader) = (h.Characteristics & IMAGE_FILE_EXECUTABLE_IMAGE) != 0
islibrary(h::COFFHeader) = (h.Characteristics & IMAGE_FILE_DLL) != 0
isrelocatable(h::COFFHeader) = !isexecutable(h) && !islibrary(h)

num_sections(h::COFFHeader) = h.NumberOfSections
num_symbols(h::COFFHeader) = h.NumberOfSymbols

function show(io::IO, header::COFFHeader)
    println(io, "COFF Header")
    println(io, "  Machine ", header.Machine)
    println(io, "  Number of Sections ", header.NumberOfSections)
    println(io, "  Time Date Stamp ", header.TimeDateStamp)
    println(io, "  Number of Symbols ", header.NumberOfSymbols)
end


text:
# This script parses and plots an abstract time series from a ROS bag where
# messages are associated with an iteration number rather than in real time.
# Bags are loaded from a given directory and *should be paired with yaml files
# of the same name.* This script then parses the yaml files to identify ranges
# of values for given params and succinctly selects and processes sets of
# repeated trials.

using PyPlot
using RosDataProcess

# Configuration
save_output = false
colors = generate_colors(10)

println("Loading bag directory")
bag_directory = "$(homedir())/bagfiles/decentralized_exploration_journal/"
bags = load_directory(bag_directory)

# show topics
println("Topic names:")
foreach(println, get_topic_names(bags[1]))

params = Dict(:num_robots => "num_robots",
              :num_decentralized =>
              "kinematic_exploration/num_decentralized_planning_rounds")

ranges = Dict(key => get_range(value, bags) for (key, value) in params)

for num_robots = ranges[:num_robots]
  figure()

  for (ii, num_decentralized) = enumerate(ranges[:num_decentralized])
    println("Plotting trials:"
            * " num_robots=$(num_robots)"
            * " num_decentralized=$(num_decentralized)")

    println("  Filtering")
    fixed = Dict(params[:num_robots] => num_robots,
                 params[:num_decentralized] => num_decentralized)
    trials = filter_bags(fixed, bags)

    println("  Reading bags")
    entropy = read_series("/kinematic_exploration/iteration",
                          "/kinematic_exploration/entropy_reduction", trials;
                          access_data=x->x.data,
                          access_time=x-> num_robots * x.data,
                          intersect=true)

    println("  Plotting data")
    plots = plot_trials(entropy, mean=true, standard_error=true,
                        trials=false, color=colors[ii])

    legend_string(n) = (n == 0 ? "G" : "D\$_$(n)\$")
    plots[:mean][1].set_label(legend_string(num_decentralized))
  end

  title("Entropy reduction: $(num_robots) robots")
  xlabel("Robot-Iteration")
  ylabel("Objective (bits)")
  legend(loc="lower right")

  if save_output
    save_latex("./fig", "entropy_reduction_$(num_robots)")
  end
end

foreach(close, bags)

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


text:
"""
    Kang07(X_cr::T=0.05, η_max::T=0.0348) where T

Efficiency model by Kang, Ryu, Cen, Ostriker 2007, http://arxiv.org/abs/0704.1521v1

## Values
- `X_cr`: P_cr / P_th defined in model for re-acceleration. Basis for interpolation between acceleration and re-acceleration efficiency.
- `η_max`: Maximum efficiency defined in the model
"""
struct Kang07{T} <: AbstractShockAccelerationEfficiency
    X_cr::T
    η_max::T

    Kang07(X_cr::T = 0.3, η_max::T = 0.57) where {T} = new{T}(X_cr, η_max)
end

"""
    η_Ms_acc(η_model::Kang07, M::Real)

Initial acceleration efficiency for model from Kang, Ryu, Cen, Ostriker 2007, http://arxiv.org/abs/0704.1521v1
"""
function η_Ms_acc(η_model::Kang07, M::Real)
    if M < 1
        return 0
    elseif M <= 2
        return 1.96e-3 * (M^2 - 1)             # eq. A3
    else
        return kr_fitting_function(M, 5.46, -9.78, 4.17, -0.334, 0.57)
    end
end

"""
    η_Ms_reacc(η_model::Kang07, M::Real)

Reacceleration efficiency for model from Kang, Ryu, Cen, Ostriker 2007, http://arxiv.org/abs/0704.1521v1
"""
function η_Ms_reacc(η_model::Kang07, M::Real)
    if M <= 1.5
        # adiabatic index of gas
        γ = 5 / 3
        m2 = M^2
        # analytic compression ratio
        xs = (γ + 1) / (γ - 1 + 2 / m2)
        δ0 = 2 * ((2γ * m2 - γ + 1.0) / (γ + 1) - (xs^γ)) /
                (γ * (γ - 1) * m2 * xs)

        return 1.025 * δ0
    else
        return kr_fitting_function(M, 0.24, -1.56, 2.8, 0.512, 0.557)
    end
end


text:
module TestConstrainingLens

include("preamble.jl")

@testset begin
    obj = (a = 1, b = 1)
    constraint = constraining() do obj
        @set obj.b = obj.a
    end
    lens = constraint ∘ @lens _.a

    @test get(obj, lens) === 1
    @test set(obj, lens, 2) === (a = 2, b = 2)
end

@testset "show" begin
    @test string(constraining(identity)) == "constraining(identity)"
    @test string(constraining(
        identity;
        onget = false,
    )) == "constraining(identity; onget=false)"
    @test string(constraining(
        identity;
        onset = Val(false),
    )) == "constraining(identity; onset=false)"
end

end  # module


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractSubProblem{T <: AbstractFloat} end
abstract type AbstractSubProblemState end

struct Gradient{T <: AbstractFloat, A <: AbstractVector}
    δQ::A
    Q::T
    id::Int

    function Gradient(δQ::AbstractVector, Q::AbstractFloat, id::Int)
        T = promote_type(eltype(δQ), Float32)
        δQ_ = convert(AbstractVector{T}, δQ)
        new{T, typeof(δQ_)}(δQ_, Q, id)
    end
end
SparseGradient{T <: AbstractFloat} = Gradient{T, SparseVector{T,Int64}}
DenseGradient{T <: AbstractFloat} = Gradient{T, Vector{T}}

"""
    SubProblem

Quasi-gradient subproblem.

"""
struct SubProblem{T <: AbstractFloat} <: AbstractSubProblem{T}
    id::Int
    probability::T
    model::JuMP.Model
    optimizer::MOI.AbstractOptimizer
    linking_constraints::Vector{MOI.ConstraintIndex}
    masterterms::Vector{Vector{Tuple{Int, Int, T}}}

    function SubProblem(model::JuMP.Model,
                        id::Integer,
                        π::AbstractFloat)
        T = typeof(π)
        # Get optimizer backend
        optimizer = backend(model)
        # Collect all constraints with known decision occurances
        constraints, terms =
            collect_linking_constraints(model,
                                        T)
        subproblem =  new{T}(id,
                             π,
                             model,
                             optimizer,
                             constraints,
                             terms)
        return subproblem
    end
end

function collect_linking_constraints(model::JuMP.Model,
                                     ::Type{T}) where T <: AbstractFloat
    linking_constraints = Vector{MOI.ConstraintIndex}()
    masterterms = Vector{Vector{Tuple{Int, Int, T}}}()
    master_indices = index.(all_known_decision_variables(model, 1))
    # Parse single rows
    F = AffineDecisionFunction{T}
    for S in [MOI.EqualTo{Float64}, MOI.LessThan{Float64}, MOI.GreaterThan{Float64}]
        for ci in MOI.get(backend(model), MOI.ListOfConstraintIndices{F,S}())
            coeffs = Vector{Tuple{Int, Int, T}}()
            f = MOI.get(backend(model), MOI.ConstraintFunction(), ci)::AffineDecisionFunction{T}
            aff = JuMP.jump_function(model, f)::DecisionAffExpr{T}
            for (coef, kvar) in linear_terms(aff.decisions)
                # Map known decisions to master decision,
                # assuming sorted order
                if state(kvar) == Known
                    col = master_indices[index(kvar).value].value
                    push!(coeffs, (1, col, T(coef)))
                end
            end
            if !isempty(coeffs)
                push!(masterterms, coeffs)
                push!(linking_constraints, ci)
            end
        end
    end
    # Parse vector rows
    F = VectorAffineDecisionFunction{T}
    for S in [MOI.Zeros, MOI.Nonpositives, MOI.Nonnegatives]
        for ci in MOI.get(backend(model), MOI.ListOfConstraintIndices{F,S}())
            coeffs = Vector{Tuple{Int, Int, T}}()
            f = MOI.get(backend(model), MOI.ConstraintFunction(), ci)::VectorAffineDecisionFunction{T}
            affs = JuMP.jump_function(model, f)::Vector{DecisionAffExpr{T}}
            for (row, aff) in enumerate(affs)
                for (coef, kvar) in linear_terms(aff.decisions)
                    # Map known decisions to master decision,
                    # assuming sorted order
                    if state(kvar) == Known
                        col = master_indices[index(kvar).value].value
                        push!(coeffs, (row, col, T(coef)))
                    end
                end
            end
            if !isempty(coeffs)
                push!(masterterms, coeffs)
                push!(linking_constraints, ci)
            end
        end
    end
    return linking_constraints, masterterms
end

function update_subproblem!(subproblem::AbstractSubProblem)
    update_known_decisions!(subproblem.optimizer)
    return nothing
end

function (subproblem::SubProblem)(x::AbstractVector)
    return solve_subproblem(subproblem, x)
end

function solve_subproblem(subproblem::SubProblem, x::AbstractVector)
    MOI.optimize!(subproblem.optimizer)
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if status ∈ AcceptableTermination
        return Subgradient(subproblem, x)
    elseif status == MOI.INFEASIBLE
        return Infeasible(subproblem)
    elseif status == MOI.DUAL_INFEASIBLE
        return Unbounded(subproblem)
    else
        error("Subproblem $(subproblem.id) was not solved properly, returned status code: $status")
    end
end

# Smooth #
# ========================== #
@with_kw mutable struct SmoothingData{T <: AbstractFloat}
    objective::AffineDecisionFunction{T} = zero(AffineDecisionFunction{T})
end

@with_kw mutable struct SmoothingParameters{T <: AbstractFloat}
    μ::T = 1.0
    objective_correction::Bool = false
end

"""
    SmoothSubProblem

Subproblem smoothed using a Moreau envelope.

...
# Parameters
- `μ::AbstractFloat = 1.0`: Moreau smoothing parameter. Controls the smoothing approximation accuracy.
...
"""
struct SmoothSubProblem{T <: AbstractFloat} <: AbstractSubProblem{T}
    data::SmoothingData{T}
    parameters::SmoothingParameters{T}

    id::Int
    probability::T
    model::JuMP.Model
    optimizer::MOI.AbstractOptimizer

    projection_targets::Vector{MOI.VariableIndex}
    penaltyterm::Quadratic

    function SmoothSubProblem(model::JuMP.Model,
                              id::Integer,
                              π::AbstractFloat,
                              ξ::AbstractVector;
                              kw...)
        T = typeof(π)
        # Get optimizer backend
        optimizer = backend(model)
        # Initialize data and params
        data = SmoothingData{T}()
        params = SmoothingParameters{T}(; kw...)
        @unpack μ = params
        # Cache objective
        F = MOI.get(backend(model), MOI.ObjectiveFunctionType())
        data.objective = MOI.get(backend(model), MOI.ObjectiveFunction{F}())
        # Create and return smoothed subproblem
        subproblem = new{T}(data,
                            params,
                            id,
                            π,
                            model,
                            optimizer,
                            Vector{MOI.VariableIndex}(),
                            Quadratic())
        # Setup Moreau envelope
        add_projection_targets!(subproblem, ξ)
        return subproblem
    end
end

function (subproblem::SmoothSubProblem)(x::AbstractVector)
    return solve_subproblem(subproblem, x)
end

function solve_subproblem(subproblem::SmoothSubProblem, x::AbstractVector)
    @unpack μ, objective_correction = subproblem.parameters
    decisions = index.(all_decision_variables(subproblem.model, 1))
    if objective_correction
        unfix.(all_decision_variables(subproblem.model, 1))
        StochasticPrograms.enable_penalty!(subproblem.penaltyterm,
                                           subproblem.optimizer,
                                           1 / (2 * μ),
                                           decisions,
                                           subproblem.projection_targets)
    end
    MOI.optimize!(subproblem.optimizer)
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if status ∈ AcceptableTermination
        return Gradient(subproblem, x)
    elseif status == MOI.INFEASIBLE
        return Infeasible(subproblem)
    elseif status == MOI.DUAL_INFEASIBLE
        return Unbounded(subproblem)
    else
        error("Subproblem $(subproblem.id) was not solved properly, returned status code: $status")
    end
end

function add_projection_targets!(subproblem::SmoothSubProblem{T}, ξ::AbstractVector) where T <: AbstractFloat
    @unpack μ = subproblem.parameters
    n = num_decisions(subproblem.model, 1)
    resize!(subproblem.projection_targets, n)
    for i in eachindex(ξ)
        name = add_subscript(:ξ, i)
        set = SingleDecisionSet(1, KnownDecision(ξ[i], T), NoSpecifiedConstraint(), false)
        empty = ScalarVariable(VariableInfo(false, NaN, false, NaN, false, NaN, false, NaN, false, false))
        dref = add_variable(subproblem.model, VariableConstrainedOnCreation(empty, set), name)
        subproblem.projection_targets[i] = index(dref)
    end
    # Initialize quadratic penalty
    x = index.(all_decision_variables(subproblem.model, 1))
    initialize_penaltyterm!(subproblem.penaltyterm,
                            subproblem.optimizer,
                            1 / (2 * μ),
                            x,
                            subproblem.projection_targets)
    return nothing
end

function restore_subproblem!(subproblem::SmoothSubProblem)
    MOI.optimize!(subproblem.optimizer)
    # Delete penalty-term
    remove_penalty!(subproblem.penaltyterm, subproblem.optimizer)
    # Delete projection targets
    for var in subproblem.projection_targets
        MOI.delete(subproblem.optimizer, var)
    end
    empty!(subproblem.projection_targets)
    return nothing
end

# Gradient #
# ========================== #
function Subgradient(subproblem::SubProblem{T}, x::AbstractVector) where T <: AbstractFloat
    π = subproblem.probability
    nterms = if isempty(subproblem.masterterms)
        nterms = 0
    else
        nterms = mapreduce(+, subproblem.masterterms) do terms
            length(terms)
        end
    end
    cols = zeros(Int, nterms)
    vals = zeros(T, nterms)
    j = 1
    for (i, ci) in enumerate(subproblem.linking_constraints)
        λ = MOI.get(subproblem.optimizer, MOI.ConstraintDual(), ci)
        for (row, col, coeff) in subproblem.masterterms[i]
            cols[j] = col
            vals[j] = π * λ[row] * coeff
            j += 1
        end
    end
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    # Create sense-corrected subgradient
    δQ = sparsevec(cols, vals, length(x))
    Q = correction * π * MOI.get(subproblem.optimizer, MOI.ObjectiveValue())
    return Gradient(δQ, Q, subproblem.id)
end

function Gradient(subproblem::SmoothSubProblem{T}, x::AbstractVector) where T <: AbstractFloat
    @unpack μ, objective_correction = subproblem.parameters
    π = subproblem.probability
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    # Get solution
    u = value.(all_decision_variables(subproblem.model, 1))
    # Sense corrected gradient
    δQ = π * (1/μ) * (x - u)
    # Sense corrected objective
    decisions = index.(all_decision_variables(subproblem.model, 1))
    if objective_correction
        fix.(all_decision_variables(subproblem.model, 1), x)
        StochasticPrograms.disable_penalty!(subproblem.penaltyterm,
                         subproblem.optimizer,
                         decisions,
                         subproblem.projection_targets)
        MOI.optimize!(subproblem.optimizer)
        status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
        if !(status ∈ AcceptableTermination)
            if status == MOI.INFEASIBLE
                return Infeasible(subproblem)
            elseif status == MOI.DUAL_INFEASIBLE
                return Unbounded(subproblem)
            end
        end
    end
    Q = correction * π * MOIU.eval_variables(subproblem.data.objective) do vi
        MOI.get(subproblem.optimizer, MOI.VariablePrimal(), vi)
    end
    return Gradient(δQ, Q, subproblem.id)
end

function Infeasible(subproblem::AbstractSubProblem)
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    return Gradient(sparsevec(Float64[]), correction * Inf, subproblem.id)
end

function Unbounded(subproblem::AbstractSubProblem)
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    return Gradient(sparsevec(Float64[]), correction * -Inf, subproblem.id)
end

# API
# ------------------------------------------------------------
"""
    Unaltered

Factory object for using regular [`SubProblem`](@ref) in the quasi-gradient algorithm. Passed by default to `subproblems` to `QuasiGradient.Optimizer`.

"""
struct Unaltered <: AbstractSubProblemState end
"""
    Smoothed

Factory object for using [`SmoothSubProblem`](@ref) through Moreau envelopes in the quasi-gradient algorithm. Pass to `subproblems ` in `QuasiGradient.Optimizer` or by setting the [`SubProblems`](@ref) attribute. See [`SmoothSubProblem`](@ref) for parameter descriptions.

"""
struct Smoothed <: AbstractSubProblemState
    parameters::SmoothingParameters{Float64}
end
Smoothed(; kw...) = Smoothed(SmoothingParameters(; kw...))


text:
using Powseeker
using POMDPs
using POMDPToolbox
using Plots

mdp = PowseekerMDP()
policy = Topout(mdp, 0.05)

hr = HistoryRecorder(max_steps=100)
# is = SkierState(1, [1200, 0], 0)
# is = SkierState(1, [-750, 1100], 0)
is = initial_state(mdp, Base.GLOBAL_RNG)
hist = simulate(hr, mdp, policy, is)
@show sum(reward_hist(hist))

plotlyjs()
plot(mdp)
plot!(hist)
gui()


text:
using MMA
using Base.Test

# write your own tests here
@test variance([1, 1, 1]) == 0
@test variance([1, 2, 4, 5]) == 2.5
@test_approx_eq(std_dev([1, 1, 1]), 0)
@test_approx_eq(std_dev([1, 3]), 1)
@test_approx_eq(dist([3, 0], [0, 4]), 5)
@test_approx_eq(rms([1, 1, 3, 5]), 3)
@test_approx_eq(MMA.angle([1, 0], [0, 1]), pi/2)


text:
module Shapes

import Base: similar, max, min, union!, empty!, isempty, convert
import LinearAlgebra: norm, dot, cross

export Shape, Empty, Space, Line, Plane, Sphere, AABB, Convex
export isvalid, transform, getnormal, getpoint, volume, union!, setplane, getintersection, intersect, outside, similar, min, max, inside, assign

iszero(x::T) where T = abs(x) < eps(T)

len2(x, y, z) = x*x + y*y + z*z
len2(x) = sum(abs2, x)
len(x, y, z) = sqrt(len2(x, y, z))
len(x) = sqrt(len2(x))

lerp(x, y, t) = x + (y-x)*t

function scale3(m::AbstractArray{T, 2}) where T
	l2 = max(len2(m[1,1], m[1,2], m[1,3]), len2(m[2,1], m[2,2], m[2,3]), len2(m[3,1], m[3,2], m[3,3]))
	sqrt(l2)
end

function transform_p3(result::AbstractArray{T, 1}, m::AbstractArray{T, 2}, p::AbstractArray{T, 1}) where T
	@simd for i = 1:3
		@inbounds result[i] = m[i,1]*p[1] + m[i,2]*p[2] + m[i,3]*p[3] + m[i,4]
	end
	nothing
end

function transform_p3(result::AbstractArray{T, 2}, m::AbstractArray{T, 2}, p::AbstractArray{T, 2}) where T
	for c = 1:size(p, 2)
		@simd for r = 1:3
			@inbounds result[r, c] = m[r,1]*p[1,c] + m[r,2]*p[2,c] + m[r,3]*p[3,c] + m[r,4]
		end
	end
	nothing
end

function transform_p4(result::AbstractArray{T, 1}, m::AbstractArray{T, 2}, p::AbstractArray{T, 1}) where T
	@simd for i = 1:3
		@inbounds result[i] = m[i,1]*p[1] + m[i,2]*p[2] + m[i,3]*p[3] + m[i,4]*p[4]
	end
	nothing
end

function transform_p4(result::AbstractArray{T, 2}, m::AbstractArray{T, 2}, p::AbstractArray{T, 2}) where T
	for c = 1:size(p, 2)
		@simd for r = 1:4
			@inbounds result[r, c] = m[r,1]*p[1,c] + m[r,2]*p[2,c] + m[r,3]*p[3,c] + m[r,4]*p[4,c]
		end
	end
	nothing
end

make_interval(a, b) = a < b ? (a, b) : (b, a)
empty_interval(i) = !(i[1] <= i[2]) # Will return true in case one of the values is NaN
infinite_interval(i) = isinf(i[1]) || isinf(i[2])
intersect_interval(min1, max1, min2, max2) = (max(min1, min2), min(max1, max2))
intersect_interval(i1, i2) = intersect_interval(i1[1], i1[2], i2[1], i2[2])
is_intersecting_interval(i1, i2) = i1[2] >= i2[1] && i1[1] <= i2[2] # will return false if either interval contains a NaN

# roots of ax^2+bx+c=0
function quadroots(a, b, c)
	d = b*b-4a*c
	if d < zero(d)
		return oftype(a, NaN), oftype(a, NaN)
	end
	sd = sqrt(d)
	x1 = (-b-sd)/2a
	x2 = (-b+sd)/2a
	return x1, x2
end

planevalue(p::Vector{T}, point::Vector{T}) where T= dot(p[1:3], point) + p[4]
planenormal(p::Vector{T}) where T = p[1:3]
function planepoint(p::Vector{T}) where T
	i = indmax(abs(p[i]) for i=1:3)
	return T[x==i ? -p[4]/p[i] : 0 for x=1:3]
end

function normalized_plane(dst, col, n, p)
	invLen = 1 / norm(n)
	dst[1, col] = n[1] * invLen
	dst[2, col] = n[2] * invLen
	dst[3, col] = n[3] * invLen
	dst[4, col] = -dot(p, n) * invLen
	nothing
end

function set_aabb_planes(dst::Matrix{T}, ptMin::Vector{T}, ptMax::Vector{T}) where T
	normalized_plane(dst, 1, T[ 1,  0,  0], ptMin)
	normalized_plane(dst, 2, T[ 0,  1,  0], ptMin)
	normalized_plane(dst, 3, T[ 0,  0,  1], ptMin)
	normalized_plane(dst, 4, T[-1,  0,  0], ptMax)
	normalized_plane(dst, 5, T[ 0, -1,  0], ptMax)
	normalized_plane(dst, 6, T[ 0,  0, -1], ptMax)
end

function normalize_plane(dst, col)
	invLen = 1 / len(dst[1, col], dst[2, col], dst[3, col])
	dst[1, col] *= invLen
	dst[2, col] *= invLen
	dst[3, col] *= invLen
	dst[4, col] *= invLen
	nothing
end

function normalize_planes(planes::Matrix{T}) where T
	for i = 1:size(planes, 2)
		normalize_plane(planes, i)
	end
end

abstract type Shape{T <: Real} end

similar(s::S) where S <: Shape = S()
function transform(m::Matrix, s::S) where S <: Shape
	dst = similar(s)
	transform(dst, m, s)
	dst
end

# empty shape
mutable struct Empty{T} <: Shape{T}
end

isvalid(e::Empty) = true
volume(e::Empty{T}) where T = zero(T)
transform(eDest::Empty, m::Matrix, e::Empty) = nothing
#union!(e::Empty{T}, s::Shape{T}) where T = s

# all space
mutable struct Space{T} <: Shape{T}
end

isvalid(s::Space) = true
volume(s::Space{T}) where T = T(Inf)
transform(sDest::Space, m::Matrix, s::Space) = nothing

mutable struct Line{T} <: Shape{T}
	p::Array{T, 2}

	Line{T}(p::Array{T, 2}) where T = new(p)
	Line{T}() where T = new(Array{T}(undef, 3, 2))
end

function Line(x1::T, y1, z1, x2, y2, z2) where T
	p = Array{T}(undef, 3, 2)
	p[1,1] = x1
	p[2,1] = y1
	p[3,1] = z1
	p[1,2] = x2
	p[2,2] = y2
	p[3,2] = z2
	Line{T}(p)
end

Line(p0, p1) = Line(p0..., p1...)
Line(v1::Vector{T}, v2::Vector{T}) where T = Line(v1[1], v1[2], v1[3], v2[1], v2[2], v2[3])

isvalid(l::Line{T}) where T = size(l.p)==(3, 2) && len2(l.p[1, 2] - l.p[1, 1], l.p[2, 2] - l.p[2, 1], l.p[3, 2] - l.p[3, 1]) >= eps(T)

function getvector(l::Line{T}) where T
	v = Vector{T}(3)
	v[1] = l.p[1, 2] - l.p[1, 1]
	v[2] = l.p[2, 2] - l.p[2, 1]
	v[3] = l.p[3, 2] - l.p[3, 1]
	v
end

getpoint(l::Line{T}, t::T) where T = lerp(l.p[:, 1], l.p[:, 2], t)

function transform(lDest::Line{T}, m::Matrix{T}, l::Line{T}) where T
	for c = 1:2
		@simd for r = 1:3
			@inbounds lDest.p[r, c] = m[r,1]*l.p[1,c] + m[r,2]*l.p[2,c] + m[r,3]*l.p[3,c] + m[r,4]
		end
	end
	nothing
end

function plane2interval(p::Vector{T}, l::Line{T}, interval::Tuple{T, T} = (T(-Inf), T(Inf))) where T
	@assert isvalid(l)
	@assert size(p) == (4,)

	int = line2plane(l.p, p)
	if !isinf(int)
		# line doesn't lie on the plane
		if isnan(int)
			# line is parallel to the plane
			if planevalue(p, l.p[:, 1]) < zero(T)
				# and on the negative side
				interval = (T(Inf), T(-Inf))
			end
		else
			rayInterval = dot(planenormal(p), getvector(l)) < 0 ? (T(-Inf), int) : (int, T(Inf))
			interval = intersect_interval(interval, rayInterval)
		end
	end
	interval
end

mutable struct Plane{T} <: Shape{T}
	p::Array{T, 1}

	Plane{T}(p::Array{T, 1}) where T = new(p)
	Plane{T}() where T = new(Array{T}(undef, 4))
end

Plane(a::T, b, c, d) where T = Plane{T}(T[a, b, c, d])
Plane(n, p) = Plane(n..., -dot(p, n))

isvalid(p::Plane{T}) where T = size(p.p) == (4,) && len2(p.p[1:3]) >= eps(T)

getnormal(p::Plane{T}) where T = planenormal(p.p)
getvalue(p::Plane{T}, point::Vector{T}) where T = planevalue(p.p, point)
getpoint(p::Plane{T}) where T = planepoint(p.p)

# multiply plane by inverse transpose of matrix
function transform(pDest::Plane{T}, m::Matrix{T}, p::Plane{T}) where T
	m_it = transpose(inv(m))
	transform_p4(pDest.p, m_it, p.p)
	nothing
end

# transform with inverse transpose pre-computed
transform_it(pDest::Plane{T}, m_it::Matrix{T}, p::Plane{T}) where T = transform_p4(pDest.p, m_it, p.p)

mutable struct Sphere{T} <: Shape{T}
	c::Array{T, 1}
	r::T

	Sphere{T}(c::Array{T, 1}, r) where T = new(c, r)
	Sphere{T}() where T = new(Array{T}(undef, 3), 0)
end

Sphere(cx::T, cy, cz, r) where T = Sphere{T}(T[cx, cy, cz], convert(T, r))
Sphere(c, z) = Sphere(c..., z)

function assign(s1::Sphere, s2::Sphere)
	s1.c[:] = s2.c
	s1.r = s2.r
	s1
end

isvalid(s::Sphere{T}) where T = size(s.c) == (3,) && !isempty(s)
volume(s::Sphere) = s.r^3*4pi/3
isempty(s::Sphere{T}) where T = s.r < zero(T)

function empty!(s::Sphere{T}) where T
	s.r = T(-Inf)
	s
end

function union!(s::Sphere{T}, p::Vector{T}) where T
	@assert isvalid(s)
	@assert length(p) == 3

	d = len(p-s.c)
	if d > s.r
		s.r = (d+s.r)/2
		s.c = lerp(p, s.c, s.r/d)
	end
	s
end

function union!(s1::Sphere{T}, s2::Sphere{T}) where T
	c12 = s2.c - s1.c
	d = norm(c12)
	if s1.r >= d + s2.r
		# s1 is the union, nothing needs changing
	elseif s2.r >= d + s1.r
		# s2 is the union
		s1.c[:] = s2.c
		s1.r = s2.r
	else
		d1 = (d + s2.r - s1.r) / 2
		s1.c += c12 * (d1 / d)
		s1.r += d1
	end
	s1
end

function transform(sDest::Sphere{T}, m::Matrix{T}, s::Sphere{T}) where T
	transform_p3(sDest.c, m, s.c)
	sDest.r = s.r * scale3(m)
	nothing
end

mutable struct AABB{T} <: Shape{T}
	p::Array{T, 2}

	AABB{T}(p::Array{T, 2}) where T = new(p)
	AABB{T}() where T = new(Array{T}(undef, 3, 2))
end

AABB(xmin::T, ymin, zmin, xmax, ymax, zmax) where T = AABB{T}(T[xmin xmax; ymin ymax; zmin zmax])
AABB(pmin, pmax) = AABB(pmin..., pmax...)

function assign(ab1::AABB, ab2::AABB)
	ab1.p[:] = ab2.p
	ab1
end

isvalid(aabb::AABB{T}) where T = size(aabb.p) == (3, 2) && !isempty(aabb)
volume(ab::AABB) = prod(ab.p[i, 2] - ab.p[i, 1] for i=1:3)
isempty(ab::AABB{T}) where T = ab.p[1, 1] > ab.p[1, 2] || ab.p[2, 1] > ab.p[2, 2] || ab.p[3, 1] > ab.p[3, 2]

function empty!(ab::AABB{T}) where T
	ab.p[:, 1] .= T(Inf)
	ab.p[:, 2] .= T(-Inf)
	ab
end

min(ab::AABB{T}) where T = ab.p[:, 1]
max(ab::AABB{T}) where T = ab.p[:, 2]

function addpoint!(minMax::Matrix{T}, p::Vector{T}) where T
	@assert size(minMax) == (3, 2)
	@assert length(p) == 3

	for i = 1:3
		if p[i] < minMax[i, 1]
			minMax[i, 1] = p[i]
		end
		if minMax[i, 2] < p[i]
			minMax[i, 2] = p[i]
		end
	end
	nothing
end

union!(ab::AABB{T}, p::Vector{T}) where T = addpoint!(ab.p, p)
function union!(ab1::AABB{T}, ab2::AABB{T}) where T
	for i = 1:3
		ab1.p[i, 1] = min(ab1.p[i, 1], ab2.p[i, 1])
		ab1.p[i, 2] = max(ab1.p[i, 2], ab2.p[i, 2])
	end
	ab1
end

function transform(abDest::AABB{T}, m::Matrix{T}, ab::AABB{T}) where T
	t = Array{T}(undef, 3)
	p = Array{T}(undef, 3)
	p[1] = ab.p[1, 1]
	p[2] = ab.p[2, 1]
	p[3] = ab.p[3, 1]
	transform_p3(t, m, p)
	abDest.p[1, 1] = t[1]
	abDest.p[2, 1] = t[2]
	abDest.p[3, 1] = t[3]
	abDest.p[1, 2] = t[1]
	abDest.p[2, 2] = t[2]
	abDest.p[3, 2] = t[3]
	for i = 1:7
		x = i & 1
		y = (i & 2)>>1
		z = (i & 4)>>2
		p[1] = ab.p[1, x+1]
		p[2] = ab.p[2, y+1]
		p[3] = ab.p[3, z+1]
		transform_p3(t, m, p)
		addpoint!(abDest.p, t)
	end
	nothing
end

# todo: add Capsule and OBB

mutable struct EdgeInfo{T <: Real}
	line::Line{T}
	interval::Tuple{T, T}
end

mutable struct Adjacency{T <: Real}
	faces::Vector{Vector{Int}}
	edges::Dict{Tuple{Int, Int}, EdgeInfo{T}}
	Adjacency{T}() where {T <: Real} = new()
end

isvalid(a::Adjacency) = isdefined(a, :faces) && isdefined(a, :edges)

function init(adj::Adjacency{T}, planes::Matrix{T}) where T
	adj.faces = [Int[] for i=1:size(planes, 2)]
	adj.edges = Dict{Tuple{Int, Int}, EdgeInfo{T}}()
	calcedges(adj, planes)
end

function init(adj::Adjacency{T}, ab::AABB{T}) where T
	adj.faces = [filter(j->i != j != (i+3-1)%6+1, 1:6) for i=1:6]
	adj.edges = Dict{Tuple{Int, Int}, EdgeInfo{T}}()

	sz = ab.p[:, 2] - ab.p[:, 1]
	for i = 1:6, j in adj.faces[i]
		iDir = (i-1)%3+1
		if i < j
			jDir = (j-1)%3+1
			pts = Array{T}(undef, 3, 2)
			for k = 1:3
				m = 1
				d = 0
				if k == iDir
					m = (i > 3) + 1
				elseif k == jDir
					m = (j > 3) + 1
				else
					d = sz[k]
				end
				pts[k, 1] = ab.p[k, m]
				pts[k, 2] = pts[k, 1] + d
			end
			l = Line{T}(pts)
			int = (T(0), T(1))
			adj.edges[(i, j)] = EdgeInfo(l, int)
		end
	end
end

function addplane(adj::Adjacency{T}, planes::Matrix{T}, index::Int, planeRange::UnitRange{Int} = 1:index-1) where T
	@assert size(planes, 1) == 4

	p = Plane{T}(planes[:, index])
	update_edges(adj, p, planeRange)
	for i = planeRange
		pi = Plane{T}(planes[:, i])
		intersection = getintersection(p, pi)
		if isa(intersection, Line)
			interval = (T(-Inf), T(Inf))
			for k = planeRange
				if k != i
					interval = plane2interval(planes[:, k], intersection, interval)
					empty_interval(interval) && break
				end
			end
			if !empty_interval(interval)
				f1 = min(i, index)
				f2 = max(i, index)
				adj.edges[(f1, f2)] = EdgeInfo(intersection, interval)
				push!(adj.faces[f1], f2)
				push!(adj.faces[f2], f1)
			end
		end
	end
end

function update_edges(adj::Adjacency{T}, p::Plane{T}, planeRange::UnitRange{Int}) where T
	for i in planeRange
		k = 1
		while k < length(adj.faces[i])
			j = adj.faces[i][k]
			if j > i && j in planeRange
				e = adj.edges[(i, j)]
				interval = plane2interval(p.p, e.line, e.interval)
				if empty_interval(interval)
					deleteat!(adj.faces[i], k)
					deleteat!(adj.faces[j], findfirst(isequal(i), adj.faces[j]))
					delete!(adj.edges, (i, j))
					continue
				else
					e.interval = interval
				end
			end
			k += 1
		end
	end
end

function combine(dst::Adjacency{T}, src1::Adjacency{T}, src2::Adjacency{T}, planes::Matrix{T}) where T
	@assert isvalid(src1)
	@assert isvalid(src2)
	@assert size(planes, 2) == length(src1.faces) + length(src2.faces)

	dst.faces = deepcopy(src1.faces)
	dst.edges = deepcopy(src2.edges)

	src2Base = length(src1.faces)
	for f in src2.faces
		push!(dst.faces, map(i->i+src2Base, f))
	end
	for (f, e) in src2.edges
		dst.edges[(f[1] + src2Base, f[2] + src2Base)] = e
	end

	calcedges(dst, planes, src2Base)
end

function calcedges(adj::Adjacency{T}, planes::Matrix{T}, splitLength::Int = 0) where T
	@assert size(planes, 1) == 4

	cols = size(planes, 2)
	range1 = 1:splitLength
	range2 = splitLength+1:cols
	for i in range1
		update_edges(adj, Plane{T}(Planes[:, i]), range2)
	end
	for i in range2
		if splitLength == 0
			range1 = 1:i-1
		end
		addplane(adj, planes, i, range1)
	end

	#=
	iRange = 1:(splitLength == 0? cols : splitLength)
	for i = iRange
		pi = Plane{T}(planes[:, i])
		jRange = (splitLength == 0? i+1 : splitLength+1):cols
		for j = i+1:cols
			pj = Plane{T}(planes[:, j])
			intersection = getintersection(pi, pj)
			if isa(intersection, Line)
				interval = (T(-Inf), T(Inf))
				for k = 1:cols
					if k != i && k != j
						interval = plane2interval(planes[:, k], intersection, interval)
						empty_interval(interval) && break
					end
				end
				if !empty_interval(interval)
					adj.edges[(i, j)] = EdgeInfo(intersection, interval)
					push!(adj.faces[i], j)
					push!(adj.faces[j], i)
				end
			end
		end
	end
	=#
	adj
end

function select_planes(adj::Adjacency{T}, keep::Vector{Bool}) where T
	@assert length(adj.faces) == length(keep)
	ind = Vector{T}(length(keep))
	nextInd = 1
	for i = 1:length(ind)
		if keep[i]
			ind[i] = nextInd
			nextInd += 1
		else
			ind[i] = 0
		end
	end
	adj.faces = adj.faces[keep]
	for i = 1:length(adj.faces)
		filter!(f->keep[f], adj.faces[i])
		map!(f->ind[f], adj.faces[i])
	end
	edges = adj.edges
	adj.edges = Dict{Tuple{Int, Int}, EdgeInfo{T}}()
	for (f, e) in edges
		if keep[f[1]] && keep[f[2]]
			adj.edges[(ind[f[1]], ind[f[2]])] = e
		end
	end
end

function transform(dst::Adjacency{T}, m::Matrix{T}, src::Adjacency{T}) where T
	dst.faces = deepcopy(src.faces)
	if !isdefined(dst.edges)
		dst.edges = Dict{Tuple{Int, Int}, EdgeInfo{T}}()
	else
		empty!(dst.edges)
	end
	for (f, e) in src.edges
		lDst = Line{T}()
		transform(lDst, m, e.line)
		dst.edges[f] = EdgeInfo(lDst, e.interval)
	end
	nothing
end

mutable struct Convex{T} <: Shape{T}
	planes::Array{T, 2}
	adj::Adjacency{T}

	Convex{T}(planes::Matrix{T}) where T= new(planes, Adjacency{T}())
end

Convex(::Type{T}, planeCount::Int) where T = Convex{T}(zeros(T, 4, planeCount))
Convex(c1::Convex{T}, c2::Convex{T}) where T = combine(Convex(T, size(c1.planes, 2) + size(c2.planes, 2)), c1, c2)
Convex(ab::AABB{T}) where T = set(Convex{T}(Array{T}(undef, 4, 6)), ab)

function Convex(planes::Plane{T}...) where T
	c = Convex(T, length(planes))
	for i in 1:length(planes)
		setplane(c, i, planes[i])
	end
	c
end

isvalid(c::Convex{T}) where T = size(c.planes, 1) == 4 && size(c.planes, 2) > 0
setplane(c::Convex{T}, planeIndex::Int, p::Vector{T}) where T = c.planes[:, planeIndex] = p / len(p[1:3])
isclosed(c::Convex{T}) where T = !(isempty(c.edges) || any(e->infinite_interval(e.interval), values(c.edges)))

function isempty(c::Convex{T}) where T
	@assert(isvalid(c))
	@assert(isvalid(c.adj))
	!isempty(c.planes) && isempty(c.adj.edges)
end

function set(c::Convex{T}, ab::AABB{T}) where T
	@assert size(c.planes) == (4, 6)
	set_aabb_planes(c.planes, ab.p[:, 1], ab.p[:, 2])
	init(c.adj, ab)
	c
end

function combine(dst::Convex{T}, src1::Convex{T}, src2::Convex{T}) where T
	src1Size = size(src1.planes, 2)
	src2Size = size(src2.planes, 2)
	if size(dst.planes, 2) != src1Size + src2Size
		dsp.planes = Array{T}(undef, 4, src1Size + src2Size)
	end
	dst.planes[:, 1:src1Size] = src1.planes
	dst.planes[:, src1Size+1:src1Size+src2Size] = src2.planes
	combine(dst.adj, src1.adj, src2.adj, dst.planes)
	dst
end

function calcedges(c::Convex{T}) where T
	@assert isvalid(c)
	normalize_planes(c.planes)
	init(c.adj, c.planes)
end

function select_planes(c::Convex{T}, keep::Vector{Bool}) where T
	c.planes = c.planes[:, keep]
	select_planes(c.adj, keep)
end

function remove_redundant(c::Convex{T}) where T
	@assert isvalid(c)
	@assert isvalid(c.adj)

	cols = size(c.planes, 2)
	keep = [true for i=1:cols]
	if isempty(c.adj.edges)
		for i = 1:cols, j = i+1:cols
			if keep[i] && keep[j] && iszero(dot(planenormal(c.planes[:, i]), planenormal(c.planes[:, j])) - 1)
				# planes are parallel, one of them is redundant
				jInFront = planevalue(c.planes[:, j], planepoint(c.planes[:, i])) < 0
				keep[i] = !jInFront
				keep[j] = jInFront
			end
		end
	else
		for i = 1:cols
			# plane does not form a common edge with any plane if it's redundant
			keep[i] = !isempty(c.adj.faces[i])
		end
	end
	if !all(keep)
		# we're actually removing planes
		select_planes(adj, keep)
	end
end

function transform(cDest::Convex{T}, m::Matrix{T}, c::Convex{T}) where T
	transform_it(cDest, transpose(inv(m)), c)
	if isvalid(c.adj)
		transform(cDest.adj, m, c.adj)
	end
	nothing
end

function transform_it(cDest::Convex{T}, m_it::Matrix{T}, c::Convex{T}) where T
	cols = size(c.planes, 2)
	if size(cDest.planes, 2) != cols
		cDest.planes = Array{T}(undef, 3, cols)
	end
	cDest.planes = m_it * c.planes
	for i = 1:cols
		invL = 1/len(cDest.planes[1, i], cDest.planes[2, i], cDest.planes[3, i])
		cDest.planes[1, i] *= invL
		cDest.planes[2, i] *= invL
		cDest.planes[3, i] *= invL
		cDest.planes[4, i] *= invL
	end
	nothing
end

# Conversion functions

convert(::Type{AABB{T}}, s::Sphere{T}) where T = AABB(s.c-s.r, s.c+s.r)
convert(::Type{AABB{T}}, e::Empty{T}) where T = empty!(AABB{T}())
convert(::Type{Sphere{T}}, ab::AABB{T}) where T = Sphere{T}(0.5(min(ab) + max(ab)), 0.5len(max(ab) - min(ab)))
convert(::Type{Sphere{T}}, e::Empty{T}) where T = empty!(Sphere{T}())

# Intersection functions

getintersection(s1::Shape{T}, s2::Shape{T}) where T = getintersection(s2, s1)
intersect(s1::Shape{T}, s2::Shape{T}) where T = intersect(s2, s1)

function line2plane(l::Matrix{T}, plane::Vector{T}) where T
	lv = Vector{T}(3)
	lv[1] = l[1, 2] - l[1, 1]
	lv[2] = l[2, 2] - l[2, 1]
	lv[3] = l[3, 2] - l[3, 1]
	pn = planenormal(plane)
	po = planepoint(plane)
	pl = Vector{T}(3)
	pl[1] = po[1] - l[1, 1]
	pl[2] = po[2] - l[2, 1]
	pl[3] = po[3] - l[3, 1]

	vn = dot(lv, pn)
	d = dot(pl, pn)
	if iszero(vn)
		# line and plane are parallel
		return iszero(d) ? T(Inf) : T(NaN)
	end

	return d / vn
end

function getintersection(l::Line{T}, p::Plane{T}) where T
	@assert isvalid(l)
	@assert isvalid(p)
	line2plane(l.p, p.p)
end

intersect(l::Line{T}, p::Plane{T}) where T = !isnan(getintersection(l, p))

function getintersection(l::Line{T}, s::Sphere{T}) where T
	@assert isvalid(l)
	@assert isvalid(s)

	po = l.p[:, 1]
	vo = l.p[:, 2] - po
	co = po - s.c
	return quadroots(dot(vo, vo), 2dot(vo, co), dot(co, co)-s.r*s.r)
end

intersect(l::Line{T}, s::Sphere{T}) where T = !empty_interval(getintersection(l, s))

function getintersection(l::Line{T}, ab::AABB{T}) where T
	@assert isvalid(l)
	@assert isvalid(ab)

	tInt = (T(-Inf), T(Inf))
	for i = 1:3
		lo = l.p[i, 1]
		d = l.p[i, 2] - lo
		if iszero(d)
			# line is parallel to this box's dimension
			if ab.p[i, 1] <= lo <= ab.p[i, 2]
				continue
			else
				tInt = (T(Inf), T(-Inf))
				break
			end
		end
		dimInt = make_interval((ab.p[i, 1] - lo)/d, (ab.p[i, 2] - lo)/d)
		tInt = intersect_interval(tInt, dimInt)
		empty_interval(tInt) && break
	end
	tInt
end

intersect(l::Line{T}, ab::AABB{T}) where T = !empty_interval(getintersection(l, ab))

function getintersection(ab1::AABB{T}, ab2::AABB{T}) where T
	result = AABB{T}()
	for i = 1:3
		result.p[i, 1] = max(ab1.p[i, 1], ab2.p[i, 1])
		result.p[i, 2] = min(ab1.p[i, 2], ab2.p[i, 2])
	end
	return result
end

intersect(ab1::AABB{T}, ab2::AABB{T}) where T = isvalid(getintersection(ab1, ab2))

# returns a line, a plane or nothing depending on the relative position of the planes
function getintersection(p1::Plane{T}, p2::Plane{T}) where T
	@assert isvalid(p1)
	@assert isvalid(p2)

	n1 = getnormal(p1)
	n2 = getnormal(p2)
	v = cross(n1, n2)
	anypt2 = getpoint(p2)
	if len2(v) < eps(T)
		#planes are parallel
		if iszero(getvalue(p1, anypt2))
			# a point on one plane lies on the other, so they are coincident
			return p1
		else
			# no intersection
			return nothing
		end
	end
	axis2 = cross(v, n2)
	lineOnP2 = Line(anypt2, anypt2+axis2)
	t = getintersection(lineOnP2, p1)
	@assert !isnan(t) && !isinf(t)
	ptIntersection = getpoint(lineOnP2, t)
	return Line(ptIntersection, ptIntersection + v)
end

intersect(p1::Plane, p2::Plane) = getintersection(p1, p2) != nothing

intersect(c::Convex, e::Empty) = false
intersect(c::Convex, s::Space) = isvalid(c)

function getintersection(c::Convex{T}, l::Line{T}, interval::Tuple{T, T} = (T(-Inf), T(Inf))) where T
	@assert isvalid(l)
	@assert isvalid(c)

	for i = 1:size(c.planes, 2)
		interval = plane2interval(c.planes[:, i], l, interval)
		empty_interval(interval) && break
	end
	interval
end

intersect(c::Convex, l::Line) = !empty_interval(getintersection(c, l))

function getintersection(c::Convex{T}, p::Plane{T}) where T
	c = Convex{T}(hcat(c.planes, p.p))
	calcedges(c)
	remove_redundant(c)
	c
end

intersect(c::Convex, p::Plane) = !isempty(getintersection(c, p))

function intersect(c::Convex{T}, s::Sphere{T}) where T
	@assert isvalid(c)
	@assert isvalid(c.adj)
	@assert isvalid(s)

	cols = size(c.planes, 2)
	incidentPlanes = Int[]
	for i = 1:cols
		dist = planevalue(c.planes[:, 4], s.c)
		dist < -s.r && return false
		if dist <= s.r
			push!(incidentPlanes, i)
		end
	end
	# if the sphere isn't intersected by any of the planes, and is not wholly on the negative side of any of them, it's inside
	isempty(incidentPlanes) && return true
	for i in incidentPlanes
		plane = c.planes[:, i]
		pn = planenormal(plane)
		po = planepoint(plane)
		projCenter = s.c - pn * dot(s.c - po, pn) # sphere center projected on the plane
		@assert iszero(planevalue(plane, projCenter))
		centerIsInside = true
		for j in c.adj.faces[i]
			# planes with indices i and j have a common edge so they are incident, check if the projected center is on the positive side of j
			if planevalue(c.planes[j], projCenter) < 0
				centerIsInside = false
				break
			end
		end
		# the projection of the center on the plane is inside the area defined by the incident planes of the convex
		centerIsInside && return true
	end
	# the sphere is incident with some planes, but its center projected on any of them doesn't lie inside a convex face
	# in this case the sphere intersects the convex only if one of the convex edges intersects the sphere
	# and the only edges the sphere can intersect are those formed between pairs of incident planes
	for i = 1:length(incidentPlanes), j = i+1:length(incidentPlanes)
		planePair = (incidentPlanes[i], incidentPlanes[j])
		if haskey(c.adj.edges, planePair)
			edge = c.adj.edges[planePair]
			interval = getintersection(edge.line, s)
			# do the line's edge interval and sphere interval intersect?
			is_intersecting_interval(interval, edge.interval) && return true
		end
	end
	return false
end

getintersection(c::Convex{T}, ab::AABB{T}) where T = getintersection(c, Convex(ab))
intersect(c::Convex{T}, ab::AABB{T}) where T = !isempty(getintersection(c, ab))

function getintersection(c1::Convex{T}, c2::Convex{T}) where T
	c = Convex(c1, c2)
	remove_redundant(c)
	c
end

intersect(c1::Convex, c2::Convex) = !isempty(getintersection(c1, c2))

outside(c::Convex, e::Empty) = true
outside(c::Convex, s::Space) = false

function outside(c::Convex{T}, s::Sphere{T}) where T
	@assert isvalid(c)
	@assert isvalid(s)
	for i = 1:size(c.planes, 2)
		if dot(c.planes[1:3, i], s.c) + c.planes[4, i] < -s.r
			return true
		end
	end
	return false
end

function outside(c::Convex{T}, ab::AABB{T}) where T
	@assert isvalid(c)
	@assert isvalid(ab)

	for i = 1:size(c.planes, 2)
		n = c.planes[1:3, i]
		# select the point on the box with the maximum value when substituted in the plane equation
		minPt = [n[j]<0 ? ab.p[j, 1] : ab.p[j, 2] for j=1:3]
		if dot(n, minPt) + c.planes[4, i] < 0
			return true
		end
	end
	return false
end

function inside(enclosing::AABB{T}, ab::AABB{T}) where T
	@assert size(enclosing.p) == (3,2)
	@assert size(ab.p) == (3,2)

	@simd for i = 1:3
		@inbounds if enclosing.p[i,1] > ab.p[i,1] || enclosing.p[i,2] < ab.p[i,2]
			return false
		end
	end
	return true
end

function inside(enclosing::AABB{T}, p::Vector{T}) where T
	@assert size(enclosing.p) == (3,2)
	@assert size(p) == (3,)

	@simd for i = 1:3
		@inbounds if !(enclosing.p[i,1] <= p[i] <= enclosing.p[i,2])
			return false
		end
	end
	return true
end

function inside(enclosing::Sphere{T}, p::Vector{T}) where T
	@assert size(enclosing.c) == (3,)
	@assert size(p) == (3,)

	len(p - enclosing.c) <= enclosing.r
end

function inside(enclosing::Convex{T}, p::Vector{T}) where T
	@assert isvalid(enclosing)
	@assert size(p) == (3,)

	for i = 1:size(enclosing.planes, 2)
		planevalue(enclosing.planes[:, i], p) < 0 && return false
	end
	true
end

end


text:
import HDF5

const HDF5_TS_ROOT_PATH = "time_series"

"""
Stores all time series data in an HDF5 file.

The file used is assumed to be temporary and will be automatically deleted when there are
no more references to the storage object.
"""
mutable struct Hdf5TimeSeriesStorage <: TimeSeriesStorage
    file_path::String
    read_only::Bool
end

"""
Constructs Hdf5TimeSeriesStorage by creating a temp file.
"""
function Hdf5TimeSeriesStorage()
    return Hdf5TimeSeriesStorage(true)
end

"""
Constructs Hdf5TimeSeriesStorage.

# Arguments
- `create_file::Bool`: create new file
- `filename=nothing`: if nothing, create a temp file, else use this name.
- `directory=nothing`: if set and filename is nothing, create a temp file in this
   directory. Use tempdir() if not set. This should be set if the time series data is larger
   than the tmp filesystem can hold.
- `read_only = false`: If true, don't allow changes to the file. Allows simultaneous read
   access.
"""
function Hdf5TimeSeriesStorage(
    create_file::Bool;
    filename = nothing,
    directory = nothing,
    read_only = false,
)
    if create_file
        if isnothing(filename)
            if isnothing(directory)
                directory = tempdir()
            end
            filename, io = mktemp(directory)
            close(io)
        end

        storage = Hdf5TimeSeriesStorage(filename, read_only)
        _make_file(storage)
    else
        storage = Hdf5TimeSeriesStorage(filename, read_only)
    end

    @debug "Constructed new Hdf5TimeSeriesStorage" storage.file_path read_only

    return storage
end

"""
Constructs Hdf5TimeSeriesStorage from an existing file.
"""
function from_file(
    ::Type{Hdf5TimeSeriesStorage},
    filename::AbstractString;
    read_only = false,
    directory = nothing,
)
    if !isfile(filename)
        error("time series storage $filename does not exist")
    end

    if read_only
        file_path = abspath(filename)
    else
        parent = isnothing(directory) ? tempdir() : directory
        file_path, io = mktemp(parent)
        close(io)
        copy_file(filename, file_path)
    end
    storage = Hdf5TimeSeriesStorage(false; filename = file_path, read_only = read_only)
    @info "Loaded time series from storage file existing=$filename new=$(storage.file_path)"
    return storage
end

get_file_path(storage::Hdf5TimeSeriesStorage) = storage.file_path

function serialize_time_series!(
    storage::Hdf5TimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
    ts::TimeSeriesData,
)
    check_read_only(storage)
    uuid = string(get_uuid(ts))
    component_name = make_component_name(component_uuid, name)

    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        if !haskey(root, uuid)
            HDF5.create_group(root, uuid)
            path = root[uuid]
            data = get_array_for_hdf(ts)
            path["data"] = data
            _write_time_series_attributes!(storage, ts, path)
            path["components"] = [component_name]
            @debug "Create new time series entry." uuid component_uuid name
        else
            path = root[uuid]
            @debug "Add reference to existing time series entry." uuid component_uuid name
            _append_item!(path, "components", component_name)
        end
    end

    return
end

"""
Return a String for the data type of the forecast data, this implementation avoids the use of `eval` on arbitrary code stored in HDF dataset.
"""
function get_data_type(ts::TimeSeriesData)
    data_type = eltype_data(ts)
    if data_type <: CONSTANT
        return "CONSTANT"
    elseif data_type == POLYNOMIAL
        return "POLYNOMIAL"
    elseif data_type == PWL
        return "PWL"
    elseif data_type <: Integer
        # We currently don't convert integers stored in TimeSeries.TimeArrays to floats.
        # This is a workaround.
        return "CONSTANT"
    else
        error("$data_type is not supported in forecast data")
    end
end

function _write_time_series_attributes!(
    storage::Hdf5TimeSeriesStorage,
    ts::T,
    path,
) where {T <: StaticTimeSeries}
    _write_time_series_attributes_common!(storage, ts, path)
end

function _write_time_series_attributes!(
    storage::Hdf5TimeSeriesStorage,
    ts::T,
    path,
) where {T <: Forecast}
    _write_time_series_attributes_common!(storage, ts, path)
    interval = get_interval(ts)
    HDF5.attributes(path)["interval"] = time_period_conversion(interval).value
end

function _write_time_series_attributes_common!(storage::Hdf5TimeSeriesStorage, ts, path)
    initial_timestamp = Dates.datetime2epochms(get_initial_timestamp(ts))
    resolution = get_resolution(ts)
    data_type = get_data_type(ts)
    HDF5.attributes(path)["module"] = string(parentmodule(typeof(ts)))
    HDF5.attributes(path)["type"] = string(nameof(typeof(ts)))
    HDF5.attributes(path)["initial_timestamp"] = initial_timestamp
    HDF5.attributes(path)["resolution"] = time_period_conversion(resolution).value
    HDF5.attributes(path)["data_type"] = data_type
end

function _read_time_series_attributes(
    storage::Hdf5TimeSeriesStorage,
    path,
    rows,
    ::Type{T},
) where {T <: StaticTimeSeries}
    return _read_time_series_attributes_common(storage, path, rows)
end

function _read_time_series_attributes(
    storage::Hdf5TimeSeriesStorage,
    path,
    rows,
    ::Type{T},
) where {T <: Forecast}
    data = _read_time_series_attributes_common(storage, path, rows)
    data["interval"] = Dates.Millisecond(HDF5.read(HDF5.attributes(path)["interval"]))
    return data
end

const _TYPE_DICT = Dict("CONSTANT" => CONSTANT, "POLYNOMIAL" => POLYNOMIAL, "PWL" => PWL)

function _read_time_series_attributes_common(storage::Hdf5TimeSeriesStorage, path, rows)
    initial_timestamp =
        Dates.epochms2datetime(HDF5.read(HDF5.attributes(path)["initial_timestamp"]))
    resolution = Dates.Millisecond(HDF5.read(HDF5.attributes(path)["resolution"]))
    data_type = _TYPE_DICT[HDF5.read(HDF5.attributes(path)["data_type"])]
    return Dict(
        "type" => _read_time_series_type(path),
        "initial_timestamp" => initial_timestamp,
        "resolution" => resolution,
        "dataset_size" => size(path["data"]),
        "start_time" => initial_timestamp + resolution * (rows.start - 1),
        "data_type" => data_type,
    )
end

function _read_time_series_type(path)
    module_str = HDF5.read(HDF5.attributes(path)["module"])
    type_str = HDF5.read(HDF5.attributes(path)["type"])
    return get_type_from_strings(module_str, type_str)
end

function add_time_series_reference!(
    storage::Hdf5TimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
    ts_uuid::UUIDs.UUID,
)
    check_read_only(storage)
    uuid = string(ts_uuid)
    component_name = make_component_name(component_uuid, name)
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        path = root[uuid]
        _append_item!(path, "components", component_name)
        @debug "Add reference to existing time series entry." uuid component_uuid name
    end
end

# TODO: This needs to change if we want to directly convert Hdf5TimeSeriesStorage to
# InMemoryTimeSeriesStorage, which is currently not supported System deserialization.
function iterate_time_series(storage::Hdf5TimeSeriesStorage)
    Channel() do channel
        HDF5.h5open(storage.file_path, "r") do file
            root = _get_root(storage, file)
            for uuid_group in root
                uuid_path = HDF5.name(uuid_group)
                uuid_str = uuid_path[(findlast("/", uuid_path).start + 1):end]
                uuid = UUIDs.UUID(uuid_str)

                data = uuid_group["data"][:]
                attributes = Dict()
                for name in keys(HDF5.attributes(uuid_group))
                    attributes[name] = HDF5.read(HDF5.attributes(uuid_group)[name])
                end
                for item in HDF5.read(uuid_group["components"])
                    component, name = deserialize_component_name(item)
                    put!(channel, (component, name, data, attributes))
                end
            end
        end
    end
end

#=
# This could be used if we deserialize the type directly from HDF.
function _make_rows_columns(dataset, ::Type{T}) where T <: StaticTimeSeries
    rows = UnitRange(1, size(dataset)[1])
    columns = UnitRange(1, 1)
    return (rows, columns)
end

function _make_rows_columns(dataset, ::Type{T}) where T <: Forecast
    rows = UnitRange(1, size(dataset)[1])
    columns = UnitRange(1, size(dataset)[2])
    return (rows, columns)
end
=#

function remove_time_series!(
    storage::Hdf5TimeSeriesStorage,
    uuid::UUIDs.UUID,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
)
    check_read_only(storage)
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        path = _get_time_series_path(root, uuid)
        if _remove_item!(path, "components", make_component_name(component_uuid, name))
            @debug "$path has no more references; delete it."
            HDF5.delete_object(path)
        end
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: StaticTimeSeries}
    # Note that all range checks must occur at a higher level.
    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert_op attributes["type"] == T
        @debug "deserializing a StaticTimeSeries" T
        data_type = attributes["data_type"]
        data = get_hdf_array(path["data"], data_type, rows)
        return T(
            ts_metadata,
            TimeSeries.TimeArray(
                range(
                    attributes["start_time"];
                    length = length(rows),
                    step = attributes["resolution"],
                ),
                data,
            ),
        )
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: AbstractDeterministic}
    # Note that all range checks must occur at a higher level.
    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        actual_type = _read_time_series_type(path)
        if actual_type == SingleTimeSeries
            last_index = size(path["data"])[1]
            return deserialize_deterministic_from_single_time_series(
                storage,
                ts_metadata,
                rows,
                columns,
                last_index,
            )
        end

        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert actual_type <: T "actual_type = $actual_type T = $T"
        @debug "deserializing a Forecast" T
        data_type = attributes["data_type"]
        data = get_hdf_array(path["data"], data_type, attributes, rows, columns)
        new_ts = actual_type(ts_metadata, data)
    end
end

function get_hdf_array(
    dataset,
    ::Type{<:CONSTANT},
    attributes::Dict{String, Any},
    rows::UnitRange{Int},
    columns::UnitRange{Int},
)
    data = SortedDict{Dates.DateTime, Vector{Float64}}()
    initial_timestamp = attributes["start_time"]
    interval = attributes["interval"]
    start_time = initial_timestamp + interval * (columns.start - 1)
    if length(columns) == 1
        data[start_time] = dataset[rows, columns.start]
    else
        data_read = dataset[rows, columns]
        for (i, it) in
            enumerate(range(start_time; length = length(columns), step = interval))
            data[it] = @view data_read[1:length(rows), i]
        end
    end
    return data
end

function get_hdf_array(
    dataset,
    type::Type{POLYNOMIAL},
    attributes::Dict{String, Any},
    rows::UnitRange{Int},
    columns::UnitRange{Int},
)
    data = SortedDict{Dates.DateTime, Vector{POLYNOMIAL}}()
    initial_timestamp = attributes["start_time"]
    interval = attributes["interval"]
    start_time = initial_timestamp + interval * (columns.start - 1)
    if length(columns) == 1
        data[start_time] = retransform_hdf_array(dataset[rows, columns.start, :], type)
    else
        data_read = retransform_hdf_array(dataset[rows, columns, :], type)
        for (i, it) in
            enumerate(range(start_time; length = length(columns), step = interval))
            data[it] = @view data_read[1:length(rows), i]
        end
    end
    return data
end

function get_hdf_array(
    dataset,
    type::Type{PWL},
    attributes::Dict{String, Any},
    rows::UnitRange{Int},
    columns::UnitRange{Int},
)
    data = SortedDict{Dates.DateTime, Vector{PWL}}()
    initial_timestamp = attributes["start_time"]
    interval = attributes["interval"]
    start_time = initial_timestamp + interval * (columns.start - 1)
    if length(columns) == 1
        data[start_time] = retransform_hdf_array(dataset[rows, columns.start, :, :], type)
    else
        data_read = retransform_hdf_array(dataset[rows, columns, :, :], type)
        for (i, it) in
            enumerate(range(start_time; length = length(columns), step = interval))
            data[it] = @view data_read[1:length(rows), i]
        end
    end
    return data
end

function get_hdf_array(dataset, type::Type{<:CONSTANT}, rows::UnitRange{Int})
    data = retransform_hdf_array(dataset[rows], type)
    return data
end

function get_hdf_array(dataset, type::Type{POLYNOMIAL}, rows::UnitRange{Int})
    data = retransform_hdf_array(dataset[rows, :, :], type)
    return data
end

function get_hdf_array(dataset, type::Type{PWL}, rows::UnitRange{Int})
    data = retransform_hdf_array(dataset[rows, :, :, :], type)
    return data
end

function retransform_hdf_array(data::Array, ::Type{<:CONSTANT})
    return data
end

function retransform_hdf_array(data::Array, T::Type{POLYNOMIAL})
    row, column, tuple_length = get_data_dims(data, T)
    if isnothing(column)
        t_data = Array{POLYNOMIAL}(undef, row)
        for r in 1:row
            t_data[r] = tuple(data[r, 1:tuple_length]...)
        end
    else
        t_data = Array{POLYNOMIAL}(undef, row, column)
        for r in 1:row, c in 1:column
            t_data[r, c] = tuple(data[r, c, 1:tuple_length]...)
        end
    end
    return t_data
end

function retransform_hdf_array(data::Array, T::Type{PWL})
    row, column, tuple_length, array_length = get_data_dims(data, T)
    if isnothing(column)
        t_data = Array{PWL}(undef, row)
        for r in 1:row
            tuple_array = Array{POLYNOMIAL}(undef, array_length)
            for l in 1:array_length
                tuple_array[l] = tuple(data[r, 1:tuple_length, l]...)
            end
            t_data[r] = tuple_array
        end
    else
        t_data = Array{PWL}(undef, row, column)
        for r in 1:row, c in 1:column
            tuple_array = Array{POLYNOMIAL}(undef, array_length)
            for l in 1:array_length
                tuple_array[l] = tuple(data[r, c, 1:tuple_length, l]...)
            end
            t_data[r, c] = tuple_array
        end
    end
    return t_data
end

function get_data_dims(data::Array, ::Type{POLYNOMIAL})
    if length(size(data)) == 2
        row, tuple_length = size(data)
        return (row, nothing, tuple_length)
    elseif length(size(data)) == 3
        return size(data)
    else
        error("Hdf data array is $(length(size(data)))-D array, expected 2-D or 3-D array.")
    end
end

function get_data_dims(data::Array, ::Type{PWL})
    if length(size(data)) == 3
        row, tuple_length, array_length = size(data)
        return (row, nothing, tuple_length, array_length)
    elseif length(size(data)) == 4
        return size(data)
    else
        error("Hdf data array is $(length(size(data)))-D array, expected 3-D or 4-D array.")
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: Probabilistic}
    # Note that all range checks must occur at a higher level.
    total_percentiles = length(get_percentiles(ts_metadata))

    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert_op attributes["type"] == T
        @assert_op length(attributes["dataset_size"]) == 3
        @debug "deserializing a Forecast" T
        data = SortedDict{Dates.DateTime, Matrix{attributes["data_type"]}}()
        initial_timestamp = attributes["start_time"]
        interval = attributes["interval"]
        start_time = initial_timestamp + interval * (first(columns) - 1)
        if length(columns) == 1
            data[start_time] =
                transpose(path["data"][1:total_percentiles, rows, first(columns)])
        else
            data_read = PermutedDimsArray(
                path["data"][1:total_percentiles, rows, columns],
                [3, 2, 1],
            )
            for (i, it) in enumerate(
                range(start_time; length = length(columns), step = attributes["interval"]),
            )
                data[it] = @view data_read[i, 1:length(rows), 1:total_percentiles]
            end
        end

        new_ts = T(ts_metadata, data)
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: Scenarios}
    # Note that all range checks must occur at a higher level.
    total_scenarios = get_scenario_count(ts_metadata)

    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert_op attributes["type"] == T
        @assert_op length(attributes["dataset_size"]) == 3
        @debug "deserializing a Forecast" T
        data = SortedDict{Dates.DateTime, Matrix{attributes["data_type"]}}()
        initial_timestamp = attributes["start_time"]
        interval = attributes["interval"]
        start_time = initial_timestamp + interval * (first(columns) - 1)
        if length(columns) == 1
            data[start_time] =
                transpose(path["data"][1:total_scenarios, rows, first(columns)])
        else
            data_read =
                PermutedDimsArray(path["data"][1:total_scenarios, rows, columns], [3, 2, 1])
            for (i, it) in enumerate(
                range(start_time; length = length(columns), step = attributes["interval"]),
            )
                data[it] = @view data_read[i, 1:length(rows), 1:total_scenarios]
            end
        end

        new_ts = T(ts_metadata, data)
    end
end

function clear_time_series!(storage::Hdf5TimeSeriesStorage)
    check_read_only(storage)
    # Re-create the file. HDF5 will not actually free up the deleted space until h5repack
    # is run on the file.
    _make_file(storage)
    @info "Cleared all time series."
end

function get_num_time_series(storage::Hdf5TimeSeriesStorage)
    num = 0

    HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        for component in root
            num += 1
        end
    end

    return num
end

function _make_file(storage::Hdf5TimeSeriesStorage)
    HDF5.h5open(storage.file_path, "w") do file
        HDF5.create_group(file, HDF5_TS_ROOT_PATH)
    end
end

_get_root(storage::Hdf5TimeSeriesStorage, file) = file[HDF5_TS_ROOT_PATH]

function _get_time_series_path(root::HDF5.Group, uuid::UUIDs.UUID)
    uuid_str = string(uuid)
    if !haskey(root, uuid_str)
        throw(ArgumentError("UUID $uuid_str does not exist"))
    end

    return root[uuid_str]
end

function _append_item!(path::HDF5.Group, name::AbstractString, value::AbstractString)
    handle = HDF5.open_object(path, name)
    values = HDF5.read(handle)
    HDF5.close(handle)

    if in(value, values)
        # Just in case any component tries to store the same reference twice.
        return nothing
    else
        push!(values, value)
    end

    # Delete and re-write.
    HDF5.delete_object(path, name)
    path[name] = values
    @debug "Appended $value to $name" values
end

"""
Removes value from the dataset called name.
Returns true if the array is empty afterwards.
"""
function _remove_item!(path::HDF5.Group, name::AbstractString, value::AbstractString)
    handle = HDF5.open_object(path, name)
    vals = HDF5.read(handle)
    HDF5.close(handle)

    orig_len = length(vals)
    filter!(x -> x != value, vals)
    exp_len = orig_len - 1
    if length(vals) != exp_len
        throw(
            ArgumentError(
                "$value wasn't stored in $name or was stored more than once. " *
                "exp_len = $exp_len actual = $(length(vals))",
            ),
        )
    end

    # Delete and rewrite.
    # This is not efficient, but this is expected to be uncommon and not to have
    # large counts.
    HDF5.delete_object(path, name)
    path[name] = vals

    @debug "Removed $value from $name" vals
    return isempty(vals)
end

function check_read_only(storage::Hdf5TimeSeriesStorage)
    if storage.read_only
        error("Operation not permitted; this time series file is read-only")
    end
end

is_read_only(storage::Hdf5TimeSeriesStorage) = storage.read_only

function compare_values(x::Hdf5TimeSeriesStorage, y::Hdf5TimeSeriesStorage)::Bool
    item_x = sort!(collect(iterate_time_series(x)), by = z -> z[1])
    item_y = sort!(collect(iterate_time_series(y)), by = z -> z[1])
    if length(item_x) != length(item_y)
        @error "lengths don't match" length(item_x) length(item_y)
        return false
    end

    for ((uuid_x, name_x, data_x, attrs_x), (uuid_y, name_y, data_y, attrs_y)) in
        zip(item_x, item_y)
        if uuid_x != uuid_y
            @error "component UUIDs don't match" uuid_x uuid_y
            return false
        end
        if name_x != name_y
            @error "names don't match" name_x name_y
            return false
        end
        if data_x != data_y
            @error "data doesn't match" data_x data_y
            return false
        end
        if sort!(collect(keys(attrs_x))) != sort!(collect(keys(attrs_y)))
            @error "attr keys don't match" attrs_x attrs_y
        end
        if collect(values(attrs_x)) != collect(values(attrs_y))
            @error "attr values don't match" attrs_x attrs_y
        end
    end

    return true
end


text:
"""
$(TYPEDEF)

One sample superiority test for proportion

Constructors
------------
* `OneSamplePropSuperior(p1::Real, p0::Real, delta::Real)`

Arguments
---------
* `p1`: True proportion

* `p0`: Proportion under null hypothesis

* `delta`: Superiority margin

Fields
------
$(FIELDS)
"""
type OneSamplePropSuperior <: TrialTest
    p1::Real
    p0::Real
    delta::Real

    # Validator
    function OneSamplePropSuperior(p1, p0, delta)

        if !((0 < p1 < 1) & (0 < p0 < 1))

            error("Proportion value must be in (0, 1)")

        end # end if

        if !(0 <= delta < 1)

            error("The superiority margin δ must be in [0, 1)")

        end # end if
        new(p1, p0, delta)

    end # end function

end # type

# One sample superiority test for proportion
function hypotheses{T <: OneSamplePropSuperior}(test::T, n::Real, std::Void, alpha::Real, side::String)

    diff = test.p1 - test.p0 - test.delta
    se = sqrt(1 / n * test.p1 * (1 - test.p1))
    z = diff / se
    p = cdf(ZDIST, z - quantile(ZDIST, 1 - alpha)) + cdf(ZDIST, -z - quantile(ZDIST, 1 - alpha))
    return p

end # end function


text:
"""
```
    AlphaNumeral
```
Numbers represented as alphabets. ex. A, B, C,... from 27 onwards AA, BB, CC etc.
"""
struct AlphaNumeral <: Integer
    val::Int
    str::String
end

"""
```
    AlphaNumeral(str::String)
    AlphaNumeral(n::Int)
```
Constructors for `AlphaNumeral`.
"""
AlphaNumeral(str::String) = parse(AlphaNumeral, str)
AlphaNumeral(n::T) where T <: Integer = convert(AlphaNumeral, n)

Base.hash(num::AlphaNumeral) = xor(hash(num.str), hash(num.val))

Base.typemax(::Type{AlphaNumeral}) = 156 # ZZZZZZ
Base.typemin(::Type{AlphaNumeral}) = 1

"""
```
    @an_str(str)
```
String decorator for `AlphaNumeral` definitions.

#Example
```
julia> an"AA"
LabelNumerals.AlphaNumeral(27, "AA")
```
"""
macro an_str(str)
    AlphaNumeral(str)
end

Base.convert(::Type{T}, num::AlphaNumeral) where T <: Integer =
    convert(T, num.val)

function Base.parse(::Type{AlphaNumeral}, str::String)
    s = uppercase(str)
    c = s[1]
    cnt = 0
    for a in s
        a != c && throw(DomainError(-2, "Characters must be repeated"))
        cnt += 1
    end
    val = 26*(cnt - 1) + (c - 'A') + 1
    typemin(AlphaNumeral) <= val <= typemax(AlphaNumeral) &&
        return AlphaNumeral(val, str)
    throw(DomainError(-1, "Value out of range"))
end

function Base.convert(::Type{AlphaNumeral}, val::Int)
    if typemin(AlphaNumeral) <= val <= typemax(AlphaNumeral)
        n = div(val, 26)
        r = rem(val, 26)
        if r == 0
            r = 26
            n -= 1
        end
        str = string(fill(Char('A' + r -1),(n+1))...)
        return AlphaNumeral(val, str)
    end
    throw(DomainError(-1, "Value out of range"))
end


text:
module Correlation_Functions
export form_atomic_time_series, autocorrelation_function, custom_acf, plot_correlation_function, vdos, plot_vdos
using FFTW, DSP, Plots, LaTeXStrings, ProgressBars

    function form_atomic_time_series(frames::AbstractArray)
        return transpose(hcat(vec.(frames)...))
    end

    function autocorrelation_function(time_series::AbstractArray, block_width::Int, num_lags_per_window::Int=20, max_blocks::Int=0)
        """
        Compute the autocorrelation function of a single time series or many time series computed
        as a matrix of time series where each column is a time series. This will average the acf over
        as many disconnected blocks of length block_width as possible, unless the optional max_blocks is specified.
        """
        signal = zeros(block_width)
        num_blocks = convert(Int, floor(size(time_series, 1) / block_width))
        if max_blocks > 0 && max_blocks < num_blocks
            num_blocks = max_blocks
        end
        block_stride = convert(Int, floor(block_width / 20))

        for col in ProgressBar(eachcol(time_series))
            for i_block in 1:num_blocks
                for i_slide in 1:block_width:block_stride
                    block = @inbounds @view col[(i_block-1)*block_width+1+(i_slide-1):i_block*block_width+(i_slide-1)]
                    signal += DSP.xcorr(block, block)[block_width:end]
                end
                signal /= num_lags_per_window
            end
            signal /= num_blocks
        end
        return signal / signal[begin]
    end

    function plot_correlation_function(time_series::AbstractArray, frame_spacing::Float64)
        plotlyjs()
        display(plot((1:length(time_series)) * frame_spacing, time_series,
        label="ACF Signal",
        xlims=(0, 2500),
        xlabel="Time (fs)",
        ylabel="Signal Correlation"))
    end

    function vdos(signal::AbstractArray, frame_spacing::Float64)
        """
        Computes the vibrational density of states from the velocity autocorrelation function.
        Args:
            signal: 1-D array representing some time series
            frame_spacing: float which is the spacing in femtoseconds of the frames
        Return: both the fourier frequencies and the signal intensity as a tuple.
        """
        freqs = FFTW.rfftfreq(length(signal), 1 / frame_spacing)[1:floor(Int64, length(signal)/2)] * 10^(15) / (2.9979*(10^10)) # fs to cm^-1
        vdos = real.(FFTW.rfft(signal)[1:floor(Int64, length(signal)/2)])
        return freqs, vdos
    end

    function plot_vdos(frequencies::AbstractArray, vdos::AbstractArray)
        plotlyjs()
        plot(frequencies, vdos,
        label="VDOS",
        xlims=(1000, 4400),
        ylims=(0, 15),
        xlabel=(L"Frequency (cm^{-1})"),
        ylabel=("Intensity"))
    end

end

max_stars_repo_name: heindelj/pyMD

max_issues_repo_name: heindelj/pyMD

max_forks_repo_name: heindelj/pyMD


text:
using NbodyGradient
import NbodyGradient: kepler_init, init_nbody

maxabs(x) = maximum(abs.(x))

# For backwards compatability
include("loglinspace.jl")
if VERSION >= v"0.7"
  using Test
  using LinearAlgebra
  using Statistics
  using DelimitedFiles
else
  using Base.Test
end

@testset "NbodyGradient" begin
    print("Initial Conditions... ")
    @testset "Initial Conditions" begin
        include("test_kepler_init.jl")
        include("test_init_nbody.jl")
        include("test_initial_conditions.jl")
    end;
    println("Finished.")

    print("Integrator... ")
    @testset "Integrator" begin
        include("test_kickfast.jl")
        include("test_phic.jl")
        include("test_kepler_driftij_gamma.jl")
        include("test_phisalpha.jl")
        include("test_integrator.jl")
    end;
    println("Finished.")

    print("Outputs... ")
    @testset "Outputs" begin
        include("test_cartesian_to_elements.jl")
    end
    println("Finished.")

    print("TTVs... ")
    @testset "TTVs" begin
        include("test_findtransit.jl")
        include("test_transit_timing.jl")
        include("test_transit_parameters.jl")
    end;
    println("Finished.")

    print("RVs... ")
    @testset "RVs" begin
        include("test_radial_velocity.jl")
    end;
    println("Finished.")
end


text:
"""
Beautiful plots by default. To plot on an existing Axes, pass it as the first keyword
argument. Keyword arguments that apply to `Line2D`s are passed to `ax.plot`. The rest are
passed to `set`.
"""
function plot(args...; kw...)
    :data in keys(kw) && error("'data' keyword not supported.")
    args = [args...]  # Tuple to Vector (so we can `pop!`)
    if first(args) isa PyObject && pyisinstance(first(args), mpl.axes.Axes)
        ax = popfirst!(args)
    else
        ax = plt.gca()
    end
    plotkw = Dict(k => v for (k, v) in kw if hasproperty(mpl.lines.Line2D, "set_$k"))
    otherkw = Dict(k => v for (k, v) in kw if k ∉ keys(plotkw))
    ax.plot((args .|> as_mpl_type)...; (plotkw |> mapvals $ as_mpl_type)...)
    _handle_units!(ax, args)  # Mutating, because `_extract_plotted_data!` peels off `args`
                              # until it's empty.
    set(ax; otherkw...)
    return ax
end

function _handle_units!(ax, plotargs)
    xs, ys = _extract_plotted_data!(plotargs)
    for (arrays, axis) in zip([xs, ys], [ax.xaxis, ax.yaxis])
        isempty(arrays) && continue
        for array in arrays
            has_mixed_dimensions(array) &&
                error("Argument has mixed dimensions: $array")
        end
        arrays_dimensions = dimension.(arrays)
        all(isequal(first(arrays_dimensions)), arrays_dimensions) ||
            error("Not all $(axis.axis_name)-axis arrays have the same dimensions: $arrays_dimensions.")
        # Store units as a new property on the array object. Note that `units` property
        # already exists in Mpl.
        axis.unitful_units = unit(eltype(first(arrays)))
    end
end

has_mixed_dimensions(x::AbstractArray{<:Quantity{T,Dims}}) where {T,Dims}  = false
has_mixed_dimensions(x::AbstractArray{<:Quantity})                         = true
has_mixed_dimensions(x::AbstractArray)                                     = false

function _extract_plotted_data!(plotargs)
    # Process `ax.plot`'s vararg by peeling off the front: [x], y, [fmt].
    # Based on https://github.com/matplotlib/matplotlib/blob/710fce/lib/matplotlib/axes/_base.py#L304-L312
    xs = []
    ys = []
    while !isempty(plotargs)
        if length(plotargs) == 1
            push!(ys, popfirst!(plotargs))
        else
            a = popfirst!(plotargs)
            b = popfirst!(plotargs)
            if b isa AbstractString  # fmt string
                push!(ys, a)
            else
                push!(xs, a)
                push!(ys, b)
                if !isempty(plotargs) && first(plotargs) isa AbstractString
                    popfirst!(plotargs)
                end
            end
        end
    end
    return asarray.(xs), asarray.(ys)
end

asarray(x::Number) = fill(x)  # → zero-dimensional array.
asarray(x::AbstractArray) = x