text:
using DustExtinction
using Test, Measurements, Unitful, UnitfulAstro, Random

Random.seed!(9994445781)

include("color_laws.jl")
include("dust_maps.jl")

@testset "redden" begin
    wave = [3372., 4404., 5428., 6509., 8090.,
            3683., 4393., 5519., 6602., 8046.,
            12660., 16732., 22152.,
            5244., 6707., 7985., 9055.,
            6993.,
            3502., 4676., 4127.,
            4861., 5479.,
            3546., 4925., 6335., 7799., 9294.,
            3047., 4711., 5498.,
            6042., 7068., 8066.,
            4814., 6571., 8183.]

    # Av=0.3, Rv=3.1, law=ccm89
    ref_values = [0.6376288197244566, 0.6937993158943373, 0.7556003129054399, 0.7959216081971665, 0.8508142213648981, 0.6517368195297566, 0.6930246485833793, 0.7596467718783736, 0.7990835963390542, 0.8493553321733043, 0.9264835715958143, 0.9524307890106202, 0.9694546832913796, 0.7467375901587614, 0.8026674093639407, 0.847310247751628, 0.8771971841412943, 0.8125473475954554, 0.6442645437827309, 0.7127773090414787, 0.6749407318003154, 0.7249166081334651, 0.7578933014813821, 0.6461535564599767, 0.7288900774623704, 0.7900092679188825, 0.8409405763429691, 0.881972061390043, 0.6091986958411011, 0.7151404275193919, 0.7587309863781261, 0.7798780926274524, 0.8151697825084553, 0.8500202342098636, 0.7219200730427752, 0.7980286586611023, 0.8538476253948861]

    flux = ones(length(wave))
    output = @inferred broadcast(redden, flux, wave, 0.3)
    @test output ≈ ref_values
    @test @inferred(broadcast(deredden, output, wave, 0.3)) ≈ flux
    map!((f, w)->redden(f, w, 0.3), flux, flux, wave)
    @test flux ≈ ref_values

    # Measurements
    flux = ones(length(wave)) .± 0.1
    output = @inferred broadcast(redden, flux, wave, 0.3)
    @test Measurements.value.(output) ≈ ref_values
    @test @inferred(broadcast(deredden, output, wave, 0.3)) ≈ flux
    map!((f, w)->redden(f, w, 0.3), flux, flux, wave)
    @test flux ≈ ref_values

    # Unitful
    flux = ones(length(wave))u"Jy"
    wave = wave * u"angstrom"
    output = @inferred broadcast(redden, flux, wave, 0.3)
    @test ustrip.(output) ≈ ref_values
    @test @inferred(broadcast(deredden, output, wave, 0.3)) ≈ flux
    map!((f, w)->redden(f, w, 0.3), flux, flux, wave)
    @test flux ≈ output
end


text:
export Day9, Day9_1, Day9_2
using DataStructures: CircularBuffer
abstract type Day9 <: Day end
abstract type Day9_1 <: Day9 end
abstract type Day9_2 <: Day9 end

line_parser(::Type{<:Day9}, line) = line |> strip |> l -> parse(Int, l)

solve(::Type{T}, lines::Vector) where T <: Day9 = solve(T, lines, 25)

solve(::Type{T}, lines::Vector, pre_length) where T <: Day9_1 = begin
    buffer = CircularBuffer{Int}((pre_length-1)*pre_length)
    # Fill with sum of first "pre_length" numbers
    for i = 1:pre_length
        for j= 1:pre_length
            if i == j
                continue
            end
            push!(buffer, lines[i]+lines[j])
        end
    end
    for i = pre_length+1:length(lines)
        if lines[i] in buffer
            for j=i-pre_length:i-1
                if i == j
                    continue
                end
                push!(buffer, lines[i]+lines[j])
            end
        else
            return lines[i]
        end
    end
end

solve(::Type{T}, lines::Vector, pre_length) where T <: Day9_2 = begin
    num_sol = solve(Day9_1, lines, pre_length)
    for i = 1:length(lines)-1
        agg = lines[i]
        for j = i+1:length(lines)
            agg = agg + lines[j]
            if agg == num_sol
                return minimum(lines[i:j]) + maximum(lines[i:j])
            elseif agg > num_sol
                break
            end
        end
    end
end

ext: jl


text:
abstract AbstractModel{T}

description(m::AbstractModel) = "DSGE Model $(spec(m)) (subspec $(subspec(m)))"

function Base.show{T<:AbstractModel}(io::IO, m::T)
    @printf io "Dynamic Stochastic General Equilibrium Model\n"
    @printf io "%s\n" T
    @printf io "description:\n %s\n"          description(m)
    @printf io "no. states:             %i\n" n_states(m)
    @printf io "no. anticipated shocks: %i\n" n_anticipated_shocks(m)
    @printf io "data vintage:           %s\n" data_vintage(m)
end

@inline function Base.getindex(m::AbstractModel, i::Integer)
    if i <= (j = length(m.parameters))
        return m.parameters[i]
    else
        return m.steady_state[i-j]
    end
end

# need to define like this so we can disable bounds checking
@inline function Base.getindex(m::AbstractModel, k::Symbol)
    i = m.keys[k]
    @inbounds if i <= (j = length(m.parameters))
        return m.parameters[i]
    else
        return m.steady_state[i-j]
    end
end

@inline function Base.setindex!{T<:Number}(m::AbstractModel, value::T, i::Integer)
    if i <= (j = length(m.parameters))
        param = m.parameters[i]
        param.value = value
        if isa(param, ScaledParameter)
            param.scaledvalue = param.scaling(value)
        end
        return param
    else
        steady_state_param = m.steady_state[i-j]
        steady_state_param.value = value
        return steady_state_param
    end
end

"""
```
setindex!{T<:AbstractParameter}(m::AbstractModel, param::T, i::Integer)
```

If `i`<length(m.parameters), overwrites m.parameters[i] with
param. Otherwise, overwrites m.steady_state[i-length(m.parameters).
"""
@inline function Base.setindex!{T<:AbstractParameter}(m::AbstractModel, param::T, i::Integer)
    if i <= (j = length(m.parameters))
        m.parameters[i] = param
    else
        m.steady_state[i-j] = param
    end
    return param
end

Base.setindex!(m::AbstractModel, value, k::Symbol) = Base.setindex!(m, value, m.keys[k])

"""
```
(<=){T}(m::AbstractModel{T}, p::AbstractParameter{T})
```

Syntax for adding a parameter to a model: m <= parameter.
NOTE: If `p` is added to `m` and length(m.steady_state) > 0, `keys(m)` will not generate the
index of `p` in `m.parameters`.
"""
function (<=){T}(m::AbstractModel{T}, p::AbstractParameter{T})

    if !in(p.key, keys(m.keys))

        new_param_index = length(m.keys) + 1

        # grow parameters and add the parameter
        push!(m.parameters, p)

        # add parameter location to dict
        setindex!(m.keys, new_param_index, p.key)
    else
        # overwrite the previous parameter with the new one
        setindex!(m, p, p.key)
    end
end

"""
```
(<=){T}(m::AbstractModel{T}, ssp::SteadyStateParameter)
```

Add a new steady-state value to the model by appending `ssp` to the `m.steady_state` and
adding `ssp.key` to `m.keys`.
"""
function (<=){T}(m::AbstractModel{T}, ssp::SteadyStateParameter)

    if !in(ssp.key, keys(m.keys))
        new_param_index = length(m.keys) + 1

        # append ssp to steady_state vector
        push!(m.steady_state, ssp)

        # add parameter location to dict
        setindex!(m.keys, new_param_index, ssp.key)
    else
        # overwrite the previous parameter with the new one
        setindex!(m, ssp, ssp.key)
    end
end

Distributions.logpdf(m::AbstractModel) = logpdf(m.parameters)
Distributions.pdf(m::AbstractModel) = exp(logpdf(m))

# Number of anticipated policy shocks
n_anticipated_shocks(m::AbstractModel) = get_setting(m, :n_anticipated_shocks)

# Padding for number of anticipated policy shocks
n_anticipated_shocks_padding(m::AbstractModel) = get_setting(m, :n_anticipated_shocks_padding)

# Index into data matrix of first period to incorporate expected rate data
function zlb_start_index(m::AbstractModel)
    zlb_start_quarter = get_setting(m, :date_zlbregime_start)
    presample_start_quarter = get_setting(m, :date_presample_start)
    quarters = subtract_quarters(zlb_start_quarter, presample_start_quarter)
    index = quarters+1
    return index
end

# Number of presample periods
function n_presample_periods(m::AbstractModel)
    mainsample_start_quarter = get_setting(m, :date_mainsample_start)
    presample_start_quarter = get_setting(m, :date_presample_start)
    quarters = subtract_quarters(mainsample_start_quarter, presample_start_quarter)
    return quarters
end

# Number of a few things that are useful
n_states(m::AbstractModel)                 = length(m.endogenous_states)
n_states_augmented(m::AbstractModel)       = n_states(m) + length(m.endogenous_states_augmented)
n_shocks_exogenous(m::AbstractModel)       = length(m.exogenous_shocks)
n_shocks_expectational(m::AbstractModel)   = length(m.expected_shocks)
n_equilibrium_conditions(m::AbstractModel) = length(m.equilibrium_conditions)
n_observables(m::AbstractModel)            = length(m.observables)
n_parameters(m::AbstractModel)             = length(m.parameters)
n_parameters_steady_state(m::AbstractModel)= length(m.steady_state)
n_parameters_free(m::AbstractModel)        = sum([!α.fixed for α in m.parameters])

# Interface for I/O settings
spec(m::AbstractModel)         = m.spec
subspec(m::AbstractModel)      = m.subspec
saveroot(m::AbstractModel)     = get_setting(m, :saveroot)
dataroot(m::AbstractModel)     = get_setting(m, :dataroot)

# Interface for data
data_vintage(m::AbstractModel) = get_setting(m, :data_vintage)
cond_vintage(m::AbstractModel) = get_setting(m, :cond_vintage)
cond_id(m::AbstractModel) = get_setting(m, :cond_id)

# Interface for general computation settings
use_parallel_workers(m::AbstractModel)    = get_setting(m, :use_parallel_workers)

# Interface for data step
use_population_forecast(m::AbstractModel) = get_setting(m, :use_population_forecast)

# Interface for estimation settings
reoptimize(m::AbstractModel)            = get_setting(m, :reoptimize)
calculate_hessian(m::AbstractModel)     = get_setting(m, :calculate_hessian)
hessian_path(m::AbstractModel)          = get_setting(m, :hessian_path)
n_hessian_test_params(m::AbstractModel) = get_setting(m, :n_hessian_test_params)

# Interface for Metropolis-Hastings settings
n_mh_blocks(m::AbstractModel)      =  get_setting(m, :n_mh_blocks)
n_mh_simulations(m::AbstractModel) =  get_setting(m, :n_mh_simulations)
n_mh_burn(m::AbstractModel)        =  get_setting(m, :n_mh_burn)
mh_thin(m::AbstractModel)          =  get_setting(m, :mh_thin)

"""
```
load_parameters_from_file(m::AbstractModel,path::AbstractString)
```
Returns a vector of parameters, read from a file, suitable for updating `m`.
"""
function load_parameters_from_file(m::AbstractModel, path::AbstractString)

    if isfile(path) && splitext(path)[2] == ".h5"
        x  = h5open(path, "r") do file
            try
                read(file, "params")
            catch
                error("$path does not contain variable params")
            end
        end
    else
        error("$path is not a valid HDF5 file.")
    end

    @assert length(x) == length(m.parameters)
    @assert eltype(x) == typeof(m.parameters[1].value)
    return x
end

"""
```
specify_mode!(m::AbstractModel, mode_file::AbstractString=""; verbose=:low)
```

Updates the values of `m.parameters` with the values from
`mode_file`. Sets `reoptimize` setting to `false`.

Usage: should be run before calling `estimate(m)`, e.g.:

    m = Model990()
    specify_mode!(m, modefile)
    estimate(m)
"""
function specify_mode!(m::AbstractModel, mode_file::AbstractString = ""; verbose=:low)

    m <= Setting(:reoptimize, false)

    if mode_file == ""
        mode_file = inpath(m, "user", "paramsmode.h5")
    end
    mode_file = normpath(mode_file)

    update!(m,load_parameters_from_file(m,mode_file))

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Loaded previous mode from $mode_file.")
    end

end

"""
```
specify_hessian(m::AbstractModel, path::AbstractString=""; verbose=:low)
```

Specify a Hessian matrix calculated at the posterior mode to use in the model estimation. If
no path is provided, will attempt to detect location.
"""
function specify_hessian(m::AbstractModel, path::AbstractString=""; verbose=:low)
    if isempty(path)
        path = inpath(m, "user", "hessian.h5")
    end

    if isfile(path) && splitext(path)[2] == ".h5"
        m <= Setting(:hessian_path, normpath(abspath(path)))
    else
        error("Invalid input Hessian file: $path",)
    end

    m <= Setting(:calculate_hessian, false)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Specified hessian from $path.")
    end
end

#=
Build paths to where input/output/results data are stored.

Description:
Creates the proper directory structure for input and output files, treating the DSGE/save
    directory as the root of a savepath directory subtree. Specifically, the following
    structure is implemented:

    dataroot/

    savepathroot/
                 output_data/<spec>/<subspec>/log/
                 output_data/<spec>/<subspec>/<out_type>/raw/
                 output_data/<spec>/<subspec>/<out_type>/work/
                 output_data/<spec>/<subspec>/<out_type>/tables/
                 output_data/<spec>/<subspec>/<out_type>/figures/

Note: we refer to the savepathroot/output_data/<spec>/<subspec>/ directory as saveroot.
=#
"""
```
logpath(model)
```
Returns path to log file. Path built as
```
<output root>/output_data/<spec>/<subspec>/log/log_<filestring>.log
```
"""
function logpath(m::AbstractModel)
    return savepath(m, "log", "log.log")
end

strs = [:work, :raw, :tables, :figures]
fns = [symbol(x, "path") for x in strs]
for (str, fn) in zip(strs, fns)
    @eval begin
        # First eval function
        function $fn{T<:AbstractString}(m::AbstractModel,
                                             out_type::T,
                                             file_name::T="",
                                             filestring_addl::Vector{T}=Vector{T}())
            return savepath(m, out_type, $(string(str)), file_name, filestring_addl)
        end

        # Then, add docstring to it
        @doc $(
        """
        ```
        $fn{T<:AbstractString}(m::AbstractModel, out_type::T, file_name::T="")
        ```

        Returns path to specific $str output file, creating containing directory as needed. If
        `file_name` not specified, creates and returns path to containing directory only. Path built
        as
        ```
        <output root>/output_data/<spec>/<subspec>/<out_type>/$str/<file_name>_<filestring>.<ext>
        ```
        """
        ) $fn
    end
end

# Not exposed to user. Actually create path and insert model string to file name.
function savepath{T<:AbstractString}(m::AbstractModel,
                                     out_type::T,
                                     sub_type::T,
                                     file_name::T="",
                                     filestring_addl::Vector{T}=Vector{T}())
    # Containing dir
    path = joinpath(saveroot(m), "output_data", spec(m), subspec(m), out_type, sub_type)
    if !isdir(path)
        mkpath(path)
    end

    # File with model string inserted
    if !isempty(file_name)
        if isempty(filestring_addl)
            myfilestring = filestring(m)
        else
            myfilestring = filestring(m, filestring_addl)
        end
        (base, ext) = splitext(file_name)
        file_name_detail = base * myfilestring * ext
        path = joinpath(path, file_name_detail)
    end

    return path
end

# Input data handled slightly differently, because it is not model-specific.
"""
```
inpath{T<:AbstractString}(m::AbstractModel, in_type::T, file_name::T="")
```

Returns path to specific input data file, creating containing directory as needed. If
`file_name` not specified, creates and returns path to containing directory only. Valid
`in_type` includes:

* `"data"`: recorded data
* `"cond"`: conditional data - nowcasts for the current forecast quarter, or related
* `"user"`: user-supplied data for starting parameter vector, hessian, or related

Path built as
```
<data root>/<in_type>/<file_name>
```
"""
function inpath{T<:AbstractString}(m::AbstractModel, in_type::T, file_name::T="")
    path = dataroot(m)
    # Normal cases.
    if in_type == "data" || in_type == "cond"
        path = joinpath(path, in_type)
    # User-provided inputs. May treat this differently in the future.
    elseif in_type == "user"
        path = joinpath(path, "user")
    else
        error("Invalid in_type: ", in_type)
    end

    # Containing dir
    if !isdir(path)
        mkpath(path)
    end

    # If file_name provided, return full path
    if !isempty(file_name)
        path = joinpath(path, file_name)
    end

    return path
end

filestring(m::AbstractModel) = filestring(m, Vector{AbstractString}())
filestring(m::AbstractModel, d::AbstractString) = filestring(m, [d])
function filestring{T<:AbstractString}(m::AbstractModel,
                                        d::Vector{T})
    if !m.testing
        filestrings = Vector{T}()
        for (skey, sval) in m.settings
            if sval.print
                push!(filestrings, to_filestring(sval))
            end
        end
        append!(filestrings, d)
        sort!(filestrings)
        return "_"*join(filestrings, "_")
    else
        return "_test"
    end
end

"""
```
transform_to_model_space!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
```

Transforms `values` from the real line to the model space, and assigns `values[i]` to
`m.parameters[i].value` for non-steady-state parameters. Recomputes the steady-state
paramter values.

### Arguments
-`m`: the model object
-`values`: the new values to assign to non-steady-state parameters.
"""
function transform_to_model_space!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
    new_values = transform_to_model_space(m.parameters, values)
    update!(m, new_values)
    steadystate!(m)
end

"""
```
update!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
```

Update `m.parameters` with `values`, recomputing the steady-state parameter values.

### Arguments:
-`m`: the model object
-`values`: the new values to assign to non-steady-state parameters.
"""
function update!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
    update!(m.parameters, values)
    steadystate!(m)
end

"""
```
rand{T<:AbstractFloat, U<:AbstractModel}(d::DegenerateMvNormal, m::U; cc::T = 1.0)
```

Generate a draw from d with variance optionally scaled by cc^2.
"""
function rand{T<:AbstractFloat, U<:AbstractModel}(d::DegenerateMvNormal, m::U; cc::T = 1.0)
    return d.μ + cc*d.σ*randn(m.rng, length(d))
end


text:
using JLD
using SparseArrays
using Recommendation

R = load("data/ml-100k.jld")["R"]
data = DataAccessor(R)

recall = cross_validation(
                          1,            # N-fold
                          Recall,       # Metric
                          5,            # Top-k
                          MostPopular,  # Recommender
                          data          # Data Accessor
                         )
println(recall)


text:
struct StaticTeam <: AbstractArray{StaticPokemon, 1}
    mon1::StaticPokemon
    mon2::StaticPokemon
    mon3::StaticPokemon
end

Base.size(s::StaticTeam) = (3,)
Base.IndexStyle(::Type{<:StaticTeam}) = IndexLinear()
Base.getindex(s::StaticTeam, i::UInt8) =
    i == 0x01 ? s.mon1 : i == 0x02 ? s.mon2 : s.mon3

StaticTeam(mons::Array{String}; league::String = "great", cup::String = "all") =
    StaticTeam(StaticPokemon(mons[1], league = league, cup = cup),
               StaticPokemon(mons[2], league = league, cup = cup),
               StaticPokemon(mons[3], league = league, cup = cup))

StaticTeam(mons::Array{Int64}; league::String = "great", cup::String = "all") =
    StaticTeam(StaticPokemon(mons[1], league = league, cup = cup),
               StaticPokemon(mons[2], league = league, cup = cup),
               StaticPokemon(mons[3], league = league, cup = cup))

struct DynamicTeam <: AbstractArray{DynamicPokemon, 1}
    #These values are initialized, but change throughout the battle
    mon1::DynamicPokemon
    mon2::DynamicPokemon
    mon3::DynamicPokemon
    switchCooldown::Int8     # Initially 0
    data::UInt8              # StatBuff info and shields
end

Base.size(d::DynamicTeam) = (3,)
Base.IndexStyle(::Type{<:DynamicTeam}) = IndexLinear()
Base.getindex(d::DynamicTeam, i::UInt8) =
    i == 0x01 ? d.mon1 : i == 0x02 ? d.mon2 : d.mon3

has_shield(team::DynamicTeam) = !iszero(team.data % 0x03)

DynamicTeam(team::StaticTeam) = DynamicTeam(
    DynamicPokemon(team[0x01]),
    DynamicPokemon(team[0x02]),
    DynamicPokemon(team[0x03]),
    Int8(0),
    UInt8(122)
)

# returns 2 if 1, 1 if 2. Note: tested some bit twiddling code that was roughly
# equivalent speed, but this is more readable.
get_other_agent(agent::UInt8) = agent == 0x01 ? 0x02 : 0x01


text:
using Documenter, CategoricalArrays

# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

makedocs(
    modules = [CategoricalArrays],
    sitename = "CategoricalArrays",
    format = Documenter.HTML(canonical = "https://juliadata.github.io/CategoricalArrays.jl/stable/"),
    pages = Any[
        "Overview" => "index.md",
        "Using CategoricalArrays" => "using.md",
        "Implementation details" => "implementation.md",
        "API index" => "apiindex.md"
        ],
    checkdocs = :exports,
    strict=true
)

deploydocs(
    repo = "github.com/JuliaData/CategoricalArrays.jl.git"
)


text:
# This files contains all the gradient definitions for the game.

using Zygote
using Zygote: @adjoint, @nograd, literal_getproperty, @showgrad, hook
import Zygote: accum
import Base: sign, floor

@nograd get_obs, reset_ball!#wall_collision!, #, paddle_collision!

@adjoint Vector_(x, y) = Vector_(x, y), ∆ -> (∆.x, ∆.y)

@adjoint literal_getproperty(v::Vector_, ::Val{:x}) =
    getproperty(v, :x), ∆ -> (Vector_(∆, 0f0), )

@adjoint literal_getproperty(v::Vector_, ::Val{:y}) =
    getproperty(v, :y), ∆ -> (Vector_(0f0, ∆), )

@adjoint Player(length, position, velocity, score) =
    Player(length, position, velocity, score), ∆ -> (∆.length, ∆.position, ∆.velocity, ∆.score)

@adjoint literal_getproperty(p::Player, ::Val{:length}) =
    getproperty(p, :length), ∆ -> (Player(∆, zero(p.position), 0, 0), )

@adjoint literal_getproperty(p::Player, ::Val{:position}) =
    getproperty(p, :position), ∆ -> (Player(0f0, ∆, 0, 0), )

@adjoint literal_getproperty(p::Player, ::Val{:velocity}) =
    getproperty(p, :velocity), ∆ -> (Player(0f0, zero(p.position), ∆, 0), )

@adjoint literal_getproperty(p::Player, ::Val{:score}) =
    getproperty(p, :score), ∆ -> (Player(0f0, zero(p.position), 0, ∆), )

@adjoint Ball(position, velocity, radius) =
    Ball(position, velocity, radius), ∆ -> (∆.position, ∆.velocity, ∆.radius)

@adjoint literal_getproperty(b::Ball, ::Val{:position}) =
    getproperty(b, :position), ∆ -> (Ball(∆, zero(b.velocity), 0f0), )

@adjoint literal_getproperty(b::Ball, ::Val{:velocity}) =
    getproperty(b, :velocity), ∆ -> (Ball(zero(b.position), ∆, 0f0), )

@adjoint literal_getproperty(b::Ball, ::Val{:radius}) =
    getproperty(b, :radius), ∆ -> (Ball(zero(b.position), zero(b.velocity), ∆), )

@adjoint Arena(dims, margin) = Arena(dims, margin), ∆ -> (∆.dims, ∆.margin)

@adjoint literal_getproperty(a::Arena, ::Val{:dims}) =
    getproperty(a, :dims), ∆ -> (Arena(∆, zero(a.margin)), )

@adjoint literal_getproperty(a::Arena, ::Val{:margin}) =
    getproperty(a, :margin), ∆ -> (Arena(zero(a.dims), ∆), )

@adjoint Env(arena, player_a, player_b, ball) =
    Env(arena, player_a, player_b, ball), ∆ -> (∆.arena, ∆.player_a, ∆.player_b, ∆.ball)

@adjoint literal_getproperty(e::Env, ::Val{:arena}) =
    getproperty(e, :arena), ∆ -> (Env(∆, zero(e.player_a), zero(e.player_b), zero(e.ball)), )

@adjoint literal_getproperty(e::Env, ::Val{:player_a}) =
    getproperty(e, :player_a), ∆ -> (Env(zero(e.arena), ∆, zero(e.player_b), zero(e.ball)), )

@adjoint literal_getproperty(e::Env, ::Val{:player_b}) =
    getproperty(e, :player_b), ∆ -> (Env(zero(e.arena), zero(e.player_a), ∆, zero(e.ball)), )

@adjoint literal_getproperty(e::Env, ::Val{:ball}) =
    getproperty(e, :ball), ∆ -> (Env(zero(e.arena), zero(e.player_a), zero(e.player_b), ∆), )

x::Player + y::Player = Player(x.length, x.position + y.position, x.velocity + y.velocity, 0)
x::Ball + y::Ball = Ball(x.position + y.position, x.velocity + y.velocity, x.radius)
x::Arena + y::Arena = Arena(x.dims + y.dims, x.margin + y.margin)
x::Env + y::Env = Env(x.arena + y.arena, x.player_a + y.player_a, x.player_b + y.player_b, x.ball + y.ball)

#=
function Zygote.accum(x::Player, y::Player)
    ans = Player(x.length, x.position + y.position, x.velocity, 0)
    println("$x   +    $y    =  $ans\n\n\n\n")
    return ans
end
function Zygote.accum(x::Ball, y::Ball)
    ans = Ball(x.position + y.position, x.velocity + y.velocity, x.radius)
    println("$x   +    $y    =  $ans\n\n\n\n")
    return ans
end

function Zygote.accum(x::Env, y::Env)
    ball = Ball(x.ball.position + y.ball.position, x.ball.velocity + y.ball.velocity,
                x.ball.radius)
    player_a = Player(x.player_a.length, x.player_a.position + y.player_a.position, x.player_a.velocity, 0)
    player_b = Player(x.player_b.length, x.player_b.position + y.player_b.position, x.player_b.velocity, 0)

    arena = Arena(x.arena.dims + y.arena.dims, x.arena.margin + y.arena.margin)

    ans = Env(arena, player_a, player_b, ball)
    println("$x   +    $y    =  $ans\n\n\n\n")
    return ans
end
Zygote.accum(x::Arena, y::Arena) = println("x = $x | y = $y\n\n\n\n")
=#

#=
@adjoint update_paddle_position(player::Player, arena::Arena, dir) =
    update_paddle_position(player, arena, dir), ∆ -> begin
        dummy = player.position.y + dir * player.velocity  # dummy holds the new position of the paddle
        # If the new position is within the acceptable range of motion, then velocity should be the gradient, wrt dir
        # Otherwise the gradient must be zero.
        0 ≤ dummy ≤ arena.dims.y - player.length ?
            (zero(player), zero(arena), ∆*Vector_(0f0, player.velocity)) :
            (zero(player), zero(arena), ∆*zero(dir))
        end

@adjoint function reset_ball!(arena::Arena, ball::Ball)
    reset_ball!(arena, ball), ∆ -> begin
        return (zero(arena), zero(ball))
    end
end

@adjoint function paddle_collision!(new_position::Vector_, player_a_pos::Vector_, player_b_pos::Vector_,
                                    paddle_length::AbstractFloat, ball::Ball, arena::Arena)
   paddle_collision(new_position, player_a_pos, player_b_pos, paddle_length, ball, arena), ∆ -> begin
    if new_position.x > arena.dims.x - ball.radius
        lower_bound_b  = max(player_b_pos.y - 3ball.radius/4f0, 0)
        higher_bound_b = min(player_b_pos.y + paddle_length + 3ball.radius/4f0, arena.dims.y)

        if lower_bound_b < new_position.y < higher_bound_b

        end
    else
        return (zero(new_position), zero(player_a_pos), zero(player_b_pos), zero(paddle_length),
                zero(ball), zero(arena))
    end
end
end

@nograd wall_collision!, get_obs#, paddle_collision!

Zygote.accum(x::Arena, y::Arena) = println("x = $x | y = $y")
Zygote.accum(x::Player, y::Player) = println("x = $x | y = $y")
=#


text:
using JuMP, EAGO, Gurobi

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[1, 2, 3, 4, 5, 6]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[5], 0.0)
JuMP.set_lower_bound(x[1], 0.0)
JuMP.set_lower_bound(x[4], 0.0)
JuMP.set_lower_bound(x[2], 0.0)
JuMP.set_lower_bound(x[6], 0.0)
JuMP.set_lower_bound(x[3], 0.0)
JuMP.set_upper_bound(x[1], 1.0)
JuMP.set_upper_bound(x[2], 2.0)
JuMP.set_upper_bound(x[3], 1.0)
JuMP.set_upper_bound(x[4], 4.0)
JuMP.set_upper_bound(x[5], 2.0)
JuMP.set_upper_bound(x[6], 6.0)

# ----- Objective ----- #
@NLobjective(m, Min, ( (2.5134-x[1]-x[3]-x[5])^2+ (2.044333373291-(exp(-0.05*x[2])*x[1]+exp(-0.05*x[4])*x[3]+exp(-0.05*x[6])*x[5]))^2+ (1.668404436564-(exp(-0.1*x[2])*x[1]+exp(-0.1*x[4])*x[3]+exp(-0.1*x[6])*x[5]))^2+ (1.366418021208-(exp(-0.15*x[2])*x[1]+exp(-0.15*x[4])*x[3]+exp(-0.15*x[6])*x[5]))^2+ (1.123232487372-(exp(-0.2*x[2])*x[1]+exp(-0.2*x[4])*x[3]+exp(-0.2*x[6])*x[5]))^2+ (0.9268897180037-(exp(-0.25*x[2])*x[1]+exp(-0.25*x[4])*x[3]+exp(-0.25*x[6])*x[5]))^2+ (0.7679338563728-(exp(-0.3*x[2])*x[1]+exp(-0.3*x[4])*x[3]+exp(-0.3*x[6])*x[5]))^2+ (0.6388775523106-(exp(-0.35*x[2])*x[1]+exp(-0.35*x[4])*x[3]+exp(-0.35*x[6])*x[5]))^2+ (0.5337835317402-(exp(-0.4*x[2])*x[1]+exp(-0.4*x[4])*x[3]+exp(-0.4*x[6])*x[5]))^2+ (0.4479363617347-(exp(-0.45*x[2])*x[1]+exp(-0.45*x[4])*x[3]+exp(-0.45*x[6])*x[5]))^2+ (0.377584788435-(exp(-0.5*x[2])*x[1]+exp(-0.5*x[4])*x[3]+exp(-0.5*x[6])*x[5]))^2+ (0.3197393199326-(exp(-0.55*x[2])*x[1]+exp(-0.55*x[4])*x[3]+exp(-0.55*x[6])*x[5]))^2+ (0.2720130773746-(exp(-0.6*x[2])*x[1]+exp(-0.6*x[4])*x[3]+exp(-0.6*x[6])*x[5]))^2+ (0.2324965529032-(exp(-0.65*x[2])*x[1]+exp(-0.65*x[4])*x[3]+exp(-0.65*x[6])*x[5]))^2+ (0.1996589546065-(exp(-0.7*x[2])*x[1]+exp(-0.7*x[4])*x[3]+exp(-0.7*x[6])*x[5]))^2+ (0.1722704126914-(exp(-0.75*x[2])*x[1]+exp(-0.75*x[4])*x[3]+exp(-0.75*x[6])*x[5]))^2+ (0.1493405660168-(exp(-0.8*x[2])*x[1]+exp(-0.8*x[4])*x[3]+exp(-0.8*x[6])*x[5]))^2+ (0.1300700206922-(exp(-0.85*x[2])*x[1]+exp(-0.85*x[4])*x[3]+exp(-0.85*x[6])*x[5]))^2+ (0.1138119324644-(exp(-0.9*x[2])*x[1]+exp(-0.9*x[4])*x[3]+exp(-0.9*x[6])*x[5]))^2+ (0.1000415587559-(exp(-0.95*x[2])*x[1]+exp(-0.95*x[4])*x[3]+exp(-0.95*x[6])*x[5]))^2+ (0.0883320908454-(exp(-x[2])*x[1]+exp(-x[4])*x[3]+exp(-x[6])*x[5]))^2+ (0.0783354401935-(exp(-1.05*x[2])*x[1]+exp(-1.05*x[4])*x[3]+exp(-1.05*x[6])*x[5]))^2+ (0.06976693743449-(exp(-1.1*x[2])*x[1]+exp(-1.1*x[4])*x[3]+exp(-1.1*x[6])*x[5]))^2+ (0.06239312536719-(exp(-1.15*x[2])*x[1]+exp(-1.15*x[4])*x[3]+exp(-1.15*x[6])*x[5]))^2))

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
println("run time: $run_time")


text:
module TestDiscrete
using CustomTest
using Base.Test
using ControlSystems

C_111 = ss([-5], [2], [3], [0])
C_212 = ss([-5 -3; 2 -9], [1; 2], eye(2), [0; 0])
C_221 = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0], [0 0])
C_222_d = ss([-5 -3; 2 -9], [1 0; 0 2], eye(2), eye(2))

macro test_c2d(ex, sys_sol, mat_sol)
    quote
        sys, mat = $ex
        @test_approx_eq sys $sys_sol
        @test_approx_eq mat $mat_sol
    end
end
@test c2d(ss(4*eye(2)), 0.5, :zoh) == (ss(4*eye(2), 0.5), zeros(0, 2))
@test c2d(ss(4*eye(2)), 0.5, :foh) == (ss(4*eye(2), 0.5), zeros(0, 2))
@test_c2d(c2d(C_111, 0.01, :zoh),
ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0])
@test_c2d(c2d(C_111, 0.01, :foh),
ss([0.951229424500714], [0.01902855227625244], [3], [0.029506188017136226], 0.01),
[1 -0.009835396005712075])
@test_c2d(c2d(C_212, 0.01, :zoh),
ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
[0.009466805409666932; 0.019219966830212765], eye(2), [0; 0], 0.01),
[1 0 0; 0 1 0])
@test_c2d(c2d(C_212, 0.01, :foh),
ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
[0.008957940478201584; 0.018468989584974498], eye(2), [0.004820885889482196;
0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675])
@test_c2d(c2d(C_221, 0.01, :zoh),
ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
[0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
[1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0])
@test_c2d(c2d(C_221, 0.01, :foh),
ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
[0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
[1 0], [0.004917457305816479 -9.657141633428213e-5], 0.01),
[1 0 -0.004917457305816479 9.657141633428213e-5;
0 1 -3.219047211142736e-5 -0.009706152723187249])
@test_c2d(c2d(C_222_d, 0.01, :zoh),
ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
[0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
eye(2), eye(2), 0.01), [1 0 0 0; 0 1 0 0])
@test_c2d(c2d(C_222_d, 0.01, :foh),
ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
[0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
eye(2), [1.0049174573058164 -9.657141633428213e-5; 3.219047211142736e-5 1.0097061527231872], 0.01),
[1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249])

# ERRORS
@test_throws ErrorException c2d(ss([1], [2], [3], [4], 0.01), 0.01)   # Already discrete
@test_throws ErrorException c2d(ss([1], [2], [3], [4], -1), 0.01)     # Already discrete

end


text:
# High-level plotting functions for use with `BlockMethod`s

"""
    showsample([backend], method, sample)

Show an unprocessed `sample` for `LearningMethod` `method` to
`backend::`[`ShowBackend`](#).

## Examples

```julia
data, blocks = loaddataset("imagenette2-160", (Image, Label))
method = ImageClassificationSingle(data)
sample = getobs(data, 1)
showsample(method, sample)  # select backend automatically
showsample(ShowText(), method, sample)
```
"""
function showsample(backend::ShowBackend, method::BlockMethod, sample)
    blocks = ("Input" => method.blocks[1], "Target" => method.blocks[2])
    showblock(backend, blocks, sample)
end
showsample(method::BlockMethod, sample) = showsample(default_showbackend(), method, sample)

"""
    showsample([backend], method, sample)

Show a vector of unprocessed `samples` for `LearningMethod` `method` to
`backend::`[`ShowBackend`](#).

## Examples

```julia
data, blocks = loaddataset("imagenette2-160", (Image, Label))
method = ImageClassificationSingle(data)
samples = [getobs(data, i) for i in 1:4]
showsamples(method, samples)  # select backend automatically
showsamples(ShowText(), method, samples)
```
"""
function showsamples(backend::ShowBackend, method::BlockMethod, samples)
    blocks = ("Input" => method.blocks[1], "Target" => method.blocks[2])
    showblocks(backend, blocks, samples)
end
showsamples(method::BlockMethod, samples) = showsamples(default_showbackend(), method, sample)

"""
    showencodedsample([backend], method, encsample)

Show an encoded sample `encsample` to `backend`.
"""
function showencodedsample(backend::ShowBackend, method::BlockMethod, encsample)
    xblock, yblock = encodedblock(method.encodings, method.blocks, true)
    showblockinterpretable(
        backend,
        method.encodings,
        ("x" => xblock, "y" => yblock),
        encsample)
end
showencodedsample(method, encsample) = showencodedsample(default_showbackend(), method, encsample)

"""
    showencodedsamples([backend], method, encsamples)

Show a vector of encoded samples `encsamples` to `backend`.
"""
function showencodedsamples(backend::ShowBackend, method::BlockMethod, encsamples::AbstractVector)
    xblock, yblock = encodedblock(method.encodings, method.blocks, true)
    showblocksinterpretable(
        backend,
        method.encodings,
        ("x" => xblock, "y" => yblock),
        encsamples)
end

"""
    showbatch([backend], method, batch)

Show a collated batch of encoded samples to `backend`.
"""
function showbatch(backend::ShowBackend, method::BlockMethod, batch)
    encsamples = collect(DataLoaders.obsslices(batch))
    showencodedsamples(backend, method, encsamples)
end
showbatch(method, batch) = showbatch(default_showbackend(), method, batch)

"""
    showprediction([backend], method, pred)
    showprediction([backend], method, sample, pred)

Show a prediction `pred`. If a `sample` is also given, show it next to
the prediction.
"""
function showprediction(backend::ShowBackend, method::BlockMethod, pred)
    predblock = decodedblock(method.encodings, method.outputblock, true)
    showblock(backend, "Prediction" => predblock, pred)
end

function showprediction(backend::ShowBackend, method::BlockMethod, sample, pred)
    predblock = decodedblock(method.encodings, method.outputblock, true)
    showblock(
        backend,
        ("Sample" => method.blocks, "Prediction" => predblock),
        (sample, pred)
    )
end

showprediction(method::BlockMethod, args...) =
    showprediction(default_showbackend(), method, args...)

"""
    showpredictions([backend], method, preds)
    showpredictions([backend], method, samples, preds)

Show predictions `pred`. If `samples` are also given, show them next to
the prediction.
"""
function showpredictions(backend::ShowBackend, method::BlockMethod, preds)
    predblock = decodedblock(method.encodings, method.outputblock, true)
    showblocks(backend, "Prediction" => predblock, preds)
end

function showpredictions(backend::ShowBackend, method::BlockMethod, samples, preds)
    predblock = decodedblock(method.encodings, method.outputblock, true)
    showblocks(
        backend,
        ("Sample" => method.blocks, "Prediction" => predblock),
        collect(zip(samples, preds)),
    )
end

showpredictions(method::BlockMethod, args...) =
    showpredictions(default_showbackend(), method, args...)

"""
    showoutput([backend], method, output)
    showoutput([backend], method, encsample, output)

Show a model output to `backend`. If an encoded sample `encsample` is also
given, show it next to the output.
"""
function showoutput(backend::ShowBackend, method::BlockMethod, output)
    showblockinterpretable(backend, method.encodings, "Output" => method.outputblock, output)
end
function showoutput(backend::ShowBackend, method::BlockMethod, encsample, output)
    encsampleblock = encodedblock(method.encodings, method.blocks, true)
    outblock = method.outputblock
    showblockinterpretable(
        backend,
        method.encodings,
        ("Encoded sample" => encsampleblock, "Output" => outblock),
        (encsample, output))
end
showoutput(method::BlockMethod, args...) = showoutput(default_showbackend(), method, args...)

"""
    showoutputs([backend], method, outputs)
    showoutputs([backend], method, encsamples, outputs)

Show model outputs to `backend`. If a vector of encoded samples `encsamples` is also
given, show them next to the outputs. Use [`showoutputbatch`](#) to show collated
batches of outputs.
"""
function showoutputs(backend::ShowBackend, method::BlockMethod, outputs)
    showblocks(backend, "Output" => method.outputblock, outputs)
end
function showoutputs(backend::ShowBackend, method::BlockMethod, encsamples, outputs)
    encsampleblock = encodedblock(method.encodings, method.blocks, true)
    outblock = method.outputblock
    showblocksinterpretable(
        backend,
        method.encodings,
        ("Encoded sample" => encsampleblock, "Output" => outblock),
        collect(zip(encsamples, outputs)))
end

showoutputs(method::BlockMethod, args...) = showoutputs(default_showbackend(), method, args...)

"""
    showoutputbatch([backend], method, outputbatch)
    showoutputbatch([backend], method, batch, outputbatch)

Show collated batch of outputs to `backend`. If a collated batch of encoded samples
`batch` is also given, show them next to the outputs. See [`showoutputs`](#) if you
have vectors of outputs and not collated batches.
"""
function showoutputbatch(backend::ShowBackend, method::BlockMethod, outputbatch)
    outputs = collect(DataLoaders.obsslices(outputbatch))
    return showoutputs(backend, method, outputs)
end
function showoutputbatch(backend::ShowBackend, method::BlockMethod, batch, outputbatch)
    encsamples = collect(DataLoaders.obsslices(batch))
    outputs = collect(DataLoaders.obsslices(outputbatch))
    return showoutputs(backend, method, encsamples, outputs)
end

showoutputbatch(method::BlockMethod, args...) = showoutputbatch(default_showbackend(), method, args...)

# Testing helper

"""
    test_method_show(method, backend::ShowBackend)

Test suite that tests that all learning method-related `show*` functions
work for `backend`

## Keyword arguments

- `sample = mockblock(method.blocks)`: Sample data to use for tests.
- `output = mockblock(method.outputblock)`: Model output data to use for tests.
"""
function test_method_show(
        method::LearningMethod, backend::ShowBackend;
        sample = mockblock(method.blocks),
        output = mockblock(method.outputblock),
        context = Training())

    encsample = encode(method, context, sample)
    pred = decodeŷ(method, context, output)

    Test.@testset "`show*` test suite for learning method $method and backend $(typeof(backend))" begin
        @test (showsample(backend, method, sample); true)
        @test (showencodedsample(backend, method, encsample); true)
        @test (showoutput(backend, method, output); true)
        @test (showoutput(backend, method, encsample, output); true)
        @test (showprediction(backend, method, pred); true)
        @test (showprediction(backend, method, sample, pred); true)
    end
end

# Deprecations


text:
"""
$(TYPEDSIGNATURES)

Gets the next prime `p`, starting from `n`, for which `(p + offset) % divisor ==
0` holds.

```jldoctest
julia> next_offset_divisible_prime(3, 1, 4, 1000)
3

julia> next_offset_divisible_prime(5, 1, 4, 1000)
7

julia> next_offset_divisible_prime(4, 1, 4, 1000)
7

```
"""
function next_offset_divisible_prime(n::Int, offset::Int, divisor::Int, tries::Int)
    for i = 1:tries
        prime = nextprime(n, i)

        if (prime + offset) % divisor == 0
            return prime
        end
    end

    error("There are no primes from $n, within $tries tries, " *
          "divisible by $divisor with offset $offset")
end

"""
$(TYPEDSIGNATURES)

To check if  a given design is  a Plackett-Burman design, we must check for the
following properties, obtained in the original Plackett-Burman paper:

1. Each component is replicated at each  of its values the same number of times,
   that is, the sum of elements in each column is zero

2. Each  pair of components  occur together at  every combination of  values the
   same number of times, that is, the sum of each pair of columns will produce a
   column with the  same number of occurrences  of ``2`` and ``-2``,  and twice that
   number of occurrences of ``0``

> Plackett,  R.L. and  Burman, J.P.,  1946. The  design of  optimum multifactorial
> experiments. Biometrika, 33(4), pp.305-325.

```jldoctest
julia> isplackettburman(plackettburman(2))
true

julia> isplackettburman(plackettburman(4))
true

julia> isplackettburman(plackettburman(16))
true

julia> isplackettburman(rand(Int, 4,4))
false

```
"""
function isplackettburman(d::Matrix{Int})
    sum(sum(d, dims = 1)) == 0 && sum(sum(d' * d, dims = 1) / size(d, 1)) == size(d, 2)
end

"""
$(TYPEDSIGNATURES)

The [Paley  construction](https://en.wikipedia.org/wiki/Paley_construction) is a
method for constructing Hadamard matrices using finite fields.

```jldoctest
julia> paley(Matrix{Int}(undef, 8, 8))
8×8 Matrix{Int64}:
 -1   1   1  -1   1   1   1  -1
  1   1  -1   1   1   1  -1  -1
  1  -1   1   1   1  -1  -1   1
 -1   1   1   1  -1  -1   1   1
  1   1   1  -1  -1   1   1  -1
  1   1  -1  -1   1   1  -1   1
  1  -1  -1   1   1  -1   1   1
 -1  -1   1   1  -1   1   1   1

```
"""
function paley(matrix::Matrix{Int})
    dimension::Tuple{Int, Int} = size(matrix)

    starting_row::Int = (dimension[1] == dimension[2] ? 1 : 2)
    ending_row::Int = dimension[1]

    prime::Int = dimension[2]

    residues_length::Int = 1 + floor((prime - 1) / 2)
    residues::Array{Int, 1} = Array{Int, 1}(undef, residues_length)

    residues[1] = 0
    residues[2] = 1

    @inbounds for i = 2:(residues_length - 1)
        residues[i + 1] = (i * i) % prime
    end

    sort!(residues)

    @inbounds for i = starting_row:ending_row
       @inbounds for j = 1:prime
           offset_value::Int = (i + j - 1) % prime
           if offset_value in residues
               matrix[i, j] = -1
           else
               matrix[i, j] = 1
           end
       end
    end

    return matrix
end

"""
$(TYPEDSIGNATURES)

Constructs a Plackett-Burman  design with size `matrix_size` if  possible, or to
the closest, largest, number for which it is possible.

```jldoctest
julia> plackettburman(4)
8×7 Matrix{Int64}:
  1   1   1   1   1   1   1
 -1   1  -1   1   1  -1  -1
  1  -1   1   1  -1  -1  -1
 -1   1   1  -1  -1  -1   1
  1   1  -1  -1  -1   1  -1
  1  -1  -1  -1   1  -1   1
 -1  -1  -1   1  -1   1   1
 -1  -1   1  -1   1   1  -1

```
"""
function plackettburman(matrix_size::Int)
    p = next_offset_divisible_prime(matrix_size, 1, 4, 1000)

    A = Matrix{Int}(undef, p + 1, p)
    A[1, :] = ones(p)
    paley(A)
end


text:
"""
    CategoricalLikelihood(l::AbstractLink=SoftMaxLink())

Categorical likelihood is to be used if we assume that the
uncertainity associated with the data follows a Categorical distribution.
```math
    p(y|f_1, f_2, \\dots, f_{n-1}) = Categorical(y | l(f_1, f_2, \\dots, f_{n-1}, 0))
```
Given an `AbstractVector` [f_1, f_2, ..., f_{n-1}], returns a `Categorical` distribution,
with probabilities given by `l(f_1, f_2, ..., f_{n-1}, 0)`.
"""
struct CategoricalLikelihood{Tl<:AbstractLink}
    invlink::Tl
end

CategoricalLikelihood() = CategoricalLikelihood(SoftMaxLink())

(l::CategoricalLikelihood)(f::AbstractVector{<:Real}) = Categorical(l.invlink(vcat(f, 0)))

function (l::CategoricalLikelihood)(fs::AbstractVector)
    return Product(Categorical.(l.invlink.(vcat.(fs, 0))))
end


text:
"""
DEMO: Generating RIRs for a microphone pair, for multiple source positions.
-----

- The source positions are selected on concentric circular grids around the array center. The grid can be uniformly spaced, or randomized.
- The RIRs and other metatata are exported as .npz file (which can then be loaded in Python scripts using numpy.load - an example is shown)

Author
------
Maja Taseska
FWO Postdoctoral Fellow, KU Leuven ESAT-STADIUS
"""

using ImageMethodReverb, Random, Plots, LinearAlgebra, NPZ
Fs = 16000                          # sampling frequency

#-------------------------
# Room parameters
#-----------------
L = Lx, Ly, Lz = 7.,8.,3.                   # room dimensions in meters
micdist = 0.15;                             # microphone distance
T60 = 0.4                                   # reverberation time in seconds
Nt = round(Int,T60*Fs);                     # RIR length

#-------------------------
# Sources
#-------------------------
d_grid = [1.5, 2, 2.5, 3];      # distance grid
theta_grid = collect(0:10:355)*pi/180;      # angle grid
ee_grid = hcat(cos.(theta_grid),sin.(theta_grid))';

# Store the randomized grid in advance
preliminary_dist = zeros(length(d_grid),length(theta_grid));
preliminary_angle = zeros(length(d_grid),length(theta_grid));
for dd = 1:length(d_grid)
    for th = 1:length(theta_grid)
        # OFF GRID FOR TESTING
        ddt = d_grid[dd] + (rand().- 0.5)*0.4
        tht = rand()*2*pi
        preliminary_dist[dd,th] = ddt
        preliminary_angle[dd,th] = tht
    end
end

#-------------------------
# Microphones
#-------------------------
C = [2.5; 2.3; 1.6];                        # array center
rec1 = C - [micdist/2, 0, 0];
rec2 = C + [micdist/2, 0, 0];
currentreceivers = [tuple(rec1...),tuple(rec2...)];

#--------------------------------
# MAIN PART: GENERATING THE RIRS
#-------------------------------

# Initialize arrays to store the different variables
sources_regulargrid = [];
distances_regulargrid = [];
rirs_regulargrid = [];

sources_randgrid = [];
distances_randgrid = [];
rirs_randgrid = [];

for dd = 1:length(d_grid)               #---- Loop over the distance grid
    for th = 1:length(theta_grid)       # ---- Loop over the source angles

        # source positions on a regular grid
        tmpsourcepos = vcat(d_grid[dd]*ee_grid[:,th] + C[1:2], C[3])
        if (0.3 < tmpsourcepos[1] < Lx-0.3 && 0.3< tmpsourcepos[2] < Ly-0.3 )  # check that minimum 30 cm from walls
            tmprir, = rim(tuple(tmpsourcepos...),currentreceivers,L,T60,Nt,Fs)
            if sources_regulargrid == []
                sources_regulargrid = tmpsourcepos
                distances_regulargrid = d_grid[dd]
                rirs_regulargrid = tmprir
            else
                sources_regulargrid = hcat(sources_regulargrid,tmpsourcepos)
                distances_regulargrid = vcat(distances_regulargrid,d_grid[dd])
                rirs_regulargrid = cat(dims = 3,rirs_regulargrid,tmprir)
            end

        end

        # source positions on a randomized grid
        ddt = preliminary_dist[dd,th]
        tht = preliminary_angle[dd,th]
        eet = hcat(cos.(tht),sin.(tht))';
        tmpt = vcat(ddt*eet + C[1:2], C[3])
        #
        if (0.3 < tmpt[1] < Lx-0.3 && 0.3< tmpt[2] < Ly-0.3 )     # check that minimum 30 cm from walls
            tmprir, = rim(tuple(tmpt...),currentreceivers,L,T60,Nt,Fs)
            if sources_randgrid == []
                sources_randgrid = tmpt
                distances_randgrid = ddt
                rirs_randgrid = tmprir
            else
                sources_randgrid = hcat(sources_randgrid,tmpt)
                distances_randgrid = vcat(distances_randgrid,ddt)
                rirs_randgrid = cat(dims = 3,rirs_randgrid,tmprir)
            end
        end
    end
end

#------------------
# Saving the data
#------------------
mydict = Dict("rirs_regulargrid"=>rirs_regulargrid,
              "rirs_randgrid"=>rirs_randgrid,
              "sources_regulargrid"=>sources_regulargrid,
              "sources_randgrid"=>sources_randgrid,
              "distances_regulargrid"=>distances_regulargrid,
              "distances_randgrid"=>distances_randgrid,
              "T60"=>[T60],
              "roomdim" => collect(L),
              "fs"=>[Fs],
              "micdist"=>micdist,
              "arrcenter"=> C);

npzwrite("saved_data.npz", mydict)

#---------------------------------------------
# Uncomment the two lines below to plot the source positions
#--------------------------------------------

#scatter!(sources_regulargrid[1,:],sources_regulargrid[2,:])
#scatter!(sources_randgrid[1,:],sources_randgrid[2,:])

#---------------------------------------------
#---------------------------------------------
#---------------------------------------------
# !!! ONLY TO BE USED IN PYTHON SCRIPTS
#---------------------------------------------
# The code excerpt below shows how to load the saved data in a python script.
# Create a python file, e.g., "check_positions.py", copy the code below, and run the file.
# Make sure that the imported packages (numpy and matplotlib.pyplot) are installed in your environment.
#---------------------------------------------
"""
import numpy as np
from matplotlib import pyplot as plt

data = np.load('saved_data.npz')

# To check what has been stored uncomment:
#data.files

s1 = data['sources_randgrid'].T
s2 = data['sources_regulargrid'].T
rir1 = data['rirs_randgrid']
rir2 = data['rirs_regulargrid']
arrcenter = data['arrcenter']
micx = [arrcenter[0]-data['micdist']/2,arrcenter[0]+data['micdist']/2]
micy = [arrcenter[1], arrcenter[1]]

plt.figure()
plt.scatter(s1[:,0],s1[:,1])
plt.scatter(s2[:,0],s2[:,1], color = 'r')
plt.scatter(micx,micy, color = 'k')

"""


text:
########################################################
# Utility functions used by the various code generators
########################################################

# Checks if the coefficient has constant value.
# If so, also returns the value.
function is_constant_coef(c)
    isit = false;
    val = 0;
    for i=1:length(coefficients)
        if c === coefficients[i].symbol
            isit = (typeof(coefficients[i].value[1]) <: Number);
            if isit
                val = coefficients[i].value[1];
            else
                val = coefficients[i].value[1].name;
            end
        end
    end

    return (isit, val);
end

# Checks the type of coefficient: constant, genfunction, or variable
# Returns: (type, val)
# special: type=-1, val=0
# number: type=0, val=number
# constant coef: type=1, val=number
# genfunction: type=2, val= index in genfunctions array
# variable: type=3, val=index in variables array
# indexed coefficient: type=4, val=index in coefficients array
function get_coef_val(c)
    if c.index == -1
        # It is a special symbol like dt
        return (-1, 0);
    elseif c.index == 0
        # It is a number
        return(0, c.name);
    end

    type = 0;
    val = 0;
    for i=1:length(coefficients)
        if c.name == string(coefficients[i].symbol)
            if typeof(c.index) == Int
                isit = (typeof(coefficients[i].value[c.index]) <: Number);
                if isit
                    type = 1; # a constant wrapped in a coefficient ... not ideal
                    val = coefficients[i].value[c.index];
                else
                    type = 2; # a function
                    name = coefficients[i].value[c.index].name;
                    for j=1:length(genfunctions)
                        if name == genfunctions[j].name
                            val = j;
                        end
                    end
                end
            else
                # indexed coefficient
                type = 4;
                val = i;
            end
        end
    end
    if type == 0 # a variable
        for i=1:length(variables)
            if c.name == string(variables[i].symbol)
                type = 3;
                val = variables[i].index;
            end
        end
    end
    if type == 0 # something special with a nonzero index
        return (-1,0);
    end

    return (type, val);
end

# Returns the index in finch's coefficient array
# or -1 if it is not there.
function get_coef_index(c)
    ind = -1;
    for i=1:length(coefficients)
        if c.name == string(coefficients[i].symbol)
            ind = coefficients[i].index;
        end
    end

    return ind;
end

# Makes a name for this value based on coefficient or variable name, derivative and other modifiers, and component index.
# c is a SymEntity
function make_entity_name(c)
    # Special case for dt or other special symbols
    if c.index == -1
        return c.name;
    end

    tag = "";
    type_label = "value_";
    for i=1:length(c.flags)
        tag = c.flags[i] * tag;
        if c.flags[i] == "NEIGHBORHOOD"
            tmp = type_label[1];
            type_label = "noBroadcast_";
        end
    end
    for i=1:length(c.derivs)
        tag = "D"*string(c.derivs[i]) * tag;
    end

    if typeof(c.index) == Int
        str = type_label*tag*"_"*string(c.name)*"_"*string(c.index);
    else
        str = type_label*tag*"_"*string(c.name)*"_";
        for i=1:length(c.index)
            str *= string(c.index[i]);
        end
    end

    return str;
end

function is_test_function(ent)
    for t in test_functions
        if string(t.symbol) == ent.name
            return true;
        end
    end
    return false;
end

function is_unknown_var(ent, vars)
    if typeof(vars) == Variable
        return ent.name == string(vars.symbol);
    elseif typeof(vars) <:Array
        for v in vars
            if string(v.symbol) == ent.name
                return true;
            end
        end
    end
    return false;
end

# They are the same if all parts of them are the same.
function is_same_entity(a, b)
    return a.name == b.name && a.index == b.index && a.derivs == b.derivs && a.flags == b.flags;
end

# Replaces math operators with broadcast(dot) versions
function broadcast_all_ops(ex)
    if typeof(ex) <: Array
        result = [];
        for i=1:length(ex)
            push!(result, broadcast_all_ops(ex[i]));
        end
        return result;

    elseif typeof(ex) == Expr
        if ex.head === :call
            bopped = false;
            if ex.args[1] === :^ ex.args[1] = :.^
            elseif ex.args[1] === :/ ex.args[1] = :./
            elseif ex.args[1] === :* ex.args[1] = :.*
            elseif ex.args[1] === :+ ex.args[1] = :.+
            elseif ex.args[1] === :- ex.args[1] = :.-
            elseif ex.args[1] in [:.+, :.-, :.*, :./, :.^] # no change
            else # named operator: abs -> abs.
                # Make sure the arguments are suitable for broadcasting.
                suitable = true;
                for i=2:length(ex.args)
                    if typeof(ex.args[i]) == SymEntity && "NEIGHBORHOOD" in ex.args[i].flags
                        # Don't broadcast
                        suitable = false;
                    end
                end
                if suitable
                    argtuple = Expr(:tuple, ex.args[2]);
                    if length(ex.args) > 2
                        append!(argtuple.args, ex.args[3:end]);
                    end
                    bop = Expr(:., ex.args[1], argtuple);
                    ex = bop;
                    bopped = true;
                end
            end
            if bopped
                for i=1:length(ex.args[2].args)
                    ex.args[2].args[i] = broadcast_all_ops(ex.args[2].args[i]);
                end
            else
                for i=2:length(ex.args)
                    ex.args[i] = broadcast_all_ops(ex.args[i]);
                end
            end
        end

        return ex;

    elseif typeof(ex) == Symbol
        # replace ^,/,+,- with .^,./,.+,.-
        if ex === :^ ex = :.^
        elseif ex === :/ ex = :./
        elseif ex === :* ex = :.*
        elseif ex === :+ ex = :.+
        elseif ex === :- ex = :.-
        end
    else
        return ex;
    end
end

# Turn entities with index=[i,j] into u[index_val_i, index_val_j]
function insert_indices(ex)
    if typeof(ex) <: Array
        result = [];
        for i=1:length(ex)
            push!(result, insert_indices(ex[i]));
        end
        return result;

    elseif typeof(ex) == SymEntity
        if typeof(ex.index) <: Array
            indices = copy(ex.index);

            for i=1:length(indices)
                indices[i] = Symbol("index_val_"*string(indices[i]));
            end
            ex = Expr(:ref, indexed_thing);
            append!(ex.args, indices);

        else
            for i=1:length(ex.args)
                ex.args[i] = insert_indices(ex.args[i]);
            end
        end

    end

    return ex;
end

# symex could be an array. If so, make a similar array containing an array of terms.
# This does any final processing of the terms. At this point symex should be a single term
# without a top level + or -. for example, a*b/c or a^2*b but not a*b+c
# What processing is done:
#   a/b -> a*(1/b)
#   Broadcast all ops: +-*/^ -> .+.-.*./.^
#   If it's a SymExpression, return the proccessed Expr in symex.tree
function process_terms(symex)
    if typeof(symex) <: Array
        terms = [];
        for i=1:length(symex)
            push!(terms, process_terms(symex[i]));
        end

    else
        if typeof(symex) == SymExpression
            newex = copy(symex.tree);
        elseif typeof(symex) == Expr
            newex = copy(symex);
        else
            newex = symex;
        end

        #println(string(newex) * " : " * string(typeof(newex)))

        if typeof(newex) == Expr && newex.head === :call
            if (newex.args[1] === :+ || newex.args[1] === :.+ ||
                (newex.args[1] === :- || newex.args[1] === :.-) && length(newex.args) > 2 )
                println("unexpected term in process terms: "*string(newex));
            else
                if newex.args[1] === :/ || newex.args[1] === :./
                    # change a/b to a*1/b
                    mulex = :(a.*b);
                    mulex.args[2] = newex.args[2];
                    newex.args[2] = 1;
                    mulex.args[3] = newex;
                    newex = mulex;
                end
            end

            # insert indices indexing_operator(u, i, j) -> u[index_val_i, index_val_j]
            newex = insert_indices(newex);

            # broadcast ops
            newex = broadcast_all_ops(newex);

            terms = newex;
        else
            # There is just one factor
            terms = newex;
        end
    end

    return terms;
end

# Assume the term looks like a*b*c or a*(b*c) or something similar.
function separate_factors(ex, var=nothing)
    test_part = nothing;
    trial_part = nothing;
    coef_part = nothing;
    test_ind = 0;
    trial_ind = 0;

    if typeof(ex) == Expr
        if (ex.args[1] === :.- || ex.args[1] === :-) && length(ex.args)==2
            # a negative sign. Stick it on one of the factors
            (test_part, trial_part, coef_part, test_ind, trial_ind) = separate_factors(ex.args[2], var);
            negex = :(-a);
            if !(coef_part === nothing)
                negex.args[2] = coef_part;
                coef_part = negex;
            elseif !(test_part === nothing)
                negex.args[2] = test_part;
                test_part = negex;
            elseif !(trial_part === nothing)
                negex.args[2] = trial_part;
                trial_part = negex;
            end

        elseif ex.args[1] === :.* || ex.args[1] === :*
            tmpcoef = [];
            # Recursively separate each arg to handle a*(b*(c*...))
            for i=2:length(ex.args)
                (testi, triali, coefi, testindi, trialindi) = separate_factors(ex.args[i], var);
                if !(testi === nothing)
                    test_part = testi;
                    test_ind = testindi;
                end
                if !(triali === nothing)
                    trial_part = triali;
                    trial_ind = trialindi;
                end
                if !(coefi === nothing)
                    push!(tmpcoef, coefi);
                end
            end
            if length(tmpcoef) > 1
                coef_part = :(a.*b);
                coef_part.args = [:.*];
                append!(coef_part.args, tmpcoef);
            elseif length(tmpcoef) == 1
                coef_part = tmpcoef[1];
            end

        else # It is an Expr, but not -() or * (note: at this point a/b is a*(1/b) )
            # This simplification may cause trouble somewhere. Revisit if needed.
            coef_part = ex;
        end
    elseif typeof(ex) == SymEntity
        numind = ex.index;
        if typeof(numind) <: Array
            numind = 1;
        end
        if is_test_function(ex)
            test_part = ex;
            test_ind = numind;
        elseif !(var === nothing) && is_unknown_var(ex, var)
            trial_part = ex;
            # find the trial index
            if typeof(var) <: Array
                tmpind = 0;
                for vi=1:length(var)
                    if is_unknown_var(ex,var[vi])
                        trial_ind = tmpind + numind;
                    else
                        tmpind += length(var[vi].symvar);
                    end
                end
            else
                trial_ind = numind;
            end
        else
            coef_part = ex;
        end
    else # a number?
        coef_part = ex;
    end

    # println(ex)
    # println(test_part)
    # println(trial_part)
    # println(coef_part)
    # println("")

    return (test_part, trial_part, coef_part, test_ind, trial_ind);
end

function replace_entities_with_symbols(ex; index=nothing)
    if typeof(ex) <: Array
        for i=1:length(ex)
            ex[i] = replace_entities_with_symbols(ex[i], index=index);
        end
    elseif typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = replace_entities_with_symbols(ex.args[i], index=index);
        end
    elseif typeof(ex) == SymEntity
        # use make_entity_name
        name = make_entity_name(ex)
        if !(index===nothing) && !(name == "dt")
            idx = Symbol(index);
            newex = :(a[$idx]);
            newex.args[1] = Symbol(name);
            return newex;
        end
        return Symbol(name);
    end

    return ex;
end

# Searches for any of the specified ents(strings matching entity.name) and adds a flag to entity.flags
function apply_flag_to_entities(ex, ents, flag; copy_ent=false)
    if typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = apply_flag_to_entities(ex.args[i], ents, flag, copy_ent=copy_ent);
        end

    elseif typeof(ex) == SymEntity
        if ex.name in ents
            if copy_ent
                newex = copy(ex);
                push!(newex.flags, flag);
                ex = newex;
            else
                push!(ex.flags, flag);
            end
        end
    end

    return ex;
end

# Apply a flag to any variable type entities only
function add_flag_to_var_entities(ex, vars, flag; nevermind="OHNEVERMIND", copy_ent=false)
    if typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = add_flag_to_var_entities(ex.args[i], vars, flag, nevermind=nevermind, copy_ent=copy_ent);
        end

    elseif typeof(ex) == SymEntity
        for fi=1:length(ex.flags)
            if occursin(nevermind, ex.flags[fi])
                return ex;
            end
        end
        for vi=1:length(vars)
            if ex.name == string(vars[vi].symbol)
                if copy_ent
                    newex = copy(ex);
                    push!(newex.flags, flag);
                    ex = newex;
                else
                    push!(ex.flags, flag);
                end
            end
        end
    end

    return ex;
end

# Traverses the Expr replacing a specific symbol
function replace_specific_symbol(ex, old_symbol, new_symbol)
    if typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = replace_specific_symbol(ex.args[i], old_symbol, new_symbol);
        end

    elseif typeof(ex) == Symbol
        if ex === old_symbol
            return new_symbol;
        end
    end

    return ex;
end

# Turns the generated code string into an Expr block to be put in a generated function.
function code_string_to_expr(s)
    lines = split(s, "\n", keepempty=false);

    # Remove all comments and empty lines
    # Wrap the whole thing in a block so the parser will make it a quote block
    clean_string = "begin\n";
    for i=1:length(lines)
        tmp = split(lines[i], "#")[1];
        if length(tmp) > 0
            clean_string *= tmp*"\n";
        end
    end
    clean_string *= "end";

    return Meta.parse(clean_string);
end

# Which side of a face is this entity on?
# Return a number: 0 = no side specified, 1/2 = side 1/2 in face2element order, 3 = average both, 4 = neighborhood
function get_face_side_info(ent)
    side = 0; # 0 means no side flag
    if typeof(ent) == SymEntity
        for flagi=1:length(ent.flags)
            if occursin("DGSIDE1", ent.flags[flagi]) || occursin("CELL1", ent.flags[flagi])
                return 1;
            elseif occursin("DGSIDE2", ent.flags[flagi]) || occursin("CELL2", ent.flags[flagi])
                return 2;
            elseif occursin("CENTRAL", ent.flags[flagi])
                return 3;
            elseif occursin("NEIGHBORHOOD", ent.flags[flagi])
                return 4;
            end
        end

    elseif typeof(ent) == Expr && length(ent.args) == 2 # This should only be allowed for negative signs :(-ent)
        side = get_face_side_info(ent.args[2]);

    # else number or ?? -> 0
    end

    return side;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
## draw: unify rand, sample

draw!(x::Distribution, o) = rand!(x, o)

draw!(x::AbstractVector, o::AbstractVector) = sample!(x, o)

function draw!(x::AbstractMatrix, o::AbstractMatrix)
    idx = sample(1:nobs(x), nobs(o))
    for (to, from) in enumerate(idx)
        o[to,:] = x[from,:]
    end
    return o
end

function draw!(x::AbstractDataFrame, o::AbstractDataFrame)
    idx = sample(1:nobs(x), nobs(o))
    for column in names(x)
        o[!,column] = x[idx,column]
    end
    return o
end

pick(x::AbstractVector, i::AbstractVector) = x[i]

pick(x::AbstractArray, i::AbstractVector) = copy(selectdim(x, 1, i))

pick(x::AbstractDataFrame, i::AbstractVector) = x[i,:]

"""
    MaximumEntropyCache{T<:Real}

A cache type for storing precomputed values for a particular input array.
This is intended to minimize memory allocations and time when drawing random samples.

# Fields
- `n::Int`: Number of elements in the input array.
- `t::Type{T}`: `eltype` of the input elements.
- `inds::Vector{Int}`: Indicies of the sorted data.
- `values::Vector{T}`: Sorted values of the input data.
- `Z::Vector{T}`: Intermediate points for the ordered values.
- `mtrm::T`: Trimmed mean of deviations.
- `m::Vector{T}`: Mean of the maximum entropy density within each interval.
- `U::Vector{T}`: Preallocated array for the random noise values.
- `quantiles::Vector{T}`: Preallocated array for sample quantiles.
- `v::Vector{T}`
"""
mutable struct MaximumEntropyCache{T <: Real}
    n::Int
    inds::Vector{Int}
    vals::Vector{T}
    Z::Vector{T}
    mtrm::T
    m::Vector{T}
    U::Vector{T}
    quantiles::Vector{T}
    v::Vector{T}
end

function MaximumEntropyCache()
    MaximumEntropyCache(0,
        Int[],
        Float64[],
        Float64[],
        0.0,
        Float64[],
        Float64[],
        Float64[],
        Float64[],
    )
end

function init!(cache::MaximumEntropyCache{T}, x::AbstractArray) where T
    cache.n = length(x)
    sorted!(cache, x)
    trimmed!(cache, x)
    intermediates!(cache)
    med!(cache)
    cache.U = zeros(T, cache.n)
    cache.quantiles = zeros(T, cache.n)
    cache.v = [y / cache.n for y in 0:cache.n]
    return nothing
end

"""
    sorted!(c::MaximumEntropyCache, x::AbstractArray)

Sets the sorted indices and values for `x`.
"""
function sorted!(c::MaximumEntropyCache, x::AbstractArray)
    c.inds = sortperm(x)
    c.vals = x[c.inds]
end

"""
    trimmed!(c::MaximumEntropyCache, x::AbstractArray)

Compute our trimmed mean of deviations.
"""
function trimmed!(c::MaximumEntropyCache, x::AbstractArray)
    c.mtrm = mean([abs(x[i] - x[i - 1]) for i in 2:c.n])
end

"""
    intermediates!(c::MaximumEntropyCache)

Compute our intermediate points for the ordered values.
"""
function intermediates!(c::MaximumEntropyCache{T}) where T
    c.Z = zeros(T, c.n + 1)
    for i in 2:c.n
        c.Z[i] = (c.vals[i - 1] + c.vals[i]) /  2
    end

    # Insert our lower and upper tails using our trimmed mean
    c.Z[1] = (c.vals[1] - c.mtrm)
    c.Z[end] = (c.vals[end] + c.mtrm)
end

"""
    med!(c::MaximumEntropyCache)

Compute the mean of the maximum entropy density within each interval such that
the ‘mean-preserving constraint’ is satisfied.
"""
function med!(c::MaximumEntropyCache{T}) where T
    c.m = zeros(T, c.n)
    c.m[1] = 0.75 * c.vals[1] + 0.25 * c.vals[1]

    for k in 2:(c.n - 1)
        c.m[k] = 0.25 * c.vals[k - 1] + 0.5 * c.vals[k] + 0.25 * c.vals[k + 1]
    end

    c.m[end] = 0.25 * c.vals[end - 1] + 0.75 * c.vals[end]
end

function draw!(cache::MaximumEntropyCache, x::T, o::T) where T <: AbstractArray
    # Generate random numbers from the [0, 1] uniform interval.
    sort!(rand!(cache.U))

    # Compute sample quantiles of the ME density at those points and sort them.
    for k in 1:cache.n
        ind = argmin(abs.(cache.v .- cache.U[k]))

        if cache.v[ind] > cache.U[k]
            ind -= 1
        end

        c = (2 * cache.m[ind] - cache.Z[ind] - cache.Z[ind + 1]) / 2
        y0 = cache.Z[ind] + c
        y1 = cache.Z[ind + 1] + c
        cache.quantiles[k] =
            y0 + (cache.U[k] - cache.v[ind]) * (y1 - y0) / (cache.v[ind + 1] - cache.v[ind])
    end

    # Reorder the sorted sample quantiles by using the ordering index.
    # This recovers the time dependence relationships of the originally observed data.
    recovered = sortperm(cache.inds)
    for i in 1:cache.n
        idx = recovered[i]
        o[i] = cache.quantiles[idx]
    end

    return o
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

initial_condition = initial_condition_convergence_test

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (0.0, 0.0, 0.0)
coordinates_max = (2.0, 2.0, 2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=10_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms_convergence_test)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)

###############################################################################
# run the simulation

sol = solve(ode, BS3(),
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-${JULIA} --color=yes --startup-file=no}"

thisdir="$(dirname "$(readlink -f "${BASH_SOURCE[0]}")")"
rootdir="$(dirname "$thisdir")"
export JULIA_PROJECT="$thisdir"
export JULIA_LOAD_PATH="@"

set -ex
${JULIA_CMD} \
    -e 'using Pkg; Pkg.develop(name="InteractiveCodeSearch", url=ARGS[1])' \
    "$rootdir"

exec ${JULIA_CMD} "${BASH_SOURCE[0]}" "$@"
=#

header = """

[![Build Status][ci-img]][ci-url]
[![codecov.io][codecov-img]][codecov-url]

![gif animation](search.gif "Searching code using @search")
"""

footer = """
[ci-img]: https://github.com/tkf/InteractiveCodeSearch.jl/actions/workflows/test.yml/badge.svg
[ci-url]: https://github.com/tkf/InteractiveCodeSearch.jl/actions/workflows/test.yml
[codecov-img]: http://codecov.io/github/tkf/InteractiveCodeSearch.jl/coverage.svg?branch=master
[codecov-url]: http://codecov.io/github/tkf/InteractiveCodeSearch.jl?branch=master

<!--
Generated by `./scripts/generate_readme.jl README.md`.
-->
"""

using Base.Docs: doc
import Markdown

using Documenter.Writers.MarkdownWriter: dropheaders

import InteractiveCodeSearch

function generate_readme(io::IO = stdout)
    text = sprint(show, "text/markdown", @doc InteractiveCodeSearch)
    lines = split(text, "\n")

    println(io, lines[1])
    println(io, header)
    for line in lines[2:end]
        println(io, line)
    end

    println(io, "## Reference")
    println(io)

    function loc(f)
        m1, = methods(f).ms
        return (m1.file, m1.line)
    end
    exports = [
        (name, getproperty(InteractiveCodeSearch, name))
        for name in names(InteractiveCodeSearch)
    ]
    exports = filter((x -> x[2] isa Function), exports)
    exports = sort(collect(exports), by=x -> loc(x[2]))

    for (name, exported) in exports
        md = dropheaders(doc(exported))
        println(io, "### `$name`")
        println(io)
        show(io, "text/markdown", md)
        println(io)
        println(io)
    end

    println(io, "### `InteractiveCodeSearch.CONFIG`")
    show(io, "text/markdown", dropheaders(@doc InteractiveCodeSearch.CONFIG))
    println(io)
    println(io)

    print(io, footer)
end

function generate_readme(filename::AbstractString)
    open(filename, "w") do io
        generate_readme(io)
    end
end

function rerender()
    io = IOBuffer()
    generate_readme(io)
    seek(io, 0)
    return Markdown.parse(io)
end

if isinteractive()
    # so that it can be called via `include("scripts/generate_readme.jl")()`
    generate_readme
else
    generate_readme(ARGS...)
end


text:
const ScreenID = UInt8
const ZIndex = Int
# ID, Area, clear, is visible, background color
const ScreenArea = Tuple{ScreenID, Scene}

abstract type GLScreen <: AbstractScreen end

mutable struct Screen <: GLScreen
    glscreen::GLFW.Window
    framebuffer::GLFramebuffer
    rendertask::RefValue{Task}
    screen2scene::Dict{WeakRef, ScreenID}
    screens::Vector{ScreenArea}
    renderlist::Vector{Tuple{ZIndex, ScreenID, RenderObject}}
    cache::Dict{UInt64, RenderObject}
    cache2plot::Dict{UInt16, AbstractPlot}
    framecache::Tuple{Matrix{RGB{N0f8}}, Matrix{RGB{N0f8}}}
    displayed_scene::Union{Scene, Nothing}
    function Screen(
            glscreen::GLFW.Window,
            framebuffer::GLFramebuffer,
            rendertask::RefValue{Task},
            screen2scene::Dict{WeakRef, ScreenID},
            screens::Vector{ScreenArea},
            renderlist::Vector{Tuple{ZIndex, ScreenID, RenderObject}},
            cache::Dict{UInt64, RenderObject},
            cache2plot::Dict{UInt16, AbstractPlot},
        )
        s = size(framebuffer)
        obj = new(
            glscreen, framebuffer, rendertask, screen2scene,
            screens, renderlist, cache, cache2plot,
            (Matrix{RGB{N0f8}}(undef, s), Matrix{RGB{N0f8}}(undef, reverse(s))),
            nothing
        )
        finalizer(obj) do obj
            # save_print("Freeing screen")
            empty!.((obj.renderlist, obj.screens, obj.cache, obj.screen2scene, obj.cache2plot))
            return
        end
        obj
    end
end
GeometryTypes.widths(x::Screen) = size(x.framebuffer.color)

Base.wait(x::Screen) = isassigned(x.rendertask) && wait(x.rendertask[])
Base.wait(scene::Scene) = wait(AbstractPlotting.getscreen(scene))
Base.show(io::IO, screen::Screen) = print(io, "GLMakie.Screen(...)")
Base.size(x::Screen) = size(x.framebuffer)

function insertplots!(screen::GLScreen, scene::Scene)
    get!(screen.screen2scene, WeakRef(scene)) do
        id = length(screen.screens) + 1
        push!(screen.screens, (id, scene))
        id
    end
    for elem in scene.plots
        insert!(screen, scene, elem)
    end
    foreach(s-> insertplots!(screen, s), scene.children)
end

function Base.delete!(screen::Screen, scene::Scene, plot::AbstractPlot)
    renderobject = get(screen.cache, objectid(plot)) do
        error("Could not find $(typeof(plot)) in current GLMakie screen!")
    end
    filter!(x-> x[3] !== renderobject, screen.renderlist)
end

function Base.empty!(screen::GLScreen)
    empty!(screen.renderlist)
    empty!(screen.screen2scene)
    empty!(screen.screens)
end

function destroy!(screen::Screen)
    empty!(screen)
    empty!(screen.cache)
    empty!(screen.cache2plot)
    destroy!(screen.glscreen)
end

function Base.resize!(window::GLFW.Window, resolution...)
    if isopen(window)
        retina_scale = retina_scaling_factor(window)
        w, h = resolution ./ retina_scale
        GLFW.SetWindowSize(window, round(Int, w), round(Int, h))
    end
end

function Base.resize!(screen::Screen, w, h)
    nw = to_native(screen)
    resize!(nw, w, h)
    fb = screen.framebuffer
    resize!(fb, (w, h))
end

function AbstractPlotting.backend_display(screen::Screen, scene::Scene)
    empty!(screen)
    register_callbacks(scene, to_native(screen))
    GLFW.PollEvents()
    insertplots!(screen, scene)
    GLFW.PollEvents()
    screen.displayed_scene = scene
    return
end

function to_jl_layout!(A, B)
    ind1, ind2 = axes(A)
    n = first(ind2) + last(ind2)
    for i in ind1
        @simd for j in ind2
            @inbounds c = A[i, j]
            c = mapc(c) do channel
                x = clamp(channel, 0.0, 1.0)
                return ifelse(isfinite(x), x, 0.0)
            end
            @inbounds B[n-j, i] = c
        end
    end
    return B
end

function fast_color_data!(dest::Array{RGB{N0f8}, 2}, source::Texture{T, 2}) where T
    GLAbstraction.bind(source)
    glPixelStorei(GL_PACK_ALIGNMENT, 1)
    glGetTexImage(source.texturetype, 0, GL_RGB, GL_UNSIGNED_BYTE, dest)
    GLAbstraction.bind(source, 0)
    nothing
end

function AbstractPlotting.colorbuffer(screen::Screen)
    if isopen(screen)
        GLFW.PollEvents()
        render_frame(screen) # let it render
        GLFW.SwapBuffers(to_native(screen))
        glFinish() # block until opengl is done rendering
        ctex = screen.framebuffer.color
        if size(ctex) != size(screen.framecache[1])
            s = size(ctex)
            screen.framecache = (Matrix{RGB{N0f8}}(undef, s), Matrix{RGB{N0f8}}(undef, reverse(s)))
        end
        fast_color_data!(screen.framecache[1], ctex)
        to_jl_layout!(screen.framecache...)
        return screen.framecache[2]
    else
        error("Screen not open!")
    end
end

Base.isopen(x::Screen) = isopen(x.glscreen)
function Base.push!(screen::GLScreen, scene::Scene, robj)
    filter!(screen.screen2scene) do (k, v)
        k.value != nothing
    end
    screenid = get!(screen.screen2scene, WeakRef(scene)) do
        id = length(screen.screens) + 1
        push!(screen.screens, (id, scene))
        id
    end
    push!(screen.renderlist, (0, screenid, robj))
    return robj
end

to_native(x::Screen) = x.glscreen
const gl_screens = GLFW.Window[]

"""
OpenGL shares all data containers between shared contexts, but not vertexarrays -.-
So to share a robjs between a context, we need to rewrap the vertexarray into a new one for that
specific context.
"""
function rewrap(robj::RenderObject{Pre}) where Pre
    RenderObject{Pre}(
        robj.main,
        robj.uniforms,
        GLVertexArray(robj.vertexarray),
        robj.prerenderfunction,
        robj.postrenderfunction,
        robj.boundingbox,
    )
end

const _global_gl_screen = Ref{Screen}()

# will get overloaded later
function renderloop end

# TODO a global is not very nice, but it's the simplest way right now to swap out
# the rendering loop
const opengl_renderloop = Ref{Function}(renderloop)

"""
Julia 1.0.3 doesn't have I:J, so we copy the implementation from 1.1 under a new name:
"""
function irange(I::CartesianIndex{N}, J::CartesianIndex{N}) where N
    CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
end

"""
Loads the makie loading icon and embedds it in an image the size of resolution
"""
function get_loading_image(resolution)
    icon = Matrix{N0f8}(undef, 192, 192)
    open(GLMakie.assetpath("loading.bin")) do io
        read!(io, icon)
    end
    img = zeros(RGBA{N0f8}, resolution...)
    center = resolution .÷ 2
    center_icon = size(icon) .÷ 2
    start = CartesianIndex(max.(center .- center_icon, 1))
    I1 = CartesianIndex(1, 1)
    stop = min(start + CartesianIndex(size(icon)) - I1, CartesianIndex(resolution))
    for idx in irange(start, stop)
        gray = icon[idx - start + I1]
        img[idx] = RGBA{N0f8}(gray, gray, gray, 1.0)
    end
    return img
end

function display_loading_image(screen::Screen)
    fb = screen.framebuffer
    fbsize = size(fb.color)
    image = get_loading_image(fbsize)
    if size(image) == fbsize
        nw = to_native(screen)
        fb.color[1:size(image, 1), 1:size(image, 2)] = image # transfer loading image to gpu framebuffer
        GLAbstraction.is_context_active(nw) || return
        w, h = fbsize
        glBindFramebuffer(GL_FRAMEBUFFER, 0) # transfer back to window
        glViewport(0, 0, w, h)
        glClearColor(0, 0, 0, 0)
        glClear(GL_COLOR_BUFFER_BIT)
        GLAbstraction.render(fb.postprocess[3]) # copy postprocess
        GLFW.SwapBuffers(nw)
    else
        error("loading_image needs to be Matrix{RGBA{N0f8}} with size(loading_image) == resolution")
    end

end

function Screen(;
        resolution = (10, 10), visible = false, title = "Makie",
        kw_args...
    )
    if !isempty(gl_screens)
        for elem in gl_screens
            isopen(elem) && destroy!(elem)
        end
        empty!(gl_screens)
    end

    window = GLFW.Window(
        name = title, resolution = (10, 10), # 10, because smaller sizes seem to error on some platforms
        windowhints = [
            (GLFW.SAMPLES,      0),
            (GLFW.DEPTH_BITS,   0),

            # SETTING THE ALPHA BIT IS REALLY IMPORTANT ON OSX, SINCE IT WILL JUST KEEP SHOWING A BLACK SCREEN
            # WITHOUT ANY ERROR -.-
            (GLFW.ALPHA_BITS,   8),
            (GLFW.RED_BITS,     8),
            (GLFW.GREEN_BITS,   8),
            (GLFW.BLUE_BITS,    8),

            (GLFW.STENCIL_BITS, 0),
            (GLFW.AUX_BUFFERS,  0),
            # (GLFW.RESIZABLE, GL_TRUE)
        ],
        visible = false,
        kw_args...
    )
    GLFW.SetWindowIcon(window , AbstractPlotting.icon())

    # tell GLAbstraction that we created a new context.
    # This is important for resource tracking, and only needed for the first context
    GLAbstraction.switch_context!(window)
    GLAbstraction.empty_shader_cache!()
    push!(gl_screens, window)

    GLFW.SwapInterval(0)

    # Retina screens on osx have a different scaling!
    retina_scale = retina_scaling_factor(window)
    resolution = round.(Int, retina_scale .* resolution)
    # Set the resolution for real now!
    GLFW.SetWindowSize(window, resolution...)
    fb = GLFramebuffer(Int.(resolution))

    screen = Screen(
        window, fb,
        RefValue{Task}(),
        Dict{WeakRef, ScreenID}(),
        ScreenArea[],
        Tuple{ZIndex, ScreenID, RenderObject}[],
        Dict{UInt64, RenderObject}(),
        Dict{UInt16, AbstractPlot}(),
    )

    GLFW.SetWindowRefreshCallback(window, window -> begin
        render_frame(screen)
        GLFW.SwapBuffers(window)
    end)

    if visible
        GLFW.ShowWindow(window)
    else
        GLFW.HideWindow(window)
    end
    screen.rendertask[] = @async((opengl_renderloop[])(screen))
    screen
end

function global_gl_screen()
    screen = if isassigned(_global_gl_screen) && isopen(_global_gl_screen[])
        _global_gl_screen[]
    else
        _global_gl_screen[] = Screen()
        _global_gl_screen[]
    end
    return screen
end

function global_gl_screen(resolution::Tuple, visibility::Bool, tries = 1)
    # ugly but easy way to find out if we create new screen.
    # could just be returned by global_gl_screen, but dont want to change the API
    isold = isassigned(_global_gl_screen) && isopen(_global_gl_screen[])
    screen = global_gl_screen()
    GLFW.set_visibility!(to_native(screen), visibility)
    resize!(screen, resolution...)
    new_size = GLFW.GetWindowSize(to_native(screen))
    # I'm not 100% sure, if there are platforms where I'm never
    # able to resize the screen (opengl might just allow that).
    # so, we guard against that with just trying another resize one time!
    if ((new_size.width, new_size.height) != resolution) && tries == 1
        # resize failed. This may happen when screen was previously
        # enlarged to fill screen. WE NEED TO DESTROY!! (I think)
        destroy!(screen)
        # try again
        return global_gl_screen(resolution, visibility, 2)
    end
    # show loading image on fresh screen
    isold || display_loading_image(screen)
    screen
end

function pick_native(screen::Screen, xy::Vec{2, Float64})
    isopen(screen) || return SelectionID{Int}(0, 0)
    sid = Base.RefValue{SelectionID{UInt16}}()
    window_size = widths(screen)
    fb = screen.framebuffer
    buff = fb.objectid
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1])
    glReadBuffer(GL_COLOR_ATTACHMENT1)
    x, y = floor.(Int, xy)
    w, h = window_size
    if x > 0 && y > 0 && x <= w && y <= h
        glReadPixels(x, y, 1, 1, buff.format, buff.pixeltype, sid)
        return convert(SelectionID{Int}, sid[])
    end
    return SelectionID{Int}(0, 0)
end

function AbstractPlotting.pick(scene::SceneLike, screen::Screen, xy::Vec{2, Float64})
    sid = pick_native(screen, xy)
    if haskey(screen.cache2plot, sid.id)
        plot = screen.cache2plot[sid.id]
        return (plot, sid.index)
    else
        return (nothing, 0)
    end
end

function AbstractPlotting.pick(screen::Screen, rect::IRect2D)
    window_size = widths(screen)
    buff = screen.framebuffer.objectid
    sid = zeros(SelectionID{UInt16}, widths(rect)...)
    glReadBuffer(GL_COLOR_ATTACHMENT1)
    x, y = minimum(rect)
    rw, rh = widths(rect)
    w, h = window_size
    if x > 0 && y > 0 && x <= w && y <= h
        glReadPixels(x, y, rw, rh, buff.format, buff.pixeltype, sid)
        return map(unique(vec(SelectionID{Int}.(sid)))) do sid
            (screen.cache2plot[sid.id], sid.index)
        end
    end
    return Tuple{AbstractPlot, Int}[]
end


text:
function compute()
	DIGITS=1000
	prev=BigInt(1)
	cur=BigInt(0)
	index=0;

	while true

		if length(string(cur))>DIGITS
			throw(RuntimeError("Not found"))
		elseif length(string(cur))==DIGITS
			return string(index)
		end

		prev,cur=cur,prev+cur
		index+=1
	end

end

println(compute())


text:
type DropConnectLayer <: Layer
    filters::Vector{Vol}
    biases::Vol
    inputActivation::Vol
    outputActivation::Vol

    keepProbability::Float64
    filterMasks::Vector{BitArray{3}}
    biasMask::BitArray{3}

    function DropConnectLayer(numOutputs::Int, bias::Float64, keepProbability::Float64, dimensions)
        filters = [Vol(dimensions) for _ in 1:numOutputs]
        masks = [BitArray(dimensions) for _ in 1:numOutputs]
        biases = Vol((1, 1, numOutputs), bias)
        new(filters, biases, Vol(dimensions, 0.0), Vol((1, 1, numOutputs), 0.0), keepProbability, masks)
    end
end

function forward(self::DropConnectLayer, inputActivation::Vol)
    @assert size(self.inputActivation) == size(inputActivation)
    self.inputActivation = inputActivation
    self.biasMask = rand(size(biases)) .< self.keepProbability

    for i in 1:length(self.filters)
        filter = self.filters[i]
        self.filterMasks[i] = rand(size(filter)) .< self.keepProbability
        @assert size(self.inputActivation) == size(filter)

        self.outputActivation.w[1, 1, i] =
            sum(filter.w .* inputActivation.w .* self.filterMasks[i]) +
            self.biases.w[1, 1, i] * self.biasMask[1, 1, i]
    end
    return self.outputActivation
end

function backward(self::DropConnectLayer)
    fill!(self.inputActivation.dw, 0.0)
    for i in 1:length(self.filters)
        filter = self.filters[i]
        chainGrad = self.outputActivation.dw[1, 1, i]
        self.inputActivation.dw += filter.w .* self.filterMasks[i] * chainGrad
        filter.dw += self.inputActivation.w .* self.filterMasks[i] * chainGrad
        self.biases.dw[1, 1, i] += chainGrad * self.biasMask[1, 1, i]
    end
end

function gradientStates(self::DropConnectLayer)
    grads = [GradientState(filter, 1.0, 1.0) for filter in self.filters]
    append!(grads, [GradientState(self.biases, 0.0, 0.0)])
    return grads
end


text:
#=

Authors: Shunsuke Hori

=#

using QuantEcon

"""
The approximate Coleman operator, updated using the endogenous grid
method.

Parameters
----------
g : Function
    The current guess of the policy function
k_grid : AbstractVector
    The set of *exogenous* grid points, for capital k = y - c
β : AbstractFloat
    The discount factor
u_prime : Function
    The derivative u'(c) of the utility function
u_prime_inv : Function
    The inverse of u' (which exists by assumption)
f : Function
    The production function f(k)
f_prime : Function
    The derivative f'(k)
shocks : AbstractVector
    An array of draws from the shock, for Monte Carlo integration (to
    compute expectations).
"""

function coleman_egm(g::Function,
                     k_grid::AbstractVector,
                     β::AbstractFloat,
                     u_prime::Function,
                     u_prime_inv::Function,
                     f::Function,
                     f_prime::Function,
                     shocks::AbstractVector)

    # Allocate memory for value of consumption on endogenous grid points
    c = similar(k_grid)

    # Solve for updated consumption value
    for (i, k) in enumerate(k_grid)
        vals = u_prime.(g.(f(k) * shocks)) .* f_prime(k) .* shocks
        c[i] = u_prime_inv(β * mean(vals))
    end

    # Determine endogenous grid
    y = k_grid + c  # y_i = k_i + c_i

    # Update policy function and return
    Kg = LinInterp(y,c)
    Kg_f(x)=Kg(x)
    return Kg_f
end


text:
module TestSolid

import Modia3D
using  Modia3D.Test
using  Modia3D.LinearAlgebra
EYE3() = Matrix(1.0I,3,3)

# Test Sphere
r = 2.0
m = 3.0
geo   = Modia3D.SolidSphere(2r)
mass  = Modia3D.MassProperties(geo,m)

rref  = [1.0,2.0,3.0]
Tref  = EYE3()
eref  = rref/norm(rref)
sp    = Modia3D.supportPoint(geo,rref,Tref,eref)
AABB1 = Modia3D.BoundingBox()
AABB2 = Modia3D.BoundingBox(rref[1]-r, rref[1]+r, rref[2]-r, rref[2]+r, rref[3]-r, rref[3]+r)
Modia3D.boundingBox!(geo,AABB1,rref,Tref)

@testset "Modia3D.Solid: Test Sphere" begin
   @test mass.rCM == [0.0, 0.0, 0.0]
   @test isapprox(Modia3D.volume(geo), 4/3*pi*r^3)
   @test isapprox(mass.I             , 2/5*m*r^2*EYE3())
   @test isapprox(sp                 , rref + r*eref)
   @test isapprox(AABB1.x_min        , AABB2.x_min)
   @test isapprox(AABB1.x_max        , AABB2.x_max)
   @test isapprox(AABB1.y_min        , AABB2.y_min)
   @test isapprox(AABB1.y_max        , AABB2.y_max)
   @test isapprox(AABB1.z_min        , AABB2.z_min)
   @test isapprox(AABB1.z_max        , AABB2.z_max)
end

# Test Ellipsoid
a = 1.0
b = 2.0
c = 3.0
m = 4.0
geo  = Modia3D.SolidEllipsoid(2a,2b,2c)
mass = Modia3D.MassProperties(geo,m)

rref  = [1.0,2.0,3.0]
Tref  = EYE3()
eref  = [0.0, 0.0, 1.0]
sp    = Modia3D.supportPoint(geo,rref,Tref,eref)
AABB1 = Modia3D.BoundingBox()
AABB2 = Modia3D.BoundingBox(rref[1]-a, rref[1]+a, rref[2]-b, rref[2]+b, rref[3]-c, rref[3]+c)
Modia3D.boundingBox!(geo,AABB1,rref,Tref)

@testset "Modia3D.Solid: Test Ellipsoid" begin
   @test mass.rCM == [0.0, 0.0, 0.0]
   @test isapprox(Modia3D.volume(geo), 4/3*pi*a*b*c)
   @test isapprox(mass.I             , m/5*Diagonal([b^2+c^2, c^2+a^2, a^2+b^2]))
   @test isapprox(sp                 , rref + c*eref)
   @test isapprox(AABB1.x_min        , AABB2.x_min)
   @test isapprox(AABB1.x_max        , AABB2.x_max)
   @test isapprox(AABB1.y_min        , AABB2.y_min)
   @test isapprox(AABB1.y_max        , AABB2.y_max)
   @test isapprox(AABB1.z_min        , AABB2.z_min)
   @test isapprox(AABB1.z_max        , AABB2.z_max)
end

# Test Box and MassProperties
a = 1.0
b = 2.0
c = 3.0
m = 4.0
d = 2700.0
solidMaterial = Modia3D.SolidMaterial(density=d)
geo   = Modia3D.SolidBox(a,b,c)
mass1 = Modia3D.MassProperties(geo,m)
mass2 = Modia3D.MassProperties(geo,solidMaterial)
mass3 = Modia3D.MassProperties(geo,"Aluminium")
V     = Modia3D.volume(geo)
m2    = mass2.m

rref  = [1.0,2.0,3.0]
Tref  = EYE3()
eref  = [0.0, 0.0, 1.0]
sp    = Modia3D.supportPoint(geo,rref,Tref,eref)
AABB1 = Modia3D.BoundingBox()
AABB2 = Modia3D.BoundingBox(rref[1]-a/2-geo.rsmall, rref[1]+a/2+geo.rsmall, rref[2]-b/2-geo.rsmall, rref[2]+b/2+geo.rsmall, rref[3]-c/2-geo.rsmall, rref[3]+c/2+geo.rsmall)
Modia3D.boundingBox!(geo,AABB1,rref,Tref)

@testset "Modia3D.Solid: Test Box and MassProperties" begin
   @test mass1.rCM == [0.0, 0.0, 0.0]
   @test isapprox(V          , a*b*c)
   @test isapprox(mass1.I    , m/12*Diagonal([b^2+c^2, c^2+a^2, a^2+b^2]))
   @test isapprox(mass2.m    , d*V)
   @test isapprox(mass2.I    , m2/12*Diagonal([b^2+c^2, c^2+a^2, a^2+b^2]))
   @test isapprox(mass3.m    , mass2.m)
   @test isapprox(mass3.I    , mass2.I)
  # @test isapprox(AABB1.x_min, AABB2.x_min)
  # @test isapprox(AABB1.x_max, AABB2.x_max)
  # @test isapprox(AABB1.y_min, AABB2.y_min)
  # @test isapprox(AABB1.y_max, AABB2.y_max)
  # @test isapprox(AABB1.z_min, AABB2.z_min)
  # @test isapprox(AABB1.z_max, AABB2.z_max)
end

# Test MassProperties
m   = 2.0
rCM = [1,2,3]
II  = [5 1 2;
       1 6 3;
       2 3 7]

mass1 = Modia3D.MassProperties(m,rCM,II)
mass2 = Modia3D.MassProperties(m=m,rCM=rCM,Ixx=5,Iyy=6,Izz=7,Ixy=1,Ixz=2,Iyz=3)

@testset "Modia3D.Solid: Test MassProperties" begin
   @test isapprox(mass1.m, m)
   @test isapprox(mass2.m, m)

   @test mass1.rCM == rCM
   @test mass2.rCM == rCM

   @test mass1.I == II
   @test mass2.I == II
end

end