text:
# This file holds the definition of the functions pertaining to the State type.

# The functions of the State type require the Attrition type.
requiredTypes = [ "attrition", "state" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export setName!,
       addRequirement!,
       removeRequirement!,
       clearRequirements!,
       setInitial!,
       setStateTarget!,
       setStateRetirementAge!,
       setStateAttritionScheme!

"""
```
setName!( state::State,
          name::String )
```
This function sets the name of the state `state` to `name`.

This function returns `nothing`.
"""
function setName!( state::State, name::String )::Void

    state.name = name
    return

end  # setName!( state, name )

"""
```
addRequirement!( state::State,
                 attribute::String,
                 value::String )
```
This function adds a requirement to the state `state`, requiring that the
attribute `attribute` has the value `value`. The spaces in the attribute are
replaced with underscores (`_`). If there is already a requirement on the
attribute, the function overwrites the requirement.

This function returns `nothing`.
"""
function addRequirement!( state::State, attribute::String, value::String )::Void

    if haskey( state.requirements, attribute )
        warn( "State already has a requirement on attribute '$attribute'." )
    end  # if haskey( state.requirements, tmpAttr )

    state.requirements[ attribute ] = [ value ]
    return

end  # addRequirement!( state, attribute, value )

"""
```
addRequirement!( state::State,
                 attribute::String,
                 values::Vector{String} )
```
This function adds a requirement to the state `state`, requiring that the
attribute `attribute` has a value in `values`. The spaces in the attribute are
replaced with underscores (`_`). If there is already a requirement on the
attribute, the function overwrites the requirement.

This function returns `nothing`.
"""
function addRequirement!( state::State, attribute::String,
    values::Vector{String} )::Void

    if haskey( state.requirements, attribute )
        warn( "State already has a requirement on attribute '$attribute'." )
    end  # if haskey( state.requirements, tmpAttr )

    state.requirements[ attribute ] = values
    return

end  # addRequirement!( state, attribute, values )

"""
```
removeRequirement!( state::State,
                    attribute::String )
```
This function removes the requirement for state `state` on attribute
`attribute`. If the attribute had no requirement on it, nothing happens.

This function returns `nothing`.
"""
function removeRequirement!( state::State, attribute::String )::Void

    delete!( state.requirements, attribute )
    return

end  # removeRequirement!( state, attribute )

"""
```
clearRequirements!( state::State )
```
This function clears all the requirements for state `state`.

This function returns `nothing`.
"""
function clearRequirements!( state::State )::Void

    empty!( state.requirements )
    return

end  # clearRequirements!( state )

"""
```
setInitial!( state::State,
             isInitial::Bool )
```
This function sets the flag of `state` which indicates the state is an initial
state to `isInitial`.

This function returns `nothing`.
"""
function setInitial!( state::State, isInitial::Bool )::Void

    state.isInitial = isInitial
    return

end  # setInitial!( state, isInitial )

"""
```
setStateTarget!( state::State,
                 target::T )
    where T <: Integer
```
This function sets the target number of personnel members in state `state` to
`target`. If the number is less than zero, it means there's no target.

This function returns `nothing`.
"""
function setStateTarget!( state::State, target::T )::Void where T <: Integer

    state.stateTarget = target < 0 ? -1 : target
    return

end  # setStateTarget!( state, target )

"""
```
setStateAttritionScheme!( state::State,
                          attrScheme::Attrition )
```
This function sets the attrition scheme of the state `state` to `attrScheme`.

This function returns `nothing.`
"""
function setStateAttritionScheme!( state::State, attrScheme::Attrition )::Void

    state.attrScheme = attrScheme
    return

end  # setStateAttritionScheme!( state, attrScheme )

"""
```
setStateAttritionScheme!( state::State,
                          attrRate::Float64,
                          attrPeriod::Float64,
                          mpSim::ManpowerSimulation )
```
This function sets the attrition scheme of the state `state` to a new attrition
scheme with an attrition rate of `attrRate` per period of `attrPeriod`, and
stores the scheme in the manpower simulation `mpSim`.

This function returns `nothing.`
"""
function setStateAttritionScheme!( state::State, attrRate::Float64,
    attrPeriod::Float64, mpSim::ManpowerSimulation )::Void

    newAttrScheme = Attrition( "Attrition:" * state.name, attrRate, attrPeriod )
    addAttritionScheme!( newAttrScheme, mpSim )
    state.attrScheme = newAttrScheme
    return

end  # setStateAttritionScheme!( state, attrScheme )

"""
```
setStateRetirementAge!( state::State,
                        retAge::T )
    where T <: Real
```
This function sets the retirement age for personnel in the state `state` to
`retAge`. If the age is set to 0, the default retirement scheme is used instead.

This function returns `nothing`.
"""
function setStateRetirementAge!( state::State, retAge::T )::Void where T <: Real

    if retAge < 0
        warn( "Retirement age must be > 0. Not making change." )
        return
    end  # if retAge < 0

    state.stateRetAge = retAge
    return

end  # setStateRetirementAge( state::State, retAge::T )

"""
```
setStateAttritionScheme!( state::State,
                          attrName::String,
                          mpSim::ManpowerSimulation )
```
This function sets the attrition scheme of the state `state` to the scheme with
name `attrName` as defined in the manpower simulation `mpSim`. If the name is
'default' or unknown, the default scheme is set.

This function returns `nothing.`
"""
function setStateAttritionScheme!( state::State, attrName::String,
    mpSim::ManpowerSimulation )::Void

    if ( lowercase( attrName ) == "default" ) ||
        !haskey( mpSim.attritionSchemes, attrName )
        state.attrScheme = mpSim.defaultAttritionScheme
        return
    end  # if ( lowercase( attrName ) == "default" ) || ...

    state.attrScheme = mpSim.attritionSchemes[ attrName ]
    return

end  # setStateAttritionScheme!( state, attrName, mpSim )

function Base.show( io::IO, state::State )

    print( io, "  State: $(state.name)" )
    print( io, "\n    Associated attrition scheme: $(state.attrScheme)" )

    if isempty( state.requirements )
        print( io, "\n    State '$(state.name)' has no requirements" )
        return
    else
        print( io, "\n    Requirements" )

        for attr in keys( state.requirements )
            print( io, "\n      $attr " )
            vals = state.requirements[ attr ]
            multival = length( vals ) != 1
            print( io, multival ? "∈ { " : "= " )
            print( io, join( map( val -> "'$val'", vals ), ", " ) )
            print( io, multival ? " }" : "" )
        end  # for attr in keys( state.requirements )
    end  # if isempty( state.requirements )

    if state.stateTarget >= 0
        print( io, "\n    Target of $(state.stateTarget) personnel members in state." )
    end  # if state.stateTarget >= 0

end  # show( io, state )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

"""
```
readState( sheet::XLSX.Worksheet,
           stateCat::XLSX.Worksheet,
           sLine::T )
```
This function reads line `sline` of the Excel sheet `sheet`, and uses the state
catalogue sheet `stateCat` to generate a single state from it.

This function returns a `Tuple{State, T}` where `T` is either a `String` or a
`Tuple{Float64, Float64}`. The first element is the state object as it is
described in the Excel sheet, the second element is either the name of the
attached attrition scheme, or the attrition period/rate pair describing
attrition in this state.
"""
function readState( sheet::XLSX.Worksheet, stateCat::XLSX.Worksheet,
    sLine::T ) where T <: Integer

    newState = State( string( sheet[ "A$sLine" ] ) )
    catLine = sheet[ "C$sLine" ]

    if isa( catLine, Missings.Missing )
        error( "State '$(newState.name)' not defined in catalogue." )
    end  # if isa( catLine, Void )

    stateTarget = sheet[ "B$sLine" ]
    setStateTarget!( newState, isa( stateTarget, Missings.Missing ) ? -1 :
        stateTarget )
    catLine += 1
    setInitial!( newState, stateCat[ "B$catLine" ] == "YES" )

    # Read attrition scheme.
    isFixedAttr = stateCat[ "C$catLine" ] == "YES"

    if isFixedAttr
        attrPar = Float64( stateCat[ "D$catLine" ] ),
            Float64( stateCat[ "E$catLine" ] )

        if any( par -> isa( par, Missings.Missing ), attrPar )
            error( "Attrition parameters for state '$(newState.name)' aren't properly defined." )
        end  # if any( par -> isa( par, Missings.Missing ), attrPar )

    else
        attrPar = stateCat[ "F$catLine" ]

        if isa( attrPar, Missings.Missing )
            attrPar = "default"
        end  # if isa( attrPar, Missings.Missing )

    end

    # Read state requirements/updates.
    nReqs = Int( stateCat[ "G$catLine" ] )

    for ii in 1:nReqs
        addRequirement!( newState,
            stateCat[ XLSX.CellRef( catLine, 6 + 2 * ii ) ],
            stateCat[ XLSX.CellRef( catLine, 7 + 2 * ii ) ] )
    end  # for ii in 1:nReqs

    return newState, attrPar

end  # readState( sheet, stateCat, sLine )

"""
```
processStateOptions( opts::String )
```
This function processes the list of options given in `opts`.

This function returns either a `String` if there was only one option, or a
`Vector{String}` if there were multiple, comma separated options.
"""
function processStateOptions( opts::String )

    vals = split( opts, "," )
    vals = map( val -> String( strip( val ) ), vals )
    return length( vals ) == 1 ? vals[ 1 ] : vals

end  # processStateOptions( opts )

"""
```
isPersonnelOfState( persAttrs::Dict{String, Any},
                    state::State )
```
This function tests if the personnel members with initialised attributes
`persAttrs` satisfies the requirements of state `state`.

This function returns a `Bool` with the result of the test.
"""
function isPersonnelOfState( persAttrs::Dict{String, Any}, state::State )::Bool

    for attr in keys( state.requirements )
        # If the attribute isn't initialised for the personnel member, the
        #   personnel doesn't satisfy it, and isn't in the state. Otherwise,
        #   the attribute's value must match with the requirements of the state.
        if !haskey( persAttrs, attr ) ||
            persAttrs[ attr ] ∉ state.requirements[ attr ]
            return false
        end  # if !haskey( persAttrs ) || ...
    end  # for attr in keys( state.requirements )

    return true

end  # isPersonnelOfState( persAttrs, state )

"""
```
orderStates( mpSim::ManpowerSimulation )
```
This function orders the states in the manpower simulation `mpSim` in such a
manner that a lower ordered state does not have any transitions to a higher
ordered state.

This function returns a `Vector{String}`, the list of ordered states.
"""
function orderStates( mpSim::ManpowerSimulation )::Vector{String}

    # Build a graph of the network.
    graph = ManpowerPlanning.buildTransitionNetwork( mpSim,
        keys( mpSim.stateList )... )[ 1 ]

    # Initialises the lists.
    sortedNodes = Vector{Int}()
    unsortedNodes = collect( vertices( graph ) )[ 1:(end-2) ]

    # For each node, find the first state in the sorted node list for which a
    #   transition state -> node exists. Insert the new node before the found
    #   transition.
    while !isempty( unsortedNodes )
        node = unsortedNodes[ 1 ]

        # Only insert the new node if it's the first one, or if it has any links
        #   to already sorted states. Otherwise, put it at the end of the
        #   unsorted node list.
        if isempty( sortedNodes ) || any( jj -> has_edge( graph, node, jj ) ||
                has_edge( graph, jj, node ), sortedNodes )
            firstFailIndex = findfirst( jj -> has_edge( graph, jj, node ),
                sortedNodes )
            lastFailIndex = findlast( jj -> has_edge( graph, node, jj ),
                sortedNodes )

            if firstFailIndex == 0
                # if lastFailIndex == 0
                    push!( sortedNodes, node )
                # else
                #     insert!( sortedNodes, lastFailIndex + 1, node )
                # end  # if lastFailIndex == 0
            else
                insert!( sortedNodes, firstFailIndex, node )
            end  # if firstFailIndex == 0

            deleteat!( unsortedNodes, 1 )
        else
            unsortedNodes = vcat( unsortedNodes[ 2:end ], unsortedNodes[ 1 ] )
        end  # if isempty( sortedNodes ) || ...

    end  # while !isempty( unsortedNodes )

    return map( ii -> get_prop( graph, ii, :state ), sortedNodes )

end  # orderStates( mpSim::ManpowerSimulation )


text:
using BinDeps
using Compat

@BinDeps.setup

windllname = "libipopt-1"
libipopt = library_dependency("libipopt", aliases=[windllname])

ipoptname = "Ipopt-3.12.1"

provides(Sources, URI("http://www.coin-or.org/download/source/Ipopt/$ipoptname.tgz"),
    libipopt, os = :Unix)

prefix=joinpath(BinDeps.depsdir(libipopt),"usr")
patchdir=BinDeps.depsdir(libipopt)
srcdir = joinpath(BinDeps.depsdir(libipopt),"src",ipoptname)

# fpu_control flag for building on ARM
if Sys.ARCH == :arm
    fpu_control = "ADD_CFLAGS=-DNO_fpu_control"
else
    fpu_control = ""
end

provides(SimpleBuild,
    (@build_steps begin
        GetSources(libipopt)
        @build_steps begin
            ChangeDirectory(srcdir)
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Blas"))
                CreateDirectory("build", true)
                `./get.Blas`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Blas","build"))
                `../configure --prefix=$prefix --disable-shared --with-pic`
                `make install`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Lapack"))
                CreateDirectory("build", true)
                `./get.Lapack`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Lapack","build"))
                `../configure --prefix=$prefix --disable-shared --with-pic
                              --with-blas="$prefix/lib/libcoinblas.a -lgfortran"`
                `make install`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","ASL"))
                `./get.ASL`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Mumps"))
                `./get.Mumps`
            end
            `./configure --prefix=$prefix coin_skip_warn_cxxflags=yes
                         --with-blas="$prefix/lib/libcoinblas.a -lgfortran"
                         --with-lapack=$prefix/lib/libcoinlapack.a
                         $fpu_control`
            `make`
            `make test`
            `make -j1 install`
        end
    end),libipopt, os = :Unix)

# OS X
if is_apple()
    using Homebrew
    provides(Homebrew.HB, "staticfloat/juliadeps/ipopt", libipopt, os = :Darwin)
end

# Windows
if is_windows()
    using WinRPM
    provides(WinRPM.RPM, "Ipopt", [libipopt], os = :Windows)
end

# Avoid Issue #62, building into /lib64 on OpenSUSE
configsite = nothing
if haskey(ENV, "CONFIG_SITE")
    configsite = pop!(ENV, "CONFIG_SITE")
end

@BinDeps.install Dict(:libipopt => :libipopt)

if configsite !== nothing
    ENV["CONFIG_SITE"] = configsite
end

ext: jl


text:
# Is left
# ---------------------------------------------------------------------------------------------
# If the point is left of the line segment in the 2D plane.
#
# The segment's direction is from 𝘅₁ to 𝘅₂. Let 𝘂 = 𝘅₂ - 𝘅₁ and 𝘃 = 𝗽 - 𝘅₁
# We may determine if the angle θ between the point and segment is in [0, π] based on the
# sign of 𝘂 × 𝘃, since 𝘂 × 𝘃 = ‖𝘂‖‖𝘃‖sin(θ).
#   𝗽    ^
#   ^   /
# 𝘃 |  / 𝘂
#   | /
#   o
# We allow points on the line (𝘂 × 𝘃 = 0) to be left, since this test is primarily
# used to determine if a point is inside a polygon. A mesh is supposed to partition
# its domain, so if we do not allow points on the line, there will exist points in the
# mesh which will not be in any face.
@inline function isleft(p::Point2D, l::LineSegment2D)
    return 0 ≤ l.𝘂 × (p - l.𝘅₁)
end

# Is left
# ---------------------------------------------------------------------------------------------
# If the point is left of the quadratic segment in the 2D plane.
#   𝗽    ^
#   ^   /
# 𝘃 |  / 𝘂
#   | /
#   o
# If the segment is straight, we can perform the isleft check with the straight
# line from the  segment's start point to the segment's end point.
# If this condition isn't met, the segment's curvature must be accounted for.
# We find the point on the curve q that is nearest point to the point of interest.
# Call this point q_near. We then perform the isleft check with the tangent vector
# of q at q_near and the vector from q_near to p, p - q_near.
function isleft(p::Point2D, q::QuadraticSegment2D)
    if isstraight(q)
        𝘃₁ = q.𝘅₂ - q.𝘅₁
        𝘃₂ = p - q.𝘅₁
        return 𝘃₁ × 𝘃₂ > 0
    else
        # See nearest_point for an explanation of the math.
        # When the cubic has 3 real roots, the point must be inside the
        # curve of the segment. Meaning:
        #   If the segment curves left, the point is right.
        #   If the segment curves right, the point is left.
        # This way we save substantial time by bypassing the complex number arithmetic
        𝘂 = q.𝘂
        𝘃 = q.𝘃
        𝘄 = p - q.𝘅₁
        # f′(r) = ar³ + br² + cr + d = 0
        a = 4(𝘂 ⋅ 𝘂)
        b = 6(𝘂 ⋅ 𝘃)
        c = 2((𝘃 ⋅ 𝘃) - 2(𝘂 ⋅𝘄))
        d = -2(𝘃 ⋅ 𝘄)
        # Lagrange's method
        e₁ = s₀ = -b/a
        e₂ = c/a
        e₃ = -d/a
        A = 2e₁^3 - 9e₁*e₂ + 27e₃
        B = e₁^2 - 3e₂
        if A^2 - 4B^3 > 0 # one real root
            s₁ = ∛((A + √(A^2 - 4B^3))/2)
            if s₁ == 0
                s₂ = s₁
            else
                s₂ = B/s₁
            end
            r = (s₀ + s₁ + s₂)/3
            𝘃₁ = 𝗗(q, r)
            𝘃₂ = p - q(r)
            return 𝘃₁ × 𝘃₂ > 0
        else # three real roots
            return (q.𝘅₂ - q.𝘅₁) × (q.𝘅₃ - q.𝘅₁) < 0
        end
    end
end

# Nearest point
# ---------------------------------------------------------------------------------------------
# Find the point on 𝗾(r) closest to the point of interest 𝘆.
#
# Note: r ∈ [0, 1] is not necessarily true for this function, since it finds the minimizer
# of the function 𝗾(r), ∀r ∈ ℝ
# Find r which minimizes ‖𝘆 - 𝗾(r)‖, where 𝗾(r) = r²𝘂 + r𝘃 + 𝘅₁.
# This r also minimizes ‖𝘆 - 𝗾(r)‖²
# It can be shown that this is equivalent to finding the minimum of the quartic function
# ‖𝘆 - 𝗾(r)‖² = f(r) = a₄r⁴ + a₃r³ + a₂r² + a₁r + a₀
# The minimum of f(r) occurs when f′(r) = ar³ + br² + cr + d = 0, where
# 𝘄 = 𝘆 - 𝘅₁, a = 4(𝘂 ⋅ 𝘂), b = 6(𝘂 ⋅ 𝘃), c = 2[(𝘃 ⋅ 𝘃) - 2(𝘂 ⋅𝘄)], d = -2(𝘃 ⋅ 𝘄)
# Lagrange's method (https://en.wikipedia.org/wiki/Cubic_equation#Lagrange's_method)
# is used to find the roots.
function nearest_point(p::Point{Dim,T}, q::QuadraticSegment) where {Dim,T}
    𝘂 = q.𝘂
    𝘃 = q.𝘃
    𝘄 = p - q.𝘅₁
    # f′(r) = ar³ + br² + cr + d = 0
    a = 4(𝘂 ⋅ 𝘂)
    b = 6(𝘂 ⋅ 𝘃)
    c = 2((𝘃 ⋅ 𝘃) - 2(𝘂 ⋅𝘄))
    d = -2(𝘃 ⋅ 𝘄)
    # Lagrange's method
    e₁ = s₀ = -b/a
    e₂ = c/a
    e₃ = -d/a
    A = 2e₁^3 - 9e₁*e₂ + 27e₃
    B = e₁^2 - 3e₂
    if A^2 - 4B^3 > 0 # one real root
        s₁ = ∛((A + √(A^2 - 4B^3))/2)
        if s₁ == 0
            s₂ = s₁
        else
            s₂ = B/s₁
        end
        r = (s₀ + s₁ + s₂)/3
        return r, q(r)
    else # three real roots
        # Complex cube root
        t₁ = exp(log((A + √(complex(A^2 - 4B^3)))/2)/3)
        if t₁ == 0
            t₂ = t₁
        else
            t₂ = B/t₁
        end
        ζ₁ = Complex{T}(-1/2, √3/2)
        ζ₂ = conj(ζ₁)
        dist_min = typemax(T)
        r_near = zero(T)
        p_near = nan(Point{Dim,T})
        # Use the real part of each root
        for rᵢ in (real((s₀ +    t₁ +    t₂)/3),
                   real((s₀ + ζ₂*t₁ + ζ₁*t₂)/3),
                   real((s₀ + ζ₁*t₁ + ζ₂*t₂)/3))
            pᵢ = q(rᵢ)
            dist = distance²(pᵢ, p)
            if dist < dist_min
                dist_min = dist
                r_near = rᵢ
                p_near = pᵢ
            end
        end
        return r_near, p_near
    end
end

# Hyperplane
# ---------------------------------------------------------------------------------------------
Base.in(p::Point, plane::Hyperplane) = p.coord ⋅ plane.𝗻 ≈ plane.d
in_halfspace(p::Point, plane::Hyperplane) = p.coord ⋅ plane.𝗻 - plane.d ≥ 0

# Given a point p and line l that lie in the plane. Check that the point is left of the line
function isleft(p::Point3D, l::LineSegment3D, plane::Hyperplane3D)
    # Since p and l ∈ plane, l.𝘂 × (p - l.𝘅₁) must either by in the exact same direction
    # as plane.𝗻 or the exact opposite direction. If the direction is the same, the point
    # is left of the line.
    return 0 ≤ (l.𝘂 × (p - l.𝘅₁)) ⋅ plane.𝗻
end

# AABox
# ---------------------------------------------------------------------------------------------
@inline Base.in(p::Point2D, aab::AABox2D) = aab.xmin ≤ p[1] ≤ aab.xmax &&
                                            aab.ymin ≤ p[2] ≤ aab.ymax
@inline Base.in(p::Point3D, aab::AABox3D) = aab.xmin ≤ p[1] ≤ aab.xmax &&
                                            aab.ymin ≤ p[2] ≤ aab.ymax &&
                                            aab.zmin ≤ p[3] ≤ aab.zmax
# Point inside polygon
# ---------------------------------------------------------------------------------------------
# Test if a point is in a polygon for 2D points/polygons
function Base.in(p::Point2D, poly::Polygon{N, 2}) where {N}
    # Test if the point is to the left of each edge.
    for i ∈ 1:N-1
        isleft(p, LineSegment2D(poly[i], poly[i + 1])) || return false
    end
    return isleft(p, LineSegment2D(poly[N], poly[1]))
end
# Not necessarily planar
#function Base.in(p::Point3D, poly::Polygon{N, 3}) where {N}
#    # Check if the point is even in the same plane as the polygon
#    plane = Hyperplane(poly[1], poly[2], poly[3])
#    p ∈ plane || return false
#    # Test that the point is to the left of each edge, oriented to the plane
#    for i = 1:N-1
#        isleft(p, LineSegment3D(poly[i], poly[i + 1]), plane) || return false
#    end
#    return isleft(p, LineSegment3D(poly[N], poly[1]), plane)
#end
#

# Point inside triangle
# ---------------------------------------------------------------------------------------------
function Base.in(p::Point3D, tri::Triangle3D)
    # P ∈ ABC iff the surface normals of CCW triangles PAB, PBC, & PCA are equal.
    𝗮 = tri[1] - p
    𝗯 = tri[2] - p
    𝗰 = tri[3] - p
    𝗻₁= 𝗮 × 𝗯
    𝗻₂= 𝗯 × 𝗰
    d₁₂ = 𝗻₁ ⋅ 𝗻₂
    # Test the normals point the same direction relative to each other
    # and that surface normals are equivalent using 𝗻̂ ⋅ 𝗻̂ = 1
    # d₁₂ > 0 is redundant if the point is in the triangle, but it is a very
    # fast check that the point is in the plane of the triangle.
    (d₁₂ > 0 && d₁₂ ≈ norm(𝗻₁)*norm(𝗻₂)) || return false
    # We need only check the direction of the norm of the last triangle to
    # prove that the point is in the triangle
    return 𝗻₂ ⋅(𝗰 × 𝗮) > 0
end

# This performs much better than the default routine, which is logically equivalent.
# Better simd this way? Chaining isleft doesn't have the same performance improvement for
# triangles.
function Base.in(p::Point2D, quad::Quadrilateral2D)
    return isleft(p, LineSegment2D(quad[1], quad[2])) &&
           isleft(p, LineSegment2D(quad[2], quad[3])) &&
           isleft(p, LineSegment2D(quad[3], quad[4])) &&
           isleft(p, LineSegment2D(quad[4], quad[1]))
end

# Test if a 2D point is in a 2D quadratic polygon
function Base.in(p::Point2D, poly::QuadraticPolygon{N, 2, T}) where {N, T}
    # Test if the point is to the left of each edge.
    bool = true
    M = N ÷ 2
    for i ∈ 1:M-1
        if !isleft(p, QuadraticSegment2D(poly[i], poly[i + 1], poly[i + M]))
            return false
        end
    end
    return isleft(p, QuadraticSegment2D(poly[M], poly[1], poly[N]))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# multicopter
abstract type Multicopter <: AbstractEnv end

"""
Common state structure of Multicopter
"""
function State(multicopter::Multicopter)
    return function (p=zeros(3), v=zeros(3), R=SMatrix{3, 3}(I), ω=zeros(3))
        ComponentArray(p=p, v=v, R=R, ω=ω)
    end
end

"""
Common input saturation (default: not applied)
"""
function saturate(multicopter::Multicopter, u)
    # u
    error("Actuator saturation not applied: `saturate`")
end

"""
Common input to force and moment transformation (default: not applied)
"""
function input_to_force_moment(multicopter::Multicopter, u)
    # ν = u
    error("Transformation of input to force and moment not defined: `input_to_force_moment`")
end

"""
# Variables
## State
p ∈ ℝ^3: (inertial) position
v ∈ ℝ^3: (inertial) velocity
R ∈ so(3): direction cosine matrix (DCM) that maps a vector read in Inertial (I)-coord.
to the same vector read in Body (B)-coord.
For example, v_B = R*v_I.
Or, it can be interpreted as "inverse rotation" from I-frame to B-frame.
For example, x̂_I = R'*[1, 0, 0] where x̂_I is the x-axis of B-frame read in I-coord.
ω ∈ ℝ^3: angular rate of body frame w.r.t. inertial frame (I to B)

## (Virtual) input
f ∈ ℝ: total thrust
M ∈ ℝ^3: moment

## BE CAREFUL
The definition of DCM, R, is the transpose of the DCM introduced in [1].
It is because many packages including `ReferenceFrameRotations.jl` follows the definition of
"DCM such that v_B = R*v_I".

# Reference
[1] T. Lee, M. Leok, and N. H. McClamroch, “Geometric Tracking Control of a Quadrotor UAV on SE(3),” in 49th IEEE Conference on Decision and Control (CDC), Atlanta, GA, Dec. 2010, pp. 5420–5425. doi: 10.1109/CDC.2010.5717652.
"""
function __Dynamics!(multicopter::Multicopter)
    @unpack m, g, J = multicopter
    J_inv = inv(J)
    e3 = [0, 0, 1]
    # skew(x): ℝ³ → ℝ⁹ such that x×y = skew(x)*y
    skew(x) = [    0 -x[3]  x[2];
                x[3]     0 -x[1];
               -x[2]  x[1]    0]
    @Loggable function dynamics!(dX, X, p, t; f, M)
        @unpack p, v, R, ω = X
        @onlylog state = X
        @nested_log :input f, M
        Ω = skew(ω)
        dX.p = v
        dX.v = -(1/m)*f*R'*e3 + g*e3
        dX.R = -Ω*R
        dX.ω = J_inv * (-Ω*J*ω + M)
    end
end

"""
A basic example of dynamics for multicopter considering rotor inputs `u`.
You can use the following closure or extend the abvoe __Dynamics! for more general
models, e.g., faulted multicopters.
"""
function _Dynamics!(multicopter::Multicopter)
    @Loggable function dynamics!(dx, x, p, t; u)
        @nested_onlylog :input u_cmd = u
        @nested_log :input u_saturated = saturate(multicopter, u)
        @nested_log :input ν = input_to_force_moment(multicopter, u_saturated)
        f, M = ν[1], ν[2:4]
        @nested_log __Dynamics!(multicopter)(dx, x, (), t; f=f, M=M)
    end
end

# multicopters
include("quadcopters.jl")
include("hexacopters.jl")


text:
type Schur_KKT_solver <: abstract_schur_solver
    # abstract_KKT_system_solver
    ls_solver::abstract_linear_system_solver # the linear system solver we wish to use, see the folder `linear_system_solvers`
    factor_it::Class_iterate # iterate where the factorization is computed
    delta_x_vec::Array{Float64,1} # the amount we perturb the Hessian
    delta_s_vec::Array{Float64,1} # MAYBE DELETE? I DON'T REALLY USE THIS ANYMORE??????
    rhs::System_rhs # the right hand side
    dir::Class_point # the direction that was computed
    kkt_err_norm::Class_kkt_error
    rhs_norm::Float64
    pars::Class_parameters
    schur_diag::Array{Float64,1} # diagonal values of the schur complement, if any of these elements is negative matrix inertia is incorrect.

    ready::Symbol # is the KKT solver correctly intialized

    # Schur_KKT_solver only
    Q::SparseMatrixCSC{Float64,Int64}
    #M::SparseMatrixCSC{Float64,Int64} # Is two linear systems really necessary????
    #K::SparseMatrixCSC{Float64,Int64}
    current_it::Class_iterate # current iterate
    reduct_factors::Class_reduction_factors # the amount that we want to reduce

    function Schur_KKT_solver()
      this = new()
      this.ready = :not_ready

      return this
    end
end

function kkt_associate_rhs!(kkt_solver::abstract_schur_solver, iter::Class_iterate, reduct_factors::Class_reduction_factors, timer::class_advanced_timer)
    start_advanced_timer(timer, "KKT/rhs");

    kkt_solver.rhs = System_rhs(iter, reduct_factors)
    kkt_solver.dir.mu = -(1.0 - reduct_factors.mu) * get_mu(iter)
    kkt_solver.dir.primal_scale = -(1.0 - reduct_factors.P) * iter.point.primal_scale

    kkt_solver.reduct_factors = reduct_factors
    kkt_solver.current_it = iter

    pause_advanced_timer(timer, "KKT/rhs");
end

function form_system!(kkt_solver::abstract_schur_solver, iter::Class_iterate, timer::class_advanced_timer)
    # form the matrix before we factorize it

    start_advanced_timer(timer, "SCHUR")
    start_advanced_timer(timer, "SCHUR/form_system");

    # TODO build specialized schur complement code.
    kkt_solver.Q = eval_J_T_J(iter, iter.point.y ./ iter.point.s) + get_lag_hess(iter);
    kkt_solver.schur_diag = diag(kkt_solver.Q)
    kkt_solver.factor_it = iter;
    kkt_solver.ready = :system_formed

    pause_advanced_timer(timer, "SCHUR/form_system");
    pause_advanced_timer(timer, "SCHUR");
end

function update_delta_vecs!(kkt_solver::abstract_schur_solver, delta_x_vec::Array{Float64,1}, delta_s_vec::Array{Float64,1}, timer::class_advanced_timer)
    start_advanced_timer(timer, "SCHUR")
    start_advanced_timer(timer, "SCHUR/delta_vecs")
    kkt_solver.delta_x_vec = delta_x_vec
    kkt_solver.delta_s_vec = delta_s_vec

    if sum(abs.(delta_s_vec)) > 0.0
        error("Not implemented")
        #kkt_solver.Q[i,i] =
        #kkt_solver.M + eval_J_T_J(kkt_solver.factor_it, delta_s_vec) + spdiagm(delta_x_vec)
    else
        for i = 1:size(kkt_solver.Q,1)
            kkt_solver.Q[i,i] = kkt_solver.schur_diag[i] + delta_x_vec[i]
        end
    end

    kkt_solver.ready = :delta_updated
    pause_advanced_timer(timer, "SCHUR/delta_vecs")
    pause_advanced_timer(timer, "SCHUR")
end

function factor_implementation!(kkt_solver::abstract_schur_solver, timer::class_advanced_timer)
    return ls_factor!(kkt_solver.ls_solver, kkt_solver.Q, dim(kkt_solver.factor_it), 0, timer)
end

function compute_direction_implementation!(kkt_solver::Schur_KKT_solver, timer::class_advanced_timer)
    # compute a search direction
    start_advanced_timer(timer, "SCHUR")

    factor_it = kkt_solver.factor_it
    #∇a_org = get_jac(factor_it);
    y_org = get_y(factor_it);
    s_org = get_s(factor_it);

    rhs = kkt_solver.rhs

    #r1 + ∇a_org' * (( r3 + (r2 .* y_org) ) ./ s_org)
    start_advanced_timer(timer, "SCHUR/rhs");
    symmetric_primal_rhs = rhs.primal_r + rhs.comp_r ./ y_org
    Σ_vec = ( y_org ./ s_org )
    y_ = (rhs.primal_r .* Σ_vec + rhs.comp_r ./ s_org )
    schur_rhs = rhs.dual_r + eval_jac_T_prod(factor_it, y_)
    pause_advanced_timer(timer, "SCHUR/rhs");

    dir = kkt_solver.dir
    dir.x = solver_schur_rhs(schur_rhs, kkt_solver, timer)

    # there are two ways to update s and y
    if true
      #J = get_jac(kkt_solver.current_it)
      dir.y = -(eval_jac_prod(factor_it,dir.x) - symmetric_primal_rhs) .* Σ_vec
      #dir.y = -(J * dir.x - symmetric_primal_rhs) .* Σ_vec
      #dir.s = ( rhs.comp_r - dir.y .* s_org ) ./ y_org
      dir.s = eval_jac_prod(factor_it,dir.x) - rhs.primal_r
    else
      dir.s = eval_jac_prod(factor_it,dir.x) - rhs.primal_r
      dir.y = ( rhs.comp_r - dir.s .* y_org ) ./ s_org # (mu_target - s_cur .* y_cur)
    end

    check_for_nan(dir)

    start_advanced_timer(timer, "SCHUR/kkt_err");
    update_kkt_error!(kkt_solver, Inf, timer)
    pause_advanced_timer(timer, "SCHUR/kkt_err");
    pause_advanced_timer(timer, "SCHUR")
end

function solver_schur_rhs(schur_rhs::Vector, kkt_solver::abstract_schur_solver, timer::class_advanced_timer)
  # Solve Ax = b where b is formed by the function compute_direction_implementation!().
  # The main purpose of this function is to do iterative refinement.

  fit = kkt_solver.factor_it
  #∇a_org = get_jac(factor_it);
  y_org = get_y(fit);
  s_org = get_s(fit);
  #(symmetric_primal_rhs .* Σ_vec);Σ
  dir = kkt_solver.dir;
  pars = kkt_solver.pars

  Σ_vec = ( y_org ./ s_org )
  sqrt_Σ_vec = sqrt.( y_org ./ s_org )

  # generalize!!!
  output_level = pars.output_level
  res_old = schur_rhs
  if output_level >= 4
    println("res", 0, " ", rd(norm(res_old,2)))
  end

  dir_x = zeros(length(dir.x));
  if pars.kkt.ItRefine_BigFloat
    dir_x = convert(Array{BigFloat,1}, dir_x)
  end

  for i = 1:pars.kkt.ItRefine_Num
      if output_level >= 4
        println("res", i, " ", rd(Float64(norm(res_old,2))))
      end

      start_advanced_timer(timer, "SCHUR/iterative_refinement");
      dir_x += ls_solve(kkt_solver.ls_solver, res_old, timer);

      start_advanced_timer(timer, "SCHUR/iterative_refinement/residual");
      jac_res = eval_jac_T_prod( fit , Σ_vec .* eval_jac_prod(fit, dir_x) )
      hess_res = hess_product(fit, dir_x) + kkt_solver.delta_x_vec .* dir_x
      res_old = schur_rhs - ( jac_res + hess_res )
      pause_advanced_timer(timer, "SCHUR/iterative_refinement/residual");

      #res_old = res
      pause_advanced_timer(timer, "SCHUR/iterative_refinement");
  end

  if output_level >= 4
    println("res ", rd(Float64(norm(res_old,2))))
  end

  return dir_x
end


text:
export EuclidDistance, PackedEuclidDistance

"""
$(TYPEDEF)

Default linear offset between two scalar variables.
"""
struct EuclidDistance{T <: SamplableBelief} <: AbstractRelativeMinimize
  Z::T
end

EuclidDistance(::UniformScaling=LinearAlgebra.I) = EuclidDistance(Normal())

# consider a different group?
getManifold(::InstanceType{EuclidDistance}) = TranslationGroup(1)
getDimension(::InstanceType{<:EuclidDistance}) = 1

# new and simplified interface for both nonparametric and parametric
(s::CalcFactor{<:EuclidDistance})(z, x1, x2) = z .- norm(x2 .- x1)

Base.convert(::Type{<:MB.AbstractManifold}, ::InstanceType{EuclidDistance}) = Manifolds.TranslationGroup(1)

"""
$(TYPEDEF)
Serialization type for `EuclidDistance` binary factor.
"""
Base.@kwdef mutable struct PackedEuclidDistance <: AbstractPackedFactor
  _type::String
  Z::PackedSamplableBelief
end

function convert(::Type{PackedEuclidDistance}, d::EuclidDistance)
  PackedEuclidDistance( "/application/JuliaLang/PackedSamplableBelief",
                        convert(PackedSamplableBelief, d.Z) )
end

function convert(::Type{<:EuclidDistance}, d::PackedEuclidDistance)
  EuclidDistance(convert(SamplableBelief, d.Z))
end

#


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using Calculus

import Base.parse

function get_reference_element_coordinates end
function eval_basis! end
function eval_dbasis! end

function calculate_interpolation_polynomials(p, V)
    basis = []
    first(p.args) == :+ || error("Use only summation between terms of polynomial")
    args = p.args[2:end]
    N = size(V, 1)
    b = zeros(N)
    for i in 1:N
        fill!(b, 0.0)
        b[i] = 1.0
        solution = V \ b
        N = Expr(:call, :+)
        for (ai, bi) in zip(solution, args)
            isapprox(ai, 0.0) && continue
            push!(N.args, simplify( :( $ai * $bi ) ))
        end
        push!(basis, N)
    end
    return basis
end

function calculate_interpolation_polynomial_derivatives(basis, D)
    vars = [:u, :v, :w]
    dbasis = []
    for N in basis
        partial_derivatives = []
        for j in 1:D
            push!(partial_derivatives, simplify(differentiate(N, vars[j])))
        end
        push!(dbasis, partial_derivatives)
    end
    return dbasis
end

function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, p::Expr) where {N, D, T}
    @debug "create basis given antsatz polynomial" name description X p
    V = vandermonde_matrix(p, X)
    basis = calculate_interpolation_polynomials(p, V)
    return create_basis(name, description, X, basis)
end

function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis) where {N, D, T}
    @debug "create basis given basis functions" name description X basis
    dbasis = calculate_interpolation_polynomial_derivatives(basis, D)
    return create_basis(name, description, X, basis, dbasis)
end

function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis, dbasis) where {N, D, T}

    @debug "create basis given basis functions and derivatives" name description X basis dbasis

    Q = Expr(:block)
    for i=1:N
        push!(Q.args, :(N[$i] = $(basis[i])))
    end

    V = Expr(:block)
    for i=1:N
        for k=1:D
            push!(V.args, :(dN[$k,$i] = $(dbasis[i][k])))
        end
    end

    if D == 1
        unpack = :((u,) = xi)
    elseif D == 2
        unpack = :((u, v) = xi)
    else
        unpack = :((u, v, w) = xi)
    end

    code = quote
        struct $name <: FEMBasis.AbstractBasis
        end

        function Base.size(::Type{$name})
            return ($D, $N)
        end

        function Base.size(::Type{$name}, j::Int)
            j == 1 && return $D
            j == 2 && return $N
        end

        function Base.length(::Type{$name})
            return $N
        end

        function FEMBasis.get_reference_element_coordinates(::Type{$name})
            return $X
        end

        function FEMBasis.eval_basis!(::Type{$name}, N::Matrix{T}, xi) where T
            $unpack
            $Q
            return N
        end

        function FEMBasis.eval_dbasis!(::Type{$name}, dN::Matrix{T}, xi) where T
            $unpack
            $V
            return dN
        end
    end

    return code
end


text:
# TODO: FD
# !!! warning "Experimental feature"
#     This is an experimental feature and may change in any future releases.

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advection_velocity = (0.2, -0.7)
equations = LinearScalarAdvectionEquation2D(advection_velocity)

initial_condition = initial_condition_convergence_test

D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
                            derivative_order=1, accuracy_order=4,
                            xmin=0.0, xmax=1.0, N=100)
solver = DG(D_SBP, nothing #= mortar =#,
            SurfaceIntegralStrongForm(flux_lax_friedrichs),
            VolumeIntegralStrongForm())

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=1,
                n_cells_max=30_000,
                periodicity=true)

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

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

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan);

summary_callback = SummaryCallback()

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

alive_callback = AliveCallback(analysis_interval=analysis_interval)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback)

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

sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-9, reltol=1.0e-9,
            save_everystep=false, callback=callbacks, maxiters=1e5)
summary_callback()


text:
layers(::Type{EarthEnv{LandCover}}) = 1:12

"""
    getraster(T::Type{EarthEnv{LandCover}}, [layer::Union{AbstractArray,Tuple,Integer}]; discover::Bool=false) => String
    getraster(T::Type{EarthEnv{LandCover}}, layer::Integer, discover::Bool) => String

Download EarthEnv landcover data, choosing `layer` from: `$(layers(EarthEnv{LandCover}))`.

Without a layer argument, all layers will be downloaded and a tuple of paths returned.
If the data is already downloaded the path will be returned.
"""
function getraster(T::Type{EarthEnv{LandCover}}, layer::Integer; discover::Bool=false)
    getraster(T, layer, discover)
end
function getraster(T::Type{EarthEnv{LandCover}}, layer::Integer, discover::Bool)
    _check_layer(T, layer)
    url = rasterurl(T, layer; discover)
    path = rasterpath(T, layer; discover)
    return _maybe_download(url, path)
end

rasterpath(T::Type{EarthEnv{LandCover}}, layer::Integer; discover::Bool=false) =
    joinpath(rasterpath(T), rastername(T, layer; discover))
function rastername(::Type{EarthEnv{LandCover}}, layer::Integer; discover::Bool=false)
    filetype = discover ? "complete" : "partial"
    "landcover_$(filetype)_$(layer).tif"
end
function rasterurl(T::Type{EarthEnv{LandCover}}, layer::Integer; discover::Bool=false)
    stem = discover ? "with_DISCover/consensus_full_class_$(layer).tif" :
                      "without_DISCover/Consensus_reduced_class_$(layer).tif"
    joinpath(rasterurl(T), stem)
end

_pathsegment(::Type{LandCover}) = "consensus_landcover"


text:
# Test variational information

using Base.Test
using Clustering

a1 = [1, 1, 1, 2, 2, 2, 3, 3, 3, 3]
a2 = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]

@test_approx_eq_eps varinfo(3, a1, 3, a1) 0.0 1.0e-12
@test_approx_eq_eps varinfo(2, a2, 2, a2) 0.0 1.0e-12

v = varinfo(3, a1, 2, a2)
v_ = varinfo(2, a2, 3, a1)
@test 0.0 < v < log(3)
@test_approx_eq v v_

a1 = [1, 2, 3, 4, 5]
a2 = [1, 1, 1, 1, 1]

@test_approx_eq varinfo(5, a1, 1, a2) log(5)
@test_approx_eq varinfo(1, a2, 5, a1) log(5)

a1 = [1, 1, 1, 2, 2, 2]
a2 = [2, 2, 2, 1, 1, 1]
@test_approx_eq_eps varinfo(2, a1, 2, a2) 0.0 1.0e-12
@test_approx_eq_eps varinfo(2, a1, 3, a2) 0.0 1.0e-12
@test_approx_eq_eps varinfo(4, a1, 3, a2) 0.0 1.0e-12


text:
export bounce!, ispinned, propagate!

@inline increment_counter(::Int, t_to_write) = 1
@inline increment_counter(::T, t_to_write) where {T<:AbstractFloat} = t_to_write

#####################################################################################
# Bounce
#####################################################################################
"""
    bounce!(p::AbstractParticle, bd::Billiard) → i, t, pos, vel
"Bounce" the particle (advance for one collision) in the billiard.
Takes care of finite-precision issues.

Return:
* index of the obstacle that the particle just collided with
* the time from the previous collision until the current collision `t`
* position and velocity of the particle at the current collision (*after* the
  collision has been resolved!). The position is given in the unit cell of
  periodic billiards. Do `pos += p.current_cell` for the position in real space.

```julia
bounce!(p, bd, raysplit) → i, t, pos, vel
```
Ray-splitting version of `bounce!`.
"""
@inline function bounce!(p::AbstractParticle{T}, bd::Billiard{T}) where {T}
    i::Int, tmin::T, cp::SV{T} = next_collision(p, bd)
    if tmin != T(Inf)
        o = bd[i]
        relocate!(p, o, tmin, cp)
        resolvecollision!(p, o)
    end
    typeof(p) <: MagneticParticle && (p.center = find_cyclotron(p))
    return i, tmin, p.pos, p.vel
end

#####################################################################################
# Relocate
#####################################################################################
"""
    relocate!(p::AbstractParticle, o::Obstacle, t, cp)
Propagate the particle to `cp` and propagate velocities for time `t`.
Check if it is on the correct side of the obstacle. If not,
change the particle position by [`distance`](@ref) along the [`normalvec`](@ref)
of the obstacle.
"""
function relocate!(p::AbstractParticle{T}, o::Obstacle{T}, tmin::T, cp::SV{T}) where {T}
    propagate!(p, cp, tmin) # propagate to collision point
    okay = _correct_pos!(p, o)
    return okay
end
function _correct_pos!(p, o)
    d = distance(p.pos, o)
    okay = _okay(d, o)
    if !okay
        n = normalvec(o, p.pos)
        p.pos -= d * n
        # due to finite precision this does not always give positive distance
        # but collision takes care of it, as it does not care about collisions
        # very close to the point.
    end
    return okay
end

_okay(d, o::Obstacle) = d ≥ 0.0
_okay(d, o::PeriodicWall) = d ≤ 0.0

"""
    propagate!(p::AbstractParticle, t)
Propagate the particle `p` for given time `t`, changing appropriately the the
`p.pos` and `p.vel` fields.

    propagate!(p, position, t)
Do the same, but take advantage of the already calculated `position` that the
particle should end up at.
"""
@inline function propagate!(p::Particle{T}, t::Real) where {T}
    p.pos += SV{T}(p.vel[1]*t, p.vel[2]*t)
end
@inline propagate!(p::Particle, newpos::SV, t::Real) = (p.pos = newpos)

@inline function propagate!(p::MagneticParticle{T}, t::Real) where {T}
    ω = p.omega; r = 1/ω
    φ0 = atan(p.vel[2], p.vel[1])
    sinωtφ, cosωtφ = sincos(ω*t + φ0)
    p.pos += SV{T}((sinωtφ - sin(φ0))*r, (-cosωtφ + cos(φ0))*r)
    p.vel = SV{T}(cosωtφ, sinωtφ)
    return
end
@inline function propagate!(
        p::MagneticParticle{T}, newpos::SV{T}, t) where {T}
    ω = p.omega; φ0 = atan(p.vel[2], p.vel[1])
    p.pos = newpos
    p.vel = SV{T}(cossin(ω*t + φ0))
    return
end

#####################################################################################
# Resolve Collisions
#####################################################################################
"""
    specular!(p::AbstractParticle, o::Obstacle)
Perform specular reflection based on the normal vector of the Obstacle.

In the case where the given obstacle is a `RandomObstacle`, the specular reflection
randomizes the velocity instead (within -π/2+ε to π/2-ε of the normal vector).
"""
@inline function specular!(p::AbstractParticle{T}, o::Obstacle{T})::Nothing where {T}
    n = normalvec(o, p.pos)
    p.vel = p.vel - 2*dot(n, p.vel)*n
    return nothing
end

@inline function specular!(p::AbstractParticle{T}, r::RandomDisk{T})::Nothing where {T}
    n = normalvec(r, p.pos)
    φ = atan(n[2], n[1]) + 0.95(π*rand() - π/2) #this cannot be exactly π/2
    p.vel = SVector{2,T}(cos(φ), sin(φ))
    return nothing
end

@inline function specular!(p::AbstractParticle{T}, r::RandomWall{T})::Nothing where {T}
    n = normalvec(r, p.pos)
    φ = atan(n[2], n[1]) + 0.95(π*rand() - π/2) #this cannot be exactly π/2
    p.vel = SVector{2,T}(cossin(φ))
    return nothing
end

"""
    periodicity!(p::AbstractParticle, w::PeriodicWall)
Perform periodicity conditions of `w` on `p`.
"""
@inline function periodicity!(p::AbstractParticle, w::PeriodicWall)::Nothing
    p.pos += w.normal
    p.current_cell -= w.normal
    return nothing
end
@inline function periodicity!(p::MagneticParticle, w::PeriodicWall)::Nothing
    p.pos += w.normal
    p.center += w.normal
    p.current_cell -= w.normal
    return nothing
end

"""
    resolvecollision!(p::AbstractParticle, o::Obstacle)
Resolve the collision between particle `p` and obstacle `o`, depending on the
type of `o` (do `specular!` or `periodicity!`).

    resolvecollision!(p, o, T::Function, θ::Function, new_ω::Function)
This is the ray-splitting implementation. The three functions given are drawn from
the ray-splitting dictionary that is passed directly to `evolve!()`. For a calculated
incidence angle φ, if T(φ) > rand(), ray-splitting occurs.
"""
@inline resolvecollision!(p::Particle, o::Obstacle) = specular!(p, o)
@inline resolvecollision!(p::Particle, o::PeriodicWall) = periodicity!(p, o)
@inline resolvecollision!(p::MagneticParticle, o::Obstacle) = specular!(p, o)

resolvecollision!(p::MagneticParticle, o::PeriodicWall) = periodicity!(p, o)

"""
    ispinned(p::MagneticParticle, bd::Billiard)
Return `true` if the particle is pinned with respect to the billiard.
Pinned particles either have no valid collisions (go in circles forever)
or all their valid collisions are with periodic walls, which again means that
they go in cirles for ever.
"""
ispinned(p::Particle, bd) = false
function _reset_ispinned(p, pos, vel, cc)
    p.pos = pos; p.vel = vel
    p.current_cell = cc; p.center = find_cyclotron(p)
end
function ispinned(p::MagneticParticle, bd::Billiard)
    pos, vel, cc = p.pos, p.vel, p.current_cell
    i, t = bounce!(p, bd)
    if t == Inf; _reset_ispinned(p, pos, vel, cc); return true; end
    if !isperiodic(bd); _reset_ispinned(p, pos, vel, cc); return false; end

    peridx = bd.peridx
    if i ∉ peridx; _reset_ispinned(p, pos, vel, cc); return false; end

    # propagate until 2π/ω
    counter = t; limit = 2π/abs(p.omega)
    while counter ≤ limit
        i, t = bounce!(p, bd)
        if i ∉ peridx; _reset_ispinned(p, pos, vel, cc); return false; end
        counter += t
    end
    _reset_ispinned(p, pos, vel, cc); return true
end


text:
# Program for plotting Basin output, getting new initial conditions, etc...
include("funcs.jl")
using DelimitedFiles
using InteractiveViz
using IterTools
let

    ### getting directory for simulation data
    done = false
    cd(string("../../erickson/output_files"))
    while !done
        v = false
        while v == false
            @printf "\nOutput directory: "
            folder_name = chomp(readline())
            try
                cd(folder_name)
                v = true
            catch
                @printf "Directory not found\n"
            end
        end
        @printf "Found directory...\n\n"

        done1 = false
        while !done1

            ### options for plotting etc..
            @printf "\t(1) Slip Plot\n"
            @printf "\t(2) Station Data\n"
            @printf "\t(3) Single Cycle\n"
            @printf "\t(4) Back\n"
            @printf "\t(5) Done\n\n"

            v = false
            flag = [0]
            while v == false
                @printf "Option: "
                try
                    flag[1] = parse(Int64, chomp(readline()))
                    if 0 < flag[1] < 7
                        v = true
                    end
                catch
                    @printf "Not an Option\n"
                end

            end

            ### basic slip plot
            if flag[1] == 1

                slip_file = open("slip.dat", "r")
                slip_data = collect(eachline(slip_file))
                y, nn, slip_data = get_y_nn(slip_data)
                num_breaks, break_indices = get_break_indices(slip_data)

                total_cycles = floor(length(break_indices)/2)
                @printf "\n%d cycles in this data\n" total_cycles
                @printf "Offset contours by # of cycles: "
                cycle_offset = parse(Int64, chomp(readline()))
                @printf "Last cycle: "
                final_cycle = parse(Int64, chomp(readline()))

                final_index, index_offset = get_plot_indices(cycle_offset,
                                                             final_cycle,
                                                             break_indices)

                slip_data = @view slip_data[index_offset:final_index, :]
                plt = plot_slip(slip_data, y, nn, "")
                display(plt)
                close(slip_file)

            ### station plots
            elseif flag[1] == 2

                @printf "Station depth: "
                depth = chomp(readline())
                station_data = @view readdlm(string("station_", depth))[2:end, :]
                t = @view station_data[ : , 1]
                V = @view station_data[ : , 3]
                τ = @view station_data[ : , 4]

                #plot(t, V)
                plot(t, τ)

            ### single cycle options
            elseif flag[1] == 3

                # parse slip file to be able to find individual cycle
                slip_file = open("slip.dat", "r")
                slip_data = collect(eachline(slip_file))
                y_depth, nn, slip_data = get_y_nn(slip_data)
                break_number, break_indices = get_break_indices(slip_data)

                begin_index = 0
                δ = zeros(nn)
                t1 = 0.0
                cycle_number = 0
                done2 = false
                while !done2
                    @printf "\n%d cycles in this data\n" floor(length(break_indices)/2)
                    @printf "Which cycle: "

                    cycle_number = parse(Int64, chomp(readline()))

                    final_index, begin_index = get_plot_indices(cycle_number,
                                                                cycle_number,
                                                                break_indices)

                    # get the single cycle slip data
                    slip_cycle_data = @view slip_data[begin_index : final_index, :]
                    # get a wider view to show location of cycle
                    temp_plot = @view slip_data[break_indices[1] : final_index, :]
                    # set offset for view
                    δ_off = get_line(temp_plot, 1)[3:end]
                    # plot the wider view
                    plt1 = plot_slip(temp_plot, y_depth, nn, "")
                    display(plt1)
                    # get slip initial slip of this cycle
                    line = get_line(slip_cycle_data, 1)
                    δ .= line[3:end]
                    t1 = line[1]
                    # plot the intial slip in green to show start of cycle
                    plt1 = plot!(δ-δ_off, y_depth, linecolor=:green, linewidth=3)

                    # check if this is the right cycle
                    @printf "Is this the right cycle (y/n):"
                    anw = chomp(readline())
                    if anw == "y"
                        done2 = true
                    end
                end

                    done3 = false
                    while !done3

                        ### get single cycle options
                        @printf "\n\t(1) Write-out cycle initial conditions\n"
                        @printf "\t(2) Plot volume variables\n"
                        @printf "\t(3) Write vtk files\n"
                        @printf "\t(4) Fault variable animations\n"
                        @printf "\t(5) Back\n"

                        @printf "\nOption: "
                        option = parse(Int64, chomp(readline()))

                        ### new initial conditions
                        if option == 1

                            # get state initial condition
                            state_file = open("state.dat", "r")
                            state_data = collect(eachline(state_file))[2:end, :]
                            ### old indexing TO BE GOTTEN RID OF!
                            #state_index = begin_index - 2*(cycle_number - 1)
                            ###
                            line = get_line(state_data, begin_index)
                            ψ = line[3:end]
                            t2 = line[1]

                            plt2 = plot(δ, y_depth, yflip = true, legend = false)
                            plt3 = plot(ψ, y_depth, yflip = true, legend = false)
                            display(plot(plt2,plt3, layout=(2,1)))

                            # check that ψ and δ are coming from the same time
                            @assert t1 == t2

                            close(state_file)
                            close(slip_file)

                            @printf "\nTime is: %f \n" t1/year_seconds
                            @printf "\nAre these initial conditions ok (y/n): "
                            anw = chomp(readline())

                            ### if initial conditions look good write them out.
                            if anw == "y"
                                write_ψδ(ψ, δ, t1)
                            else
                                done3 = false
                            end

                            ### make the 4 pannel volume plot
                        elseif option == 2

                            @printf "Getting volume data...\n"
                            flush(stdout)
                            u_file, x, y, u_iter = parse_volume_file("us.dat", cycle_number)
                            v_file, _, _, v_iter = parse_volume_file("vs.dat", cycle_number)

                            @printf "Getting stress data...\n"
                            flush(stdout)
                            # get inter and coseismic stress data
                            stress_file = open("stress.dat", "r")

                            # this indexing only works for new "break" writting
                            stress_data = @view collect(eachline(stress_file))[2:end, :]
                            inter_stress_data = @view stress_data[begin_index : break_indices[2*cycle_number] - 2, :]
                            co_stress_data = @view stress_data[break_indices[2*cycle_number] : final_index, :]
                            @printf "Getting slip data...\n"
                            flush(stdout)
                            # get inter and coeseismic slip data
                            inter_slip_data = @view slip_data[begin_index : break_indices[2*cycle_number] - 2, :]
                            co_slip_data = @view slip_data[break_indices[2*cycle_number] : final_index, :]
                            δ_off = get_line(inter_slip_data, 1)[3:end]

                            # again checking to make sure data is aligned
                            t1 = get_line(co_stress_data, size(co_stress_data)[1])[1]
                            t2 = get_line(co_slip_data, size(co_slip_data)[1])[1]
                            @assert t1 == t2

                            slip_plot = plot(legend=false, yflip = true, ylabel="Depth(Km)", xlabel="Slip(m)")
                            stress_plot = plot(legend=false, yflip = true, ylabel="Depth(Km)", xlabel="Stress(MPa)")
                            u_plot = heatmap(c=cgrad([:blue, :white,:red]), yflip=true)
                            v_plot = heatmap(clim = (0.0, .75), c=cgrad([:blue, :white,:red]), yflip=true)

                            @printf "Beginning plotting...\n"
                            for i in 1:size(inter_stress_data)[1]

                                @printf "\r %f" i/size(inter_stress_data)[1]
                                δ = get_line(inter_slip_data, i)[3:end]
                                τ = get_line(inter_stress_data, i)[3:end]

                                slip_plot = plot!(slip_plot,
                                                  δ .- δ_off,
                                                  y_depth,
                                                  linecolor=:blue,
                                                  linewidth=.1)

                                stress_plot = plot(τ,
                                                   y_depth,
                                                   legend=false,
                                                   linecolor=:black,
                                                   flip = true,
                                                   ylabel="Depth(Km)",
                                                   xlabel="Stress(MPa)")

                                plot(slip_plot, v_plot, u_plot, stress_plot, layout=(2,2))
                                gui()

                            end

                            # loop over iterators
                            for (i, (u_line, v_line)) in enumerate(zip(u_iter, v_iter))
                                @printf "\r %f" i/size(co_stress_data)[1]
                                u = reshape(parse_line(u_line), length(x), length(y))'
                                v = reshape(parse_line(v_line), length(x), length(y))'

                                δ = get_line(co_slip_data, i)[3:end]
                                τ = get_line(co_stress_data, i)[3:end]

                                slip_plot = plot!(slip_plot,
                                                  δ .- δ_off,
                                                  y_depth,
                                                  linecolor=:red,
                                                  linewidth=.1)

                                stress_plot = plot(τ,
                                                   y_depth,
                                                   legend=false,
                                                   linecolor=:black,
                                                   flip = true,
                                                   ylabel="Depth(Km)",
                                                   xlabel="Stress(MPa)")

                                v_plot = heatmap(x, y, v, clim = (0.0, .75),
                                                 c=cgrad([:blue, :white,:red]), yflip=true)
                                u_plot = heatmap(x, y, u,
                                                 c=cgrad([:blue, :white,:red]), yflip=true)

                                plot(slip_plot, v_plot, u_plot, stress_plot, layout=(2,2))
                                gui()

                            end

                            close(u_file)
                            close(v_file)
                            close(slip_file)
                            close(stress_file)
                        elseif option == 3

                            # do animations for velocity, stress, and state
                        elseif option == 4
                            ### old indexing TO BE GOTTEN RID OF!
                            begin_index = begin_index - 2*(cycle_number - 1)
                            final_index = final_index - 2*(cycle_number - 2)
                            ###
                            files, plot_data = open_files("stress.dat",
                                                          "state.dat",
                                                          "slip_rate.dat",
                                                          begin_index = begin_index,
                                                          final_index = final_index)
                            fault_animations(plot_data, y_depth)

                        elseif option == 5
                            done2 = true
                        end
                    end

                    ### go back to look at different data
                elseif flag[1] == 4
                    done1 = true
                    cd("..")

                elseif flag[1] == 5
                    exit()
                end
            end
        end
        cd("/home/tharvey2/Basin_Simulations/utils")
    end


text:
import Base: >>, <<, &, and_int, lshr_int, shl_int, or_int, |, bswap,
            zero
import SortingAlgorithms: uint_mapping

primitive type  Bits24 24 end
primitive type  Bits192 192 end
primitive type  Bits256 256 end

>>(x::Bits24, y) = Base.lshr_int(x, y)
<<(x::Bits24, y) = Base.shl_int(x, y)

>>(x::Bits192, y) = Base.lshr_int(x, y)
<<(x::Bits192, y) = Base.shl_int(x, y)

(&)(x::Bits192, y::Bits192) = Base.and_int(x,y)
(|)(x::Bits192, y::Bits192) = Base.or_int(x,y)

function Bits192(x)
    # it loads from the end
    z = unsafe_load(Ptr{Bits192}(pointer_from_objref(x)))

    if sizeof(Bits192) > sizeof(x)
        lzx = leading_zeros(x)
        shift_n = sizeof(Bits192)*8 - sizeof(eltype(x))*8
        z = z << shift_n >> shift_n
    end
    z
end

const u192_mask = Bits192(2^16-1)

mask16bit(::Type{Bits192}) = u192_mask
mask16bit(::Type) = 0xffff

function make_mask(::Type{Bits192})
    x = UInt(2)^(sizeof(UInt)*8) - 1

    yy = unsafe_load(Ptr{Bits192}(pointer_from_objref(x)))
    for i = 2:3
        y = unsafe_load(Ptr{Bits192}(pointer_from_objref(x)))
        # get rid of zeros in front
        y = y << 8*(sizeof(Bits192) - sizeof(UInt)) >> (sizeof(Bits192) - sizeof(UInt))*8
        yy = (yy << 8*sizeof(UInt)) | y
    end
    yy
end

# UInt16(x::Bits192) = unsafe_load(Ptr{UInt16}(pointer_from_objref(x)))

bswap(x::Bits192) = Base.bswap_int(x)

uint_mapping(::Base.Order.ForwardOrdering, x::Bits192) = x

zero(::Type{Bits192}) = Bits192(0)

Int(x::Bits192) = Base.Intrinsics.trunc_int(Int, x)

convert(::Type{UInt16}, x::Bits192) = Base.Intrinsics.trunc_int(UInt16, x)

# promote_rule(::Type{Bits192}, ::Type{Integer}) = Bits192
(&)(x::Bits192, y::UInt16) = UInt16(x) & y

#########################################################################
#Bits 256
#########################################################################
>>(x::Bits256, y) = Base.lshr_int(x, y)
<<(x::Bits256, y) = Base.shl_int(x, y)

(&)(x::Bits256, y::Bits256) = Base.and_int(x,y)
(|)(x::Bits256, y::Bits256) = Base.or_int(x,y)

function Bits256(x)
    # it loads from the end
    z = unsafe_load(Ptr{Bits256}(pointer_from_objref(x)))

    if sizeof(Bits256) > sizeof(x)
        lzx = leading_zeros(x)
        shift_n = sizeof(Bits256)*8 - sizeof(eltype(x))*8
        z = z << shift_n >> shift_n
    end
    z
end

const u1256_mask = Bits256(2^16-1)

mask16bit(::Type{Bits256}) = u1256_mask

function make_mask(::Type{Bits256})
    x = UInt(2)^(sizeof(UInt)*8) - 1

    yy = unsafe_load(Ptr{Bits256}(pointer_from_objref(x)))
    for i = 2:3
        y = unsafe_load(Ptr{Bits256}(pointer_from_objref(x)))
        # get rid of zeros in front
        y = y << 8*(sizeof(Bits256) - sizeof(UInt)) >> (sizeof(Bits256) - sizeof(UInt))*8
        yy = (yy << 8*sizeof(UInt)) | y
    end
    yy
end

# UInt16(x::Bits256) = unsafe_load(Ptr{UInt16}(pointer_from_objref(x)))

bswap(x::Bits256) = Base.bswap_int(x)

uint_mapping(::Base.Order.ForwardOrdering, x::Bits256) = x

zero(::Type{Bits256}) = Bits256(0)

Int(x::Bits256) = Base.Intrinsics.trunc_int(Int, x)
UInt16(x::Bits256) = Base.Intrinsics.trunc_int(UInt16, x)

# promote_rule(::Type{Bits256}, ::Type{Integer}) = Bits256
(&)(x::Bits256, y::UInt16) = (&)(UInt16(x), y)


text:
## Test models ##
module Test_model

using MultipleTesting
using Base.Test
using Distributions

println(" ** ", "models")
bum = BetaUniformMixtureModel(0.7)
bum0 = BetaUniformMixtureModel(1.0) ## uniform, all H_0
bum1 = BetaUniformMixtureModel(0.0) ## beta, all H_1

x = collect(0:0.01:1)

@test_approx_eq cdf(bum0, x) cdf(Uniform(), x)
@test_approx_eq cdf(bum1, x) cdf(Beta(0.5, 1.0), x)

@test_approx_eq pdf(bum0, x) pdf(Uniform(), x)
@test_approx_eq pdf(bum1, x) pdf(Beta(0.5, 1.0), x)

@test isin(rand(bum, 100), 0.0, 1.0)
@test isin(rand(bum0, 100), 0.0, 1.0)
@test isin(rand(bum1, 100), 0.0, 1.0)

@test_throws MethodError BetaUniformMixtureModel()
@test_throws DomainError BetaUniformMixtureModel(-0.2)
@test_throws DomainError BetaUniformMixtureModel(1.2, 0.2, 2.0)

end


text:
using StatsBase: countmap

export OneClassClassifier, regions, centroid, fit, predict

mutable struct OneClassClassifier{T}
    centers::Vector{T}
    freqs::Vector{Int}
    n::Int
    epsilon::Float64
end

function fit(::Type{OneClassClassifier}, dist::Function, X::AbstractVector{T}, m::Int; centroids=true) where T
    Q = fftclustering(dist, X, m)
    C = X[Q.irefs]
    P = Dict(Q.irefs[i] => i for i in eachindex(Q.irefs))
    freqs = zeros(Int, length(Q.irefs))
    for nn in Q.NN
        freqs[P[first(nn).objID]] += 1
    end

    if centroids
        CC = centroid_correction(dist, X, C)
        OneClassClassifier(CC, freqs, length(X), Q.dmax)
    else
        OneClassClassifier(C, freqs, length(X), Q.dmax)
    end

end

function regions(dist::Function, X, refs::Index)
    I = invindex(dist, X, refs, k=1)
    (freqs=[length(lst) for lst in I], regions=I)
end

function regions(dist::Function, X, refs)
    regions(dist, X, fit(Sequential, refs))
end

function centroid(D)
    sum(D) ./ length(D)
end

function centroid_correction(dist::Function, X, C)
    [centroid(X[lst]) for lst in regions(dist, X, C).regions if length(lst) > 0]
end

function predict(occ::OneClassClassifier{T}, dist::Function, q::T) where T
    seq = fit(Sequential, occ.centers)
    res = search(seq, dist, q, KnnResult(1))
    #1.0 - first(res).dist  / occ.epsilon
    (similarity=max(0.0, 1.0 - first(res).dist  / occ.epsilon), freq=occ.freqs[first(res).objID])
    #occ.freqs[first(res).objID] / occ.n
end


text:
module Stats_module

push!(LOAD_PATH, "../src")

using JLD, Aircraft_module, Constants_module
include("check_status.jl")

export Stats, run_stats

mutable struct Stats
	# airspace_controller::Airspace_Controller
	# cas_ctrl_type::DataType
	# sensor_ctrl::Sensor_Property_Controller
	# dyn_update_rule::Update_Rule

	generated_time::Float64

	# Metrics within simulation wrt time steps:
	t::Array{Int64,1}
	num_ac::Array{Int64,1}
	num_NMAC_t::Array{Int64,1}
	num_NMAC_cum::Array{Int64,1}
	average_flight_time::Array{Float64,1}
	num_ac_arrived::Array{Int64,1}
	num_ac_created::Array{Int64,1}
	traffic_density::Array{Float64,1}
	avg_ac_encountered::Array{Float64,1}
	ac_encountered_per_step_arr::Array{Float64,1}

	# Simulation-wise metrics:
	# Safety:
	total_NMACs::Int64
	total_flight_hours::Int64
	final_traffic_density::Float64

	# Efficiency:
	landed_over_created::Float64
	final_avg_flight_hours::Float64
	total_route_len::Float64
	total_nominal_route_len::Float64
	encounter_distribution::Vector{Float64}

	function Stats(
		generated_time::Float64)
		return new(generated_time,
			[], [], [], [], [], [], [], [], [], [],
			0, 0, 0., 0., 0., 0., 0., zeros(6))
	end
end

function run_stats(All_AC::Array{Aircraft,1}, stats::Stats, t::Int64)
	push!(stats.t, t)
	push!(stats.num_ac, length(All_AC))
	push!(stats.num_NMAC_t, check_NMAC(All_AC))
	push!(stats.num_NMAC_cum, sum(stats.num_NMAC_t))
	avg_flight_time = 0
	if isempty(stats.num_ac_arrived)
		num_ac_arrived = 0
	else
		num_ac_arrived = stats.num_ac_arrived[t-1]
	end
	avg_ac_encountered = 0.
	temp_num_intruders_per_step = 0.
	for ac in All_AC
		avg_flight_time += (t - ac.created_time)
		if check_desti(ac) == true
			num_ac_arrived += 1
		end
		ac.t_to_nominal_t = (t - ac.created_time) / (ac.nominal_route_len / ac.dyn.v)
		ac.d_to_nominal_d = sqrt((ac.dyn.x - ac.x_dest)^2 + (ac.dyn.y - ac.y_dest)^2) / ac.nominal_route_len
		avg_ac_encountered += ac.cumulative_num_ac_encountered / (t - ac.created_time)
		temp_num_intruders_per_step += ac.num_ac_encountered_arr[end]
		if ac.num_ac_encountered_arr[end] > 5
			stats.encounter_distribution[end] += 1
		elseif ac.num_ac_encountered_arr[end] >= 1 && ac.num_ac_encountered_arr[end] <= 5
			stats.encounter_distribution[ac.num_ac_encountered_arr[end]] += 1
		end
	end
	avg_flight_time /= length(All_AC)
	push!(stats.average_flight_time, avg_flight_time)
	push!(stats.num_ac_arrived, num_ac_arrived)
	push!(stats.num_ac_created, num_ac_arrived + length(All_AC))
	push!(stats.traffic_density, length(All_AC) / airspace_dim^2)
	stats.total_flight_hours += length(All_AC) * dt
	stats.total_NMACs = stats.num_NMAC_cum[end]
	stats.final_traffic_density = stats.traffic_density[end]
	stats.landed_over_created = stats.num_ac_arrived[end] / stats.num_ac_created[end]
	stats.final_avg_flight_hours = stats.total_flight_hours / stats.num_ac_created[end]
	avg_ac_encountered /= length(All_AC)
	push!(stats.avg_ac_encountered, avg_ac_encountered)
	temp_num_intruders_per_step /= length(All_AC)
	push!(stats.ac_encountered_per_step_arr, temp_num_intruders_per_step)
end

end # module

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
### CONCRETE TYPE: ITERATIVE PROX EVALUATION
# prox! is computed using CG on a system with matrix A'A + I/(lambda*gamma)
# or AA' + I/(lambda*gamma), according to which matrix is smaller.

using LinearAlgebra

struct LeastSquaresIterative{R <: Real, RC <: RealOrComplex{R}, M, V <: AbstractVector{RC}, O} <: LeastSquares
    A::M # m-by-n operator
    b::V
    Atb::V
    lambda::R
    shape::Symbol
    S::O
    res::Vector{RC} # m-sized buffer
    res2::Vector{RC} # m-sized buffer
    q::Vector{RC} # n-sized buffer
end

is_prox_accurate(f::LeastSquaresIterative) = false

function LeastSquaresIterative(A::M, b::V, lambda::R) where {R <: Real, RC <: RealOrComplex{R}, M, V <: AbstractVector{RC}}
    if size(A, 1) != length(b)
        error("A and b have incompatible dimensions")
    end
    if lambda <= 0
        error("lambda must be positive")
    end
    m, n = size(A)
    if m >= n
        shape = :Tall
        S = AcA(A)
        LeastSquaresIterative{R, RC, M, V, AcA}(A, b, A'*b, lambda, shape, S, zeros(RC, m), [], zeros(RC, n))
    else
        shape = :Fat
        S = AAc(A)
        LeastSquaresIterative{R, RC, M, V, AAc}(A, b, A'*b, lambda, shape, S, zeros(RC, m), zeros(RC, m), zeros(RC, n))
    end
end

function (f::LeastSquaresIterative{R, RC, M, V})(x::AbstractVector{RC}) where {R, RC, M, V}
    mul!(f.res, f.A, x)
    f.res .-= f.b
    return (f.lambda/2)*norm(f.res, 2)^2
end

function prox!(y::AbstractVector{D}, f::LeastSquaresIterative{R, RC, M, V}, x::AbstractVector{D}, gamma::R=R(1)) where {R, RC, M, V, D <: RealOrComplex{R}}
    lamgam = f.lambda*gamma
    f.q .= f.Atb .+ x./lamgam
    # two cases: (1) tall A, (2) fat A
    if f.shape == :Tall
        y .= x
        op = Shift(f.S, R(1)/lamgam)
        IterativeSolvers.cg!(y, op, f.q)
    else # f.shape == :Fat
        # y .= lamgam*(f.q - (f.A'*(f.fact\(f.A*f.q))))
        mul!(f.res, f.A, f.q)
        op = Shift(f.S, R(1)/lamgam)
        IterativeSolvers.cg!(f.res2, op, f.res)
        mul!(y, adjoint(f.A), f.res2)
        y .-= f.q
        y .*= -lamgam
    end
    mul!(f.res, f.A, y)
    f.res .-= f.b
    return (f.lambda/2)*norm(f.res, 2)^2
end

function gradient!(y::AbstractVector{D}, f::LeastSquaresIterative{R, RC, M, V}, x::AbstractVector{D}) where {R, RC, M, V, D <: Union{R, Complex{R}}}
    mul!(f.res, f.A, x)
    f.res .-= f.b
    mul!(y, adjoint(f.A), f.res)
    y .*= f.lambda
    fy = (f.lambda/2)*dot(f.res, f.res)
end

function prox_naive(f::LeastSquaresIterative, x::AbstractVector{D}, gamma::R=R(1)) where {R, D <: RealOrComplex{R}}
    lamgam = f.lambda*gamma
    y = IterativeSolvers.cg(f.A'*f.A + I/lamgam, f.Atb + x/lamgam)
    fy = (f.lambda/2)*norm(f.A*y-f.b)^2
    return y, fy
end


text:
module LagrangianRefFEsTests

using Test
using Gridap.ReferenceFEs
using Gridap.Polynomials
using Gridap.Fields

D = 2
T = Float64
order = 1
prebasis = MonomialBasis{D}(T,order)

polytope = QUAD
x = get_vertex_coordinates(polytope)
dofs = LagrangianDofBasis(T,x)

ndofs = 4

face_own_dofs = [[1],[2],[3],[4],Int[],Int[],Int[],Int[],Int[]]
a = [1]
b = Int[]
face_own_dofs_permutations = [[a],[a],[a],[a],[b,b],[b,b],[b,b],[b,b],fill(b,8)]
face_dofs = [[1],[2],[3],[4],[1,2],[3,4],[1,3],[2,4],[1,2,3,4]]
metadata = nothing

struct MockConformity <: Conformity end

conf = MockConformity()
reffe = GenericRefFE{typeof(conf)}(
  ndofs, polytope, prebasis, dofs,
  conf, metadata, face_dofs)

node_coordinates = Point{2,Float64}[(0,0),(1,0),(0,1),(1,1)]
node_and_comp_to_dof = [1,2,3,4]
face_own_nodes = face_own_dofs
face_own_nodes_permutations = face_own_dofs_permutations
face_nodes = face_dofs

reffe = GenericLagrangianRefFE(reffe,face_nodes)

function ReferenceFEs.get_face_own_nodes(reffe::GenericLagrangianRefFE{MockConformity},::MockConformity)
  face_own_nodes
end

function ReferenceFEs.get_face_own_nodes_permutations(reffe::GenericLagrangianRefFE{MockConformity},::MockConformity)
  face_own_nodes_permutations
end

test_lagrangian_reference_fe(reffe)

@test get_dof_to_node(reffe) == [1, 2, 3, 4]
@test get_own_nodes_permutations(reffe) == get_face_own_nodes_permutations(reffe)[end]

@test get_face_own_dofs(reffe,L2Conformity()) == [[], [], [], [], [], [], [], [], [1, 2, 3, 4]]
@test get_face_own_nodes(reffe,L2Conformity()) == [[], [], [], [], [], [], [], [], [1, 2, 3, 4]]
@test get_face_own_dofs_permutations(reffe,L2Conformity()) == [[[]], [[]], [[]], [[]], [[]], [[]], [[]], [[]], [[1, 2, 3, 4]]]
@test get_face_own_nodes_permutations(reffe,L2Conformity()) == [[[]], [[]], [[]], [[]], [[]], [[]], [[]], [[]], [[1, 2, 3, 4]]]

end # module


text:
@testset "13.roman-to-integer.jl" begin
    @test roman_to_integer("III") == 3
    @test roman_to_integer("IV") == 4
    @test roman_to_integer("IX") == 9
    @test roman_to_integer("LVIII") == 58
    @test roman_to_integer("MCMXCIV") == 1994
end


text:
"""

```
initcells!(x::AbstractVector{T}, box::Box, lc::LinkedCells) where T
```

Function that initializes the linked cells by computing to each cell each atom
belongs and filling up the firstatom and nexatom arrays.

Modifies the data in the lc structure

"""
function initcells!(x::AbstractVector{T}, box::Box, lc::LinkedCells) where T

  # Count the number of boxes and checks if there is a problem with dimensions
  nboxes = box.nc[1]*box.nc[2]*box.nc[3]
  if length(lc.firstatom) < nboxes
    resize!(lc.firstatom,nboxes)
  end

  # Reset arrays
  for i in 1:nboxes
    lc.firstatom[i] = 0
  end
  @. lc.nextatom = 0

  # Initialize cell, firstatom and nexatom
  for iat in eachindex(x)
    ic, jc, kc = icell3D(x[iat],box)
    icell = icell1D(box.nc,ic,jc,kc)
    lc.nextatom[iat] = lc.firstatom[icell]
    lc.firstatom[icell] = iat
  end

  nothing
end


text:
using Test

@time @testset "Unit Tests" begin
    @testset "Small functions" begin

        @info "Bincount"
        @test include("bincount.jl")

        @info "Configuration"
        @test include("configuration.jl")

        @info "dfs"
        @test include("dfs.jl")

        @info "graph"
        @test include("graph.jl")

        @info "forgettingweights"
        @test include("forgettingweights.jl")

        # Needs translation from Python
        # include("adaptive_parzen_normal_orig.jl")

        @info "Trials"
        @test include("trials.jl")

        @info "GMM1"
        @test include("gmm.jl")

        @info "GMM1 Math and QGMM1 Math"
        @test include("gmm_math.jl")

        @info "LGMM1"
        @test include("lgmm.jl")

        @info "Resolve"
        include("resolvenodes.jl")

        @info "Spaces"
        include("spaces.jl")
    end

    @testset "Larger tests" begin
        @info "Basic"
        @test include("basic.jl")

        @info "bjkomer/Squared"
        @test include("bjkomer/squared.jl")

        @info "bjkomer/Function fitting"
        @test include("bjkomer/function_fitting.jl")

        @info "Official Cases"
        @test include("official_cases.jl")

        @info "fmin/Quadratic"
        @test include("fmin/quadratic.jl")

        @info "fmin/Return Inf"
        @test include("fmin/return_inf.jl")

        @info "fmin/Submit points to Trial"
        @test include("fmin/points.jl")

        @info "Silvrback"
        @test include("silvrback.jl")

        @info "Vooban/Basic"
        @test include("vooban/basic.jl")

        @info "Vooban/Find min"
        @test include("vooban/find_min.jl")

        @info "Vooban/Status Fail skip"
        @test include("vooban/status_fail_skip.jl")
    end

    @testset "Samplers" begin
        @info "hp_pchoice"
        @test include("hp.jl")

        @info "LogQuantNormal"
        @test include("logquantnormal.jl")

        @info "QuanLogNormal"
        @test include("quantlognormal.jl")

        @info "LogUniform"
        @test include("loguniform.jl")

        @info "QuantUniform"
        @test include("quantuniform.jl")

        @info "QuantNormal"
        @test include("quantnormal.jl")

        @info "LogQuantUniform"
        @test include("logquantuniform.jl")

        @info "QuantLogUniform"
        @test include("quantloguniform.jl")

        @info "Uniform"
        include("uniform.jl")
    end

    @testset "MLJ" begin
        @info "MLJ Unit tests"
        @test include("MLJ/unit.jl")

        @info "MLJ integration"
        @test include("MLJ/integration.jl")
    end

    @info "API"
    @test include("api.jl")

    # Run this test last so that the print output is just above the test report
    @info "SpacePrint"
    @test include("spaceprint.jl")
end


text:
using Base: AbstractFloat
########################################################################################################
## Details
# This file will create the type structures for the UQ inference problems.
#
#########################################################################################################
abstract type Potential{T} end

mutable struct SNAP{T<:Base.AbstractFloat} <: Potential{T}
    A :: Matrix{T}
    b :: Vector{T}
    β :: Vector{T}
end

abstract type ArbitraryDistribution{T} end

mutable struct PotentialDistribution{T<:AbstractFloat} <: ArbitraryDistribution{T}
    potential::Potential{T}       # potential defined in Potential.jl
    x                             # Container of learnable parameters
    p                             # Container of nonlearnable parameters
    distribution                  # -log likelihood
end

function PotentialDistribution(potential :: Potential{T}, distribution) where T<: AbstractFloat
    pd = PotentialDistribution(potential, p.x, p.q, distribution)
    return pd
end

#########################################################################################################
## SNAP

# Nonlearnable Q
mutable struct SNAPDistribution{T<:AbstractFloat} <: ArbitraryDistribution{T}
    snap::SNAP{T}               # SNAP potential, assume has A, b, β
    x                        # Container of learnable parameters
    p                        # Container of nonlearnable parameters
    distribution             # -log likelihood
end

function SNAPDistribution(snap :: SNAP{T}, Q) where T<: AbstractFloat
    x = (β = snap.β, )
    p = (A = snap.A, b = snap.b, Q = Q)

    distribution(x, p) = 0.5*tr( transpose(p.A*x.β-p.b)*(p.Q\(p.A*x.β-p.b)) + logdet(p.Q))
    pd = SNAPDistribution{T}(snap, x, p, distribution)
    return pd
end

# Learnable Q
mutable struct SNAPQDistribution{T<:AbstractFloat} <: ArbitraryDistribution{T}
    snap::SNAP{T}               # SNAP potential, assume has A, b, β
    x                           # Container of learnable parameters
    p                           # Container of nonlearnable parameters
    distribution                # SNAP likelihood(p)
end

function SNAPQDistribution(snap :: SNAP{T}, Q) where T<: AbstractFloat
    x = (β = snap.β, Q = Q)
    p = (A = snap.A, b = snap.b)
    distribution(x, p) = 0.5*tr( transpose(p.A*x.β-p.b)*(x.Q\(p.A*x.β-p.b)) + logdet(x.Q))
    pd = SNAPQDistribution{T}(snap, x, p, distribution)
    return pd
end


text:
@checked struct LLJITBuilder
    ref::API.LLVMOrcLLJITBuilderRef
end
Base.unsafe_convert(::Type{API.LLVMOrcLLJITBuilderRef}, builder::LLJITBuilder) = builder.ref

@checked mutable struct LLJIT
    ref::API.LLVMOrcLLJITRef
end
Base.unsafe_convert(::Type{API.LLVMOrcLLJITRef}, lljit::LLJIT) = lljit.ref

function LLJITBuilder()
    ref = API.LLVMOrcCreateLLJITBuilder()
    LLJITBuilder(ref)
end

function dispose(builder::LLJITBuilder)
    API.LLVMOrcDisposeLLJITBuilder(builder)
end

function targetmachinebuilder!(builder::LLJITBuilder, tmb::TargetMachineBuilder)
    API.LLVMOrcLLJITBuilderSetJITTargetMachineBuilder(builder, tmb)
end

function linkinglayercreator!(builder::LLJITBuilder, callback, ctx)
    API.LLVMOrcLLJITBuilderSetObjectLinkingLayerCreator(builder, callback, ctx)
end

"""
    LLJIT(::LLJITBuilder)

Creates a LLJIT stack based on the provided builder.

!!! note
    Takes ownership of the provided builder.
"""
function LLJIT(builder::LLJITBuilder)
    ref = Ref{API.LLVMOrcLLJITRef}()
    @check API.LLVMOrcCreateLLJIT(ref, builder)
    LLJIT(ref[])
end

function dispose(lljit::LLJIT)
    API.LLVMOrcDisposeLLJIT(lljit)
end

"""
    LLJIT(;tm::Union{Nothing, TargetMachine})

Use the provided TargetMachine and construct an LLJIT from it.
"""
function LLJIT(; tm::Union{Nothing, TargetMachine} = nothing)
    builder = LLJITBuilder()
    if tm === nothing
        tmb = TargetMachineBuilder()
    else
        tmb = TargetMachineBuilder(tm)
    end
    targetmachinebuilder!(builder, tmb)
    LLJIT(builder)
end

function LLJIT(f::Core.Function, args...; kwargs...)
    lljit = LLJIT(args...; kwargs...)
    try
        f(lljit)
    finally
        dispose(lljit)
    end
end

function triple(lljit::LLJIT)
    cstr = API.LLVMOrcLLJITGetTripleString(lljit)
    Base.unsafe_string(cstr)
end

function apply_datalayout!(lljit::LLJIT, mod::LLVM.Module)
    LLVM.API.LLVMOrcLLJITApplyDataLayout(lljit, mod)
end

function get_prefix(lljit::LLJIT)
    return API.LLVMOrcLLJITGetGlobalPrefix(lljit)
end


text:
#####
##### `Samples`
#####

"""
    Samples(signal::Signal, encoded::Bool, data::AbstractMatrix)

Return a `Samples` instance with the following fields:

  - `signal::Signal`: The `Signal` object that describes the `Samples` instance.

  - `encoded::Bool`: If `true`, the values in `data` are LPCM-encoded as
    prescribed by the `Samples` instance's `signal`. If `false`, the values in
    `data` have been decoded into the `signal`'s canonical units.
  - `data::AbstractMatrix`: A matrix of sample data. The `i` th row of the matrix
    corresponds to the `i`th channel in `signal.channel_names`, while the `j`th
    column corresponds to the `j`th multichannel sample.

Note that `getindex` and `view` are defined on `Samples` to accept normal integer
indices, but also accept channel names for row indices and [`TimeSpan`](@ref)
values for column indices; see `Onda/examples/tour.jl` for a comprehensive
set of indexing examples.

See also: [`encode`](@ref), [`encode!`](@ref), [`decode`](@ref), [`decode!`](@ref)
"""
struct Samples{D<:AbstractMatrix}
    signal::Signal
    encoded::Bool
    data::D
    function Samples(signal::Signal, encoded::Bool, data::AbstractMatrix)
        n_channels = channel_count(signal)
        n_rows = size(data, 1)
        if n_channels != n_rows
            throw(DimensionMismatch("number of channels in signal ($n_channels) " *
                                    "does not match number of rows in data matrix " *
                                    "($n_rows)"))
        end
        if encoded && !(eltype(data) <: signal.sample_type)
            throw(ArgumentError("signal and encoded data matrix have mismatched element types"))
        end
        return new{typeof(data)}(signal, encoded, data)
    end
end

for f in (:getindex, :view)
    @eval begin
        import Base: $f  # required for deprecating methods for functions extended from Base
        @inline function Base.$f(samples::Samples, rows, columns)
            rows = row_arguments(samples, rows)
            columns = column_arguments(samples, columns)
            signal = rows isa Colon ? samples.signal :
                     signal_from_template(samples.signal;
                                          channel_names=samples.signal.channel_names[rows])
            return Samples(signal, samples.encoded, $f(samples.data, rows, columns))
        end
        Base.@deprecate $f(samples::Samples, columns) $f(samples, :, columns)
    end
end

_rangify(i) = i
_rangify(i::Integer) = i:i

function _indices_fallback(f, samples::Samples, i::Union{Colon,AbstractRange,Integer})
    return i
end
_indices_fallback(f, samples::Samples, args) = map(x -> f(samples, x), args)

row_arguments(samples::Samples, args) = _rangify(_row_arguments(samples, args))
function _row_arguments(samples::Samples, args)
    return _indices_fallback(_row_arguments, samples, args)
end
_row_arguments(samples::Samples, name::Symbol) = channel(samples, name)

function column_arguments(samples::Samples, args)
    return _rangify(_column_arguments(samples, args))
end
function _column_arguments(samples::Samples, args)
    return _indices_fallback(_column_arguments, samples, args)
end
function _column_arguments(samples::Samples, t::Period)
    return _column_arguments(samples, TimeSpan(t))
end

function _column_arguments(samples::Samples, span::AbstractTimeSpan)
    return index_from_time(samples.signal.sample_rate, span)
end

"""
    channel(samples::Samples, name::Symbol)

Return `channel(samples.signal, name)`.

This function is useful for indexing rows of `samples.data` by channel names.
"""
channel(samples::Samples, name::Symbol) = channel(samples.signal, name)

"""
    channel(samples::Samples, i::Integer)

Return `channel(samples.signal, i)`.
"""
channel(samples::Samples, i::Integer) = channel(samples.signal, i)

"""
    duration(samples::Samples)

Returns the `Nanosecond` value for which `samples[TimeSpan(0, duration(samples))] == samples.data`.

!!! warning

    `duration(samples)` is not generally equivalent to `duration(samples.signal)`;
    the former is the duration of the entire original signal in the context of its
    parent recording, whereas the latter is the actual duration of `samples.data`
    given `samples.signal.sample_rate`.
"""
function duration(samples::Samples)
    return time_from_index(samples.signal.sample_rate, size(samples.data, 2) + 1)
end

"""
    channel_count(samples::Samples)

Return `channel_count(samples.signal)`.
"""
channel_count(samples::Samples) = channel_count(samples.signal)

"""
    sample_count(samples::Samples)

Return the number of multichannel samples in `samples` (i.e. `size(samples.data, 2)`)

!!! warning

    `sample_count(samples)` is not generally equivalent to `sample_count(samples.signal)`;
    the former is the sample count of the entire original signal in the context of its parent
    recording, whereas the latter is actual number of multichannel samples in `samples.data`.
"""
sample_count(samples::Samples) = size(samples.data, 2)

#####
##### encoding utilities
#####

const VALID_SAMPLE_TYPE_UNION = Union{Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64}

function encode_sample(::Type{S}, resolution_in_unit, offset_in_unit, sample_in_unit,
                       noise=zero(sample_in_unit)) where {S<:VALID_SAMPLE_TYPE_UNION}
    sample_in_unit += noise
    isnan(sample_in_unit) && return typemax(S)
    from_unit = clamp((sample_in_unit - offset_in_unit) / resolution_in_unit, typemin(S),
                      typemax(S))
    return round(S, from_unit)
end

function dither_noise!(rng::AbstractRNG, storage, step)
    rand!(rng, storage)
    broadcast!(_dither_noise, storage, storage, step + step, step)
    return storage
end

dither_noise!(storage, step) = dither_noise!(Random.GLOBAL_RNG, storage, step)

function _dither_noise(x, range, step)
    rs = range * x
    return if rs < step
        return sqrt(rs * step) - step
    else
        return step - sqrt(range * (1 - x) * step)
    end
end

#####
##### `encode`/`encode!`
#####

"""
    encode(sample_type::DataType, sample_resolution_in_unit, sample_offset_in_unit,
           samples, dither_storage=nothing)

Return a copy of `samples` quantized according to `sample_type`, `sample_resolution_in_unit`,
and `sample_offset_in_unit`. `sample_type` must be a concrete subtype of `Onda.VALID_SAMPLE_TYPE_UNION`.
Quantization of an individual sample `s` is performed via:

    round(S, (s - sample_offset_in_unit) / sample_resolution_in_unit)

with additional special casing to clip values exceeding the encoding's dynamic range.

If `dither_storage isa Nothing`, no dithering is applied before quantization.

If `dither_storage isa Missing`, dither storage is allocated automatically and
triangular dithering is applied to the signal prior to quantization.

Otherwise, `dither_storage` must be a container of similar shape and type to
`samples`. This container is then used to store the random noise needed for the
triangular dithering process, which is applied to the signal prior to quantization.
"""
function encode(::Type{S}, sample_resolution_in_unit, sample_offset_in_unit, samples,
                dither_storage=nothing) where {S}
    return encode!(similar(samples, S), S, sample_resolution_in_unit, sample_offset_in_unit,
                   samples, dither_storage)
end

"""
    encode!(result_storage, sample_type::DataType, sample_resolution_in_unit,
            sample_offset_in_unit, samples, dither_storage=nothing)
    encode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit,
            samples, dither_storage=nothing)

Similar to `encode(sample_type, sample_resolution_in_unit, sample_offset_in_unit, samples, dither_storage)`,
but write encoded values to `result_storage` rather than allocating new storage.

`sample_type` defaults to `eltype(result_storage)` if it is not provided.
"""
function encode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit, samples,
                 dither_storage=nothing)
    return encode!(result_storage, eltype(result_storage), sample_resolution_in_unit,
                   sample_offset_in_unit, samples; dither_storage=nothing)
end

function encode!(result_storage, ::Type{S}, sample_resolution_in_unit,
                 sample_offset_in_unit, samples, dither_storage=nothing) where {S}
    if dither_storage isa Nothing
        broadcast!(encode_sample, result_storage, S, sample_resolution_in_unit,
                   sample_offset_in_unit, samples)
    else
        if dither_storage isa Missing
            dither_storage = similar(samples)
        elseif size(dither_storage) != size(samples)
            throw(DimensionMismatch("dithering storage container does not match shape of samples"))
        end
        dither_noise!(dither_storage, sample_resolution_in_unit)
        broadcast!(encode_sample, result_storage, S, sample_resolution_in_unit,
                   sample_offset_in_unit, samples, dither_storage)
    end
    return result_storage
end

"""
    encode(samples::Samples, dither_storage=nothing)

If `samples.encoded` is `false`, return a `Samples` instance that wraps:

    encode(samples.signal.sample_type,
           samples.signal.sample_resolution_in_unit,
           samples.signal.sample_offset_in_unit,
           samples.data, dither_storage)

If `samples.encoded` is `true`, this function is the identity.
"""
function encode(samples::Samples, dither_storage=nothing)
    samples.encoded && return samples
    data = encode(samples.signal.sample_type, samples.signal.sample_resolution_in_unit,
                  samples.signal.sample_offset_in_unit, samples.data, dither_storage)
    return Samples(samples.signal, true, data)
end

"""
    encode!(result_storage, samples::Samples, dither_storage=nothing)

If `samples.encoded` is `false`, return a `Samples` instance that wraps:

    encode!(result_storage,
            samples.signal.sample_type,
            samples.signal.sample_resolution_in_unit,
            samples.signal.sample_offset_in_unit,
            samples.data, dither_storage)`.

If `samples.encoded` is `true`, return a `Samples` instance that wraps
`copyto!(result_storage, samples.data)`.
"""
function encode!(result_storage, samples::Samples, dither_storage=nothing)
    if samples.encoded
        copyto!(result_storage, samples.data)
        return Samples(samples.signal, samples.encoded, result_storage)
    end
    encode!(result_storage, samples.signal.sample_type,
            samples.signal.sample_resolution_in_unit, samples.signal.sample_offset_in_unit,
            samples.data, dither_storage)
    return Samples(samples.signal, true, result_storage)
end

#####
##### `decode`/`decode!`
#####

"""
    decode(sample_resolution_in_unit, sample_offset_in_unit, samples)

Return `sample_resolution_in_unit .* samples .+ sample_offset_in_unit`
"""
function decode(sample_resolution_in_unit, sample_offset_in_unit, samples)
    return sample_resolution_in_unit .* samples .+ sample_offset_in_unit
end

"""
    decode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit, samples)

Similar to `decode(sample_resolution_in_unit, sample_offset_in_unit, samples)`, but
write decoded values to `result_storage` rather than allocating new storage.
"""
function decode!(result_storage, sample_resolution_in_unit, sample_offset_in_unit, samples)
    f = x -> sample_resolution_in_unit * x + sample_offset_in_unit
    return broadcast!(f, result_storage, samples)
end

"""
    decode(samples::Samples)

If `samples.encoded` is `true`, return a `Samples` instance that wraps
`decode(samples.signal.sample_resolution_in_unit, samples.signal.sample_offset_in_unit, samples.data)`.

If `samples.encoded` is `false`, this function is the identity.
"""
function decode(samples::Samples)
    samples.encoded || return samples
    data = decode(samples.signal.sample_resolution_in_unit,
                  samples.signal.sample_offset_in_unit, samples.data)
    return Samples(samples.signal, false, data)
end

"""
    decode!(result_storage, samples::Samples)

If `samples.encoded` is `true`, return a `Samples` instance that wraps
`decode!(result_storage, samples.signal.sample_resolution_in_unit, samples.signal.sample_offset_in_unit, samples.data)`.

If `samples.encoded` is `false`, return a `Samples` instance that wraps
`copyto!(result_storage, samples.data)`.
"""
function decode!(result_storage, samples::Samples)
    if samples.encoded
        decode!(result_storage, samples.signal.sample_resolution_in_unit,
                samples.signal.sample_offset_in_unit, samples.data)
        return Samples(samples.signal, false, result_storage)
    end
    copyto!(result_storage, samples.data)
    return Samples(samples.signal, samples.encoded, result_storage)
end

#####
##### load/store
#####

function load_samples(file_path::AbstractString, signal::Signal;
                      serializer=serializer(signal))
    return Samples(signal, true, deserialize_lpcm(read(file_path), serializer))
end

function load_samples(file_path::AbstractString, signal::Signal, span::AbstractTimeSpan;
                      serializer=serializer(signal))
    sample_range = index_from_time(signal.sample_rate, span)
    offset, n = first(sample_range) - 1, length(sample_range)
    data = open(io -> deserialize_lpcm(io, serializer, offset, n), file_path, "r")
    return Samples(signal, true, data)
end

function store_samples!(file_path::AbstractString, samples::Samples; overwrite::Bool=true,
                        serializer=serializer(samples.signal))
    overwrite || (isfile(file_path) &&
     error("overwrite disabled but file path already exists: $(file_path)"))
    samples = encode(samples)
    open(io -> serialize_lpcm(io, samples.data, serializer), file_path, "w")
    return file_path
end