text:
__precompile__()

module ArbFloats

import Base: STDOUT,
    hash, convert, promote_rule, isa,
    string, show, showcompact, showall, parse,
    finalizer, decompose, precision, setprecision,
    serialize, deserialize,
    typemin, typemax, realmin, realmax,
    copy, deepcopy,
    zero, one, isinteger,
    ldexp, frexp, eps,
    isequal, isless, (==),(!=),(<),(<=),(>=),(>), contains,
    min, max, minmax,
    typemax, typemin, realmax, realmin,
    float, nextfloat, prevfloat,
    isnan, isinf, isfinite, issubnormal,
    signbit, sign, flipsign, copysign, abs, abs2,
    (+),(-),(*),(/),(\),(%),(^), inv, sqrt, hypot,
    (.+),(.-),(.*),(./),
    trunc, round, ceil, floor,
    fld, cld, div, mod, rem, divrem, fldmod,
    muladd, fma,
    exp, expm1, log, log1p, log2, log10,
    sin, cos, tan, csc, sec, cot, asin, acos, atan, atan2,
    sinh, cosh, tanh, csch, sech, coth, asinh, acosh, atanh,
    sinc, gamma, lgamma, digamma, zeta, factorial,
    in, union, intersect,
    BigInt, BigFloat, Float64, Float32, Int128, Int64, Int32, Rational,
    Cint

export ArbFloat,      # co-matched decimal rounding, n | round(hi,n,10) == round(lo,n,10)
       @ArbFloat,     # converts string form of argument, precision is optional first arg
       setprecisionAugmented,
       simeq, nsime, prec, preceq, succ, succeq, # non-strict total ordering comparisons
       (≃), (≄), (≺), (⪯), (≻), (⪰),           #    matched binary operators
       upperbound, lowerbound, bounds,
       midpoint, radius, midpoint_radius,
       bounding_midpoint, bounding_radius, bounding_midpoint_radius,
       stringsmall, stringcompact, stringmedium, stringlarge, stringall,
       stringsmall_pm, stringcompact_pm, string_pm,
       stringlarge_pm, stringall_pm,
       showsmall, showcompact, showlarge, showallcompact, showpretty,
       showsmall_pm, showcompact_pm, show_pm,
       showlarge_pm, showall_pm,
       stringpretty, smartvalue, smartstring, showsmart,
       two, three, four, copymidpoint, copyradius, deepcopyradius,
       get_emax, get_emin, bounded, boundedrange,
       decompose, isexact, notexact,
       isposinf, isneginf,
       notnan, notinf, notposinf, notneginf, notfinite,
       iszero, notzero, nonzero, isone, notone, notinteger,
       ispositive, notpositive, isnegative, notnegative,
       includesAnInteger, excludesIntegers, includesZero, excludesZero,
       includesPositive, excludesPositive, includesNegative, excludesNegative,
       includesNonpositive,  includesNonnegative,
       areequal, notequal, approxeq, (≊),
       narrow, overlap, donotoverlap,
       contains, iscontainedby, doesnotcontain, isnotcontainedby,
       absz, invsqrt, pow, root, tanpi, cotpi, logbase, sincos, sincospi, sinhcosh,
       doublefactorial, risingfactorial, rgamma, agm, polylog,
       relativeError, relativeAccuracy, midpointPrecision, trimmed,
       PI,SQRTPI,LOG2,LOG10,EXP1,EULER,CATALAN,KHINCHIN,GLAISHER,APERY, # constants
       get_midpoint_digits_shown, get_radius_digits_shown,  # some interface control
       set_midpoint_digits_shown, set_radius_digits_shown,
       isolate_nonnegative_content, isolate_positive_content, # for interval algorithms
       force_nonnegative_content, force_positive_conent

NotImplemented(info::AbstractString="") = error(string("this is not implemented\n\t",info,"\n"))

include("support/NemoLibs.jl")                 # for precompiled libraries
include("support/ReadableNumbers.jl")          # digit subsequence separators

include("type/ArbCstructs.jl")
include("type/MagFloat.jl")
include("type/ArfFloat.jl")
include("type/ArbFloat.jl")
include("type/ArbInterval.jl")

include("basics/primitive.jl")
include("basics/IEEEfp.jl")

include("basics/predicates.jl")
include("basics/convert.jl")
include("basics/compare.jl")

include("basics/string.jl")
include("basics/smartstring.jl")
include("basics/show.jl")
include("basics/serialize.jl")

include("math/arith.jl")
include("math/round.jl")
include("math/elementary.jl")
include("math/constants.jl")
include("math/special.jl")

include("math/rounding.jl")

include("math/arrayops.jl")

#=
# precision is significand precision, significand_bits(FloatNN) + 1, for the hidden bit
typealias ArbFloat16  ArbFloat{ 11}  # read   2 ? 3 or fewer decimal digits to write the same digits ( 16bit Float)
typealias ArbFloat32  ArbFloat{ 24}  # read   6 ? 7 or fewer decimal digits to write the same digits ( 32bit Float)
typealias ArbFloat64  ArbFloat{ 53}  # read  15 ?15 or fewer decimal digits to write the same digits ( 64bit Float)
typealias ArbFloat128 ArbFloat{113}  # read  33 ?34 or fewer decimal digits to write the same digits (128bit Float)
typealias ArbFloat256 ArbFloat{237}  # read  71 ?71 or fewer decimal digits to write the same digits (256bit Float)
typealias ArbFloat512 ArbFloat{496}  # read 148?149 or fewer decimal digits to write the same digits (512bit Float)
=#

end # ArbFloats


text:
struct ConstPotential{T} <: Potential end
ConstPotential() = ConstPotential{Float64}()

UU(phi,  ::ConstPotential) = phi
UUp(phi, ::ConstPotential{T}) where {T} = one(T)

parameters(potential::ConstPotential) = ()

struct SquarePotential{T} <: Potential end
SquarePotential() = SquarePotential{Float64}()

UU(phi,  ::SquarePotential) = -1 + phi*phi / 2
UUp(phi, ::SquarePotential) = phi

parameters(potential::SquarePotential) = ()

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule librttopo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("librttopo")
JLLWrappers.@generate_main_file("librttopo", UUID("4a3ffa5b-425c-5e62-932c-e04bcaef15ee"))
end  # module librttopo_jll


text:
using Flux

"""
    vgg19

returns the vgg19 model
"""
vgg19 = Chain(
  Conv((3, 3), 6 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  MaxPool((2,2)),
  Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  MaxPool((2,2)),
  Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(8192, 4096, relu),
  Dropout(0.3),
  Dense(4096, 4096, relu),
  Dropout(0.3),
  Dense(4096, 30)
)

vgg19_large = Chain(
  Conv((3, 3), 12 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  MaxPool((2,2)),
  Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  MaxPool((2,2)),
  Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  MaxPool((2,2)),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(4608, 4096, relu),
  Dropout(0.3),
  Dense(4096, 4096, relu),
  Dropout(0.3),
  Dense(4096, 30)
)


text:
#################### Geweke Diagnostic ####################

"""
    gewekediag(x::Vector{<:Real}; first, last, etype)
    gewekediag(chains::Chains; sections, first, last, etype, kwargs...)

Geweke diagnostic.
"""
function gewekediag(x::Vector{<:Real}; first::Real=0.1, last::Real=0.5,
                             etype=:imse, args...)
  if !(0.0 < first < 1.0)
    throw(ArgumentError("first is not in (0, 1)"))
  elseif !(0.0 < last < 1.0)
    throw(ArgumentError("last is not in (0, 1)"))
  elseif first + last > 1.0
    throw(ArgumentError("first and last proportions overlap"))
  end
  n = length(x)
  x1 = x[1:round(Int, first * n)]
  x2 = x[round(Int, n - last * n + 1):n]
  z = (mean(x1) - mean(x2)) /
      sqrt(mcse(x1, etype; args...)^2 + mcse(x2, etype; args...)^2)
  [z, 1 - erf(abs(z) / sqrt(2))]
end

function gewekediag(
    chains::Chains;
    sections = _default_sections(chains),
    first::Real = 0.1,
    last::Real = 0.5,
    etype = :imse,
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    _, p, m = size(_chains.value)
    diags = [Array{Float64}(undef, p, 2) for _ in 1:m]
    for j in 1:p, k in 1:m
        diags[k][j,:] = gewekediag(
                            collect(skipmissing(_chains.value[:, j, k])),
                            first=first,
                            last=last,
                            etype=etype;
                            kwargs...
                        )
    end

    # Obtain names of parameters.
    names_of_params = names(_chains)

    # Compute data frames.
    vector_of_df = [
        ChainDataFrame(
            "Geweke Diagnostic - Chain $i",
            (parameters = names_of_params, zscore = d[:, 1], pvalue = d[:, 2])
        )
        for (i, d) in enumerate(diags)
    ]

    return vector_of_df
end


text:
using Push
using Base.Test

cfg_path = joinpath(dirname(@__FILE__), "configuration/integer.cfg")
cfg = Push.load_configuration(cfg_path)

# INTEGER.*
s = Push.run("(1 INTEGER.*)", cfg)
@test s.integer == [1]
s = Push.run("(7 3 2 10 INTEGER.*)", cfg)
@test s.integer == [7, 3, 20]

# INTEGER.+
s = Push.run("(1 INTEGER.+)", cfg)
@test s.integer == [1]
s = Push.run("(7 3 2 10 INTEGER.+)", cfg)
@test s.integer == [7, 3, 12]
s = Push.run("(7 3 2 -10 INTEGER.+)", cfg)
@test s.integer == [7, 3, -8]

# INTEGER.-
s = Push.run("(1 INTEGER.-)", cfg)
@test s.integer == [1]
s = Push.run("(7 3 10 2 INTEGER.-)", cfg)
@test s.integer == [7, 3, 8]
s = Push.run("(7 3 2 10 INTEGER.-)", cfg)
@test s.integer == [7, 3, -8]

# INTEGER./
s = Push.run("(1 INTEGER./)", cfg)
@test s.integer == [1]
s = Push.run("(7 3 10 2 INTEGER./)", cfg)
@test s.integer == [7, 3, 5]
s = Push.run("(7 3 -10 4 INTEGER./)", cfg)
@test s.integer == [7, 3, -2]
s = Push.run("(7 3 -10 0 INTEGER./)", cfg)
@test s.integer == [7, 3, -10, 0]

# INTEGER.<
s = Push.run("(1 INTEGER.<)", cfg)
@test s.integer == [1] && isempty(s.boolean)
s = Push.run("(100 900 2 3 INTEGER.<)", cfg)
@test s.integer == [100, 900] && s.boolean == [true]
s = Push.run("(100 900 5 -10 INTEGER.<)", cfg)
@test s.integer == [100, 900] && s.boolean == [false]
s = Push.run("(100 900 5 5 INTEGER.<)", cfg)
@test s.integer == [100, 900] && s.boolean == [false]

# INTEGER.=
s = Push.run("(1 INTEGER.=)", cfg)
@test s.integer == [1] && isempty(s.boolean)
s = Push.run("(100 900 -10 90 INTEGER.=)", cfg)
@test s.integer == [100, 900] && s.boolean == [false]
s = Push.run("(100 900 -10 -10 INTEGER.=)", cfg)
@test s.integer == [100, 900] && s.boolean == [true]

# INTEGER.>
s = Push.run("(1 INTEGER.>)", cfg)
@test s.integer == [1] && isempty(s.boolean)
s = Push.run("(100 900 2 3 INTEGER.>)", cfg)
@test s.integer == [100, 900] && s.boolean == [false]
s = Push.run("(100 900 5 -10 INTEGER.>)", cfg)
@test s.integer == [100, 900] && s.boolean == [true]
s = Push.run("(100 900 5 5 INTEGER.>)", cfg)
@test s.integer == [100, 900] && s.boolean == [false]

# INTEGER.DUP
s = Push.run("(10 20 30 INTEGER.DUP)", cfg)
@test s.integer == [10, 20, 30, 30]

# INTEGER.FLUSH
s = Push.run("(10 20 30 INTEGER.FLUSH)", cfg)
@test s.integer == []

# INTEGER.FROMBOOLEAN
s = Push.run("(10 20 TRUE FALSE INTEGER.FROMBOOLEAN)", cfg)
@test s.integer == [10, 20, 0] && s.boolean == [true]
s = Push.run("(10 20 FALSE TRUE INTEGER.FROMBOOLEAN)", cfg)
@test s.integer == [10, 20, 1] && s.boolean == [false]

# INTEGER.FROMFLOAT
s = Push.run("(1.0 2.0 3.0 INTEGER.FROMFLOAT)", cfg)
@test s.integer == [3] && s.float == [1.0, 2.0]

# INTEGER.MAX
s = Push.run("(10 10 INTEGER.MAX)", cfg)
@test s.integer == [10]
s = Push.run("(20 10 INTEGER.MAX)", cfg)
@test s.integer == [20]
s = Push.run("(10 20 INTEGER.MAX)", cfg)
@test s.integer == [20]

# INTEGER.MIN
s = Push.run("(10 10 INTEGER.MIN)", cfg)
@test s.integer == [10]
s = Push.run("(20 10 INTEGER.MIN)", cfg)
@test s.integer == [10]
s = Push.run("(10 20 INTEGER.MIN)", cfg)
@test s.integer == [10]

# INTEGER.POP
s = Push.run("(10 20 30 INTEGER.POP)", cfg)
@test s.integer == [10, 20]

# INTEGER.ROT
s = Push.run("(1 INTEGER.ROT)", cfg)
@test s.integer == [1]
s = Push.run("(0 10 20 30 INTEGER.ROT)", cfg)
@test s.integer == [0, 30, 10, 20]

# INTEGER.SHOVE
s = Push.run("(1 INTEGER.SHOVE)", cfg)
@test s.integer == [1]
s = Push.run("(2 0 INTEGER.SHOVE)", cfg)
@test s.integer == [2]
s = Push.run("(2 9000 INTEGER.SHOVE)", cfg)
@test s.integer == [2]
s = Push.run("(1 2 3 4 5 6 7 8 900 -20 INTEGER.SHOVE)", cfg)
@test s.integer == [1,2,3,4,5,6,7,8,900]
s = Push.run("(1 2 3 4 5 6 7 8 900 0 INTEGER.SHOVE)", cfg)
@test s.integer == [1,2,3,4,5,6,7,8,900]
s = Push.run("(1 2 3 4 5 6 7 8 900 1 INTEGER.SHOVE)", cfg)
@test s.integer == [1,2,3,4,5,6,7,900,8]
s = Push.run("(1 2 3 4 5 6 7 8 900 2 INTEGER.SHOVE)", cfg)
@test s.integer == [1,2,3,4,5,6,900,7,8]
s = Push.run("(1 2 3 4 5 6 7 8 900 9 INTEGER.SHOVE)", cfg)
@test s.integer == [900,1,2,3,4,5,6,7,8]
s = Push.run("(1 2 3 4 5 6 7 8 900 1000 INTEGER.SHOVE)", cfg)
@test s.integer == [900,1,2,3,4,5,6,7,8]

# INTEGER.STACKDEPTH
s = Push.run("(INTEGER.STACKDEPTH)", cfg)
@test s.integer == [0]
s = Push.run("(10 20 30 40 INTEGER.STACKDEPTH)", cfg)
@test s.integer == [10,20,30,40, 4]

# INTEGER.SWAP
s = Push.run("(1 INTEGER.SWAP)", cfg)
@test s.integer == [1]
s = Push.run("(1 2 3 4 5 6 INTEGER.SWAP)", cfg)
@test s.integer == [1,2,3,4,6,5]

# INTEGER.YANK
s = Push.run("(1 INTEGER.YANK)", cfg)
@test s.integer == [1]
s = Push.run("(10 0 INTEGER.YANK)", cfg)
@test s.integer == [10]
s = Push.run("(10 20 30 40 50 -90 INTEGER.YANK)", cfg)
@test s.integer == [10, 20, 30, 40, 50]
s = Push.run("(10 20 30 40 50 0 INTEGER.YANK)", cfg)
@test s.integer == [10, 20, 30, 40, 50]
s = Push.run("(10 20 30 40 50 1 INTEGER.YANK)", cfg)
@test s.integer == [10, 20, 30, 50, 40]
s = Push.run("(10 20 30 40 50 2 INTEGER.YANK)", cfg)
@test s.integer == [10, 20, 40, 50, 30]
s = Push.run("(10 20 30 40 50 5 INTEGER.YANK)", cfg)
@test s.integer == [20, 30, 40, 50, 10]
s = Push.run("(10 20 30 40 50 987 INTEGER.YANK)", cfg)
@test s.integer == [20, 30, 40, 50, 10]

# INTEGER.YANKDUP
s = Push.run("(1 INTEGER.YANKDUP)", cfg)
@test s.integer == [1]
s = Push.run("(10 20 30 40 50 0 INTEGER.YANKDUP)", cfg)
@test s.integer == [10, 20, 30, 40, 50, 50]
s = Push.run("(10 20 30 40 50 1 INTEGER.YANKDUP)", cfg)
@test s.integer == [10, 20, 30, 40, 50, 40]
s = Push.run("(10 20 30 40 50 2 INTEGER.YANKDUP)", cfg)
@test s.integer == [10, 20, 30, 40, 50, 30]
s = Push.run("(10 20 30 40 50 3 INTEGER.YANKDUP)", cfg)
@test s.integer == [10, 20, 30, 40, 50, 20]
s = Push.run("(10 20 30 40 50 4 INTEGER.YANKDUP)", cfg)
@test s.integer == [10, 20, 30, 40, 50, 10]
s = Push.run("(10 20 30 40 50 987 INTEGER.YANKDUP)", cfg)
@test s.integer == [10, 20, 30, 40, 50, 10]
s = Push.run("(10 20 30 40 50 -987 INTEGER.YANKDUP)", cfg)
@test s.integer == [10, 20, 30, 40, 50, 50]

# INTEGER.DEFINE
s = Push.run("(X INTEGER.DEFINE X)", cfg)
@test isempty(s.integer) && s.name == [:X, :X]
s = Push.run("(3 X INTEGER.DEFINE X)", cfg)
@test isempty(s.name) && s.integer == [3]

# INTEGER.RAND
s = Push.run("(INTEGER.RAND)", cfg)
@test length(s.integer) == 1

# INTEGER.%


text:
abstract type AbstractEmbedding end

struct Embedding{𝒯} <: AbstractEmbedding
    manifold::𝒯
end

function CartesianEmbedding(; manifold)
    return Embedding(manifold)
end

function (embedding::CartesianEmbedding{Interval})(x)
    a = embedding.manifold.a
    b = embedding.manifold.b
    @assert a ≤ x ≤ b

    return [x]
end

function (embedding::CartesianEmbedding{Rectangle})(x, y)
    a = embedding.manifold.a
    b = embedding.manifold.b
    c = embedding.manifold.c
    d = embedding.manifold.d
    @assert a ≤ x ≤ b
    @assert c ≤ y ≤ d

    return [x, y]
end

function (embedding::SphericalEmbedding{Sphere})(λ, ϕ)
    @assert -π ≤ λ ≤ π
    @assert -π/2 ≤ ϕ ≤ π/2
    a = embedding.manifold.radius

    return [a*cos(λ)*sin(ϕ), ]
end


text:
@testset "Signature" begin
    @test Signature("+++") == Signature(3, 0)
    @test Signature("++-") == Signature(2, 1)
    @test Signature("+-𝟎") == Signature(1, 1, 1)

    @test !is_degenerate(ℝ⁴)
    @test !is_degenerate(spacetime)
    @test is_degenerate(sig_111)

    @test metric(SIGNATURE, Val(1), Val(1)) == 1
    @test metric(SIGNATURE, Val(1), Val(2)) == 0
    @test metric(SIGNATURE, Val(3), Val(3)) == -1
end


text:
##< work
work!( world::World, d::Dog, dt ) = work!( world, Val(d.state), d, dt )

function work!( world, ::Val{ACTIVE}, dog::Dog, dt )
    i,j = getSheepCoords( world, dog.position )
    n_sheep = getNSheep( world.sheeps, i, j )
    n_sheep > 0 && @show n_sheep
    n_sheep == 0 && return dog
    if rand() < n_sheep * world.business_rate * dt
        # pick up a sheep
        setNSheep!( world.sheeps, i, j, n_sheep - 1 )
        return Dog( dog, state = OCCUPIED )
    end # if
    return dog
end # function

function work!( world, ::Val{OCCUPIED}, dog::Dog, dt )
    i,j = getSheepCoords( world, dog.position )
    n_sheep = getNSheep( world.sheeps, i, j )
    n_sheep == 0 && return dog
    if rand() < n_sheep * world.freedom_rate * dt
        # drop a sheep
        setNSheep!( world.sheeps, i, j, n_sheep + 1 )
        return Dog( dog;
                state = SLEEPY,
                sleepyness = rand( Exponential( world.meanSleepyness ) ) ÷ dt
            )
    end # if
    return dog
end # function

function work!( world, ::Val{SLEEPY}, dog::Dog, dt )
    if dog.sleepyness == 0
        return Dog( dog; state = ACTIVE )
    end # if
    return Dog( dog; sleepyness = dog.sleepyness - 1 )
end # function
##> work


text:
"""
    LinearOettliPrager <: AbstractDirectSolver

Type for the OettliPrager solver of the interval linear system ``Ax=b``. The solver first
converts the system of interval equalities into a system of real inequalities using
Oettli-Präger theorem [[OET64]](@ref) and then finds the feasible set by solving a LP
problem in each orthant using `LazySets.jl`.

### Notes
- You need to import `LazySets.jl` to use this functionality.
- An object of type `LinearOettliPrager` is a function with methods

        (op::LinearOettliPrager)(A::AbstractMatrix{T},
                                 b::AbstractVector{T}) where {T<:Interval}

    #### Input
    - `A`   -- N×N interval matrix
    - `b`   -- interval vector of length N

### Examples

```julia-repl
julia> A = [2..4 -2..1;-1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> b = [-2..2, -2..2]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-2, 2]

julia> polytopes = solve(A, b, LinearOettliPrager());

julia> typeof(polytopes)
Vector{HPolytope{Float64, SparseArrays.SparseVector{Float64, Int64}}}
```
"""
struct LinearOettliPrager <: AbstractDirectSolver end

"""
    NonLinearOettliPrager <: AbstractIterativeSolver

Type for the OettliPrager solver of the interval linear system ``Ax=b``. The solver first
converts the system of interval equalities into a system of real inequalities using
Oettli-Präger theorem [[OET64]](@ref) and then finds the feasible set using
the forward-backward contractor method [[JAU14]](@ref) implemented in
`IntervalConstraintProgramming.jl`.

### Fields
- `tol` -- tolerance for the paving, default 0.01.

### Notes
- You need to import `IntervalConstraintProgramming.jl` to use this functionality.
- An object of type `NonLinearOettliPrager` is a function with methods

        (op::NonLinearOettliPrager)(A::AbstractMatrix{T},
                                    b::AbstractVector{T},
                                    [X]::AbstractVector{T}=enclose(A, b)) where {T<:Interval}

        (op::NonLinearOettliPrager)(A::AbstractMatrix{T},
                                    b::AbstractVector{T},
                                    X::IntervalBox) where {T<:Interval}

    #### Input
    - `A`   -- N×N interval matrix
    - `b`   -- interval vector of length N
    - `X`   -- (optional) initial enclosure for the solution of ``Ax = b``. If not given,
               it is automatically computed using [`enclose`](@ref enclose)

### Examples

```julia-repl
julia> A = [2..4 -2..1;-1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> b = [-2..2, -2..2]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-2, 2]

julia> solve(A, b, NonLinearOettliPrager(0.1))
Paving:
- tolerance ϵ = 0.1
- inner approx. of length 1195
- boundary approx. of length 823
```
"""
struct NonLinearOettliPrager <: AbstractIterativeSolver
    tol::Float64
end
NonLinearOettliPrager() = NonLinearOettliPrager(0.01)


text:
const DEFAULT_RUG_SPACE = 0.05

kde(x; kwargs...) = ArviZ.arviz.stats.density_utils.kde(x; kwargs...)

function kde(x, y; kwargs...)
    density, xmin, xmax, ymin, ymax = ArviZ.arviz.stats.density_utils._fast_kde_2d(
        x, y; kwargs...
    )
    nx, ny = size(density)
    return range(xmin, xmax; length=nx), range(ymin, ymax; length=ny), density
end

"""
    kdeplot(x[, y]; kwargs...)
    kdeplot!(x[, y]; kwargs...)

Plot a 1D or 2D kernel density estimate (KDE) taking into account boundary conditions.

# Positional Arguments
- `x`: Values to plot.
- `y`: Values to plot. If present, a 2D KDE will be estimated. See [`kde2dplot`](@ref).

# Keyword Arguments
- `bw=:default`: If numeric, indicates the bandwidth and must be positive. If `Symbol`,
    indicates the method to estimate the bandwidth and must be one of `:scott`,
    `:silverman`, `:isj` or `:experimental` when `circular=false` and `:taylor` (for now)
    when `circular=true`. `:default` means `:experimental` when variable is not circular
    and `:taylor` when it is.
- `circular=false`: Select input type in `(:radians, :degrees)` for circular KDE plot. If
    `true`, default input type is `:radians`. Pass `projection=:polar` to also plot on a
    polar axis; however, this will not plot correctly when `circular=:degrees`.
- `cumulative=false`: If `true`, plot the estimated cumulative distribution function
- `showrug=false`: If `true`, adds a rugplot
- `rugspace=$DEFAULT_RUG_SPACE`: control the `y` position of the rugplot. The larger this
    number the lower the rugplot
- `kwargs`: Additional attributes understood by Plots.jl. In particular, see keyword
    arguments for [`kde2dplot`](@ref) and [`rugplot`](@ref).
"""
RecipesBase.@shorthands kdeplot

RecipesBase.@recipe function f(
    ::Type{Val{:kdeplot}},
    x,
    y,
    z;
    bw=:default,
    circular=false,
    cumulative=false,
    showrug=false,
    rugspace=DEFAULT_RUG_SPACE,
)
    if x === eachindex(y) # x was (probably) auto-generated from y
        if bw === :default
            bw = circular === false ? :experimental : :taylor
        end
        xnew, ynew = kde(y; bw=bw, circular=(circular !== false))
        if cumulative
            ynew = cumsum(ynew)
            ynew ./= ynew[end]
        end
        if showrug
            ymax = maximum(ynew)
            if haskey(plotattributes, :rugsize)
                rugspace *= plotattributes[:rugsize] / DEFAULT_RUG_SIZE
            end
            RecipesBase.@series begin
                primary := false
                seriestype := :rugplot
                rugposition := -rugspace * ymax
                y := y
                ()
            end
        else
            # add dummy point to widen density axis to include density of 0
            RecipesBase.@series begin
                seriestype := :path
                primary := false
                seriesalpha := 0
                label := ""
                xpt, ypt = isrotated(plotattributes) ? (0, xnew[1]) : (xnew[1], 0)
                x := [xpt]
                y := [ypt]
                ()
            end
        end
        seriestype := :path
        if isrotated(plotattributes)
            xnew, ynew = ynew, xnew
        end
        x := xnew
        y := ynew
    else
        seriestype := :kde2dplot
        x := x
        y := y
    end
    return ()
end

"""
    kde2dplot(x, y; kwargs...)
    kde2dplot!(x, y; kwargs...)

Plot a 2D kernel density estimate (KDE) taking into account boundary conditions.

# Positional Arguments
- `x`: `x` coordinates of the values
- `y`: `y` coordinates of the values

# Keyword Arguments
- `contour=true`: If `true`, plot the 2D KDE using contours, otherwise plot a heatmap.
- `kwargs`: Additional attributes understood by Plots.jl.
"""
RecipesBase.@shorthands kde2dplot

RecipesBase.@recipe function f(::Type{Val{:kde2dplot}}, x, y, z; contour=true)
    xnew, ynew, density = kde(x, y)
    seriestype := contour ? :contour : :heatmap
    if isrotated(plotattributes)
        xnew, ynew = ynew, xnew
    end
    x := xnew
    y := ynew
    z := RecipesPipeline.Surface(density')
    colorbar --> false
    grid --> false
    return ()
end


text:
"""
    AbstractLatent
"""
abstract type AbstractLatent end

include("1D.jl")
include("ND.jl")


text:
using Gadfly
using TSne

"""
Normalize `A` columns, so that the mean and standard deviation
of each column are 0 and 1, resp.
"""
function rescale(A, dim::Integer=1)
    res = A .- mean(A, dim)
    res ./= map!(x -> x > 0.0 ? x : 1.0, std(A, dim))
    res
end

if length(ARGS)==0
    println("usage:\n\tjulia demo.jl iris\n\tjulia demo.jl mnist")
    exit(0)
end

if ARGS[1] == "iris"
    using RDatasets
    println("Using Iris dataset.")
    iris = dataset("datasets","iris")
    X = convert(Matrix{Float64}, iris[:, 1:4])
    labels = iris[:, 5]
    plotname = "iris"
    initial_dims = -1
    iterations = 1500
    perplexity = 15
elseif ARGS[1] == "mnist"
    using MLDatasets
    println("Using MNIST dataset.")
    X, labels = MNIST.traindata(Float64);
    npts = min(2500, size(X, 2), size(labels))
    labels = labels[1:npts]
    X = rescale(X[:, 1:npts]')
    plotname = "mnist"
    initial_dims = 50
    iterations = 1000
    perplexity = 20
else
    error("Unknown dataset \"", ARGS[1], "\"")
end

println("X dimensions are: ", size(X))
Y = tsne(X, 2, initial_dims, iterations, perplexity)
println("Y dimensions are: ", size(Y))

writecsv(plotname*"_tsne_out.csv", Y)
open("labels.txt", "w") do io
    writedlm(io, labels)
end

theplot = plot(x=Y[:,1], y=Y[:,2], color=labels)
draw(PDF(plotname*".pdf", 4inch, 3inch), theplot)
#draw(SVG(plotname*".svg", 4inch, 3inch), theplot)


text:
__precompile__()

module StarMatch

using Logging
using LinearAlgebra
using StaticArrays
using ProgressMeter
using PushVectors

export generatespd, solve
export CoordinateVector, Camera, CatalogStar

const N_NEAREST = 5
const VOTE_THRESHOLD = 3
const VECTOR_TOLERANCE = 4

struct CoordinateVector{T<:Real} <: AbstractVector{T}
    p::Vector{SVector{2,T}}
end
CoordinateVector(::Type{T}, dims::Int) where {T} = CoordinateVector{T}(Vector{SVector{2,T}}(undef, dims));
Base.size(A::CoordinateVector) = size(A.p)
Base.IndexStyle(::Type{<:CoordinateVector}) = IndexLinear()
Base.getindex(A::CoordinateVector, i::Int) = A.p[i]
Base.setindex!(A::CoordinateVector{T}, v, i::Int) where {T} = (A.p[i] = v)
Base.similar(A::CoordinateVector{T}) where {T} = CoordinateVector(T, size(A))
Base.similar(A::CoordinateVector, ::Type{T}, dims::Dims) where {T} = CoordinateVector(T, dims)
Base.eltype(A::CoordinateVector{T}) where {T} = T
Base.deleteat!(A::CoordinateVector{T}, inds) where {T} = deleteat!(A.p, inds)

struct CatalogStar
  id::Int
  ra::Float64
  dec::Float64
  mag::Float64
end

abstract type ImageStar end

struct UnknownStar <: ImageStar
    xy::SVector{2, Float64}
    vec::SVector{2, Float64}
    uvec::SVector{2, Float64}
    d::Float64
end

struct KnownStar <: ImageStar
    catalogidx::Int
    xy::SVector{2, Float64}
    vec::SVector{2, Float64}
    uvec::SVector{2, Float64}
    d::Float64
end
Base.zero(::Type{KnownStar}) = KnownStar(0, SVector(0.0,0.0), SVector(0.0,0.0), SVector(0.0,0.0), 0.0)

struct SPDEntry
    starOidx::Int  # catalog idx
    SPidx::Int  # catalog idx
    d::Float64
    path::CoordinateVector{Float64}
    pathidxs::Vector{Int}  # catalog idx
end

struct Camera
    img_width::Int
    img_height::Int
    img_center::SVector{2,Float64}
    pixelsize::Float64
    focallength::Float64
    fov::Float64
end
function Camera(w, h, pixelsize, focallength)
    imdim = max(h, w)
    fov = atand(imdim*pixelsize/2/focallength)*2
    center = SVector(w/2, h/2)
    Camera(w, h, center, pixelsize, focallength, fov)
end

"""
Transform RA/DEC to ECI unit vector.

RA and DEC in degrees.

See Vallado 4-1.
"""
function radec2eci(ra, dec)
    sinδ, cosδ = sincos(deg2rad(dec))
    sinα, cosα = sincos(deg2rad(ra))
    return SVector(cosα*cosδ, sinα*cosδ, sinδ)
end

"""
Transform from RA/DEC to image plane XY.

See Vallado chapter 4 and the SPD code from Samirbhai.
"""
function cameraattitude(Veci)
    tmp = 1/sqrt(Veci[1]^2 + Veci[2]^2)
    c1x = Veci[2]*tmp
    c1y = -Veci[1]*tmp
    c1z = 0
    c1 = SVector(c1x, c1y, c1z)

    c2 = cross(Veci, c1)

    C = @SMatrix [c1[1] c1[2] c1[3];
                  c2[1] c2[2] c2[3];
                  Veci[1] Veci[2] Veci[3]]

    return C
end

"""
    camera2image(Vcam, camera)

Convert vector `Vcam` in camera frame to image plane coordinates.
"""
function camera2image(Vcam, camera::Camera)
    f = camera.focallength
    ρ = camera.pixelsize
    w = camera.img_width
    h = camera.img_height

    x = f/ρ*Vcam[1]/Vcam[3] + w/2
    y = f/ρ*Vcam[2]/Vcam[3] + h/2

    return SVector(x, y)
end

"""
    buildpath(neighbors, V1, V1u)

Construct a rotation invariant vector path along `neighbors` with vector `V1` and unit
vector `V1u`.
"""
function buildpath(neighbors::AbstractVector{<:ImageStar}, V1, V1u)
    path = CoordinateVector(Float64, length(neighbors)-1)
    Vbsum = zero(V1)
    @inbounds for j in 1:length(neighbors)-1
        k = j + 1

        Vj = neighbors[j].vec
        Vk = neighbors[k].vec
        Vku = neighbors[k].uvec

        Vbsum += (Vk - Vj)

        path[j] = SVector(dot(Vbsum, V1u), dot(Vbsum, Vku))
    end
    return path
end

"""
    vote(candidateindices, neighborpaths, spd, vectortolerance)

Give a vote if norm of the difference for `candidateindices` of `neighborpaths` and `spd`
are within `vectortolerance`.
"""
function vote(
    candidateindices::AbstractVector{Tuple{Int,Int}},
    neighborpaths::AbstractVector{CoordinateVector{T}},
    spd::AbstractVector{SPDEntry},
    vectortolerance=VECTOR_TOLERANCE) where T

    # TODO: Reuse `votes` when there are multiple images
    votes = zeros(UInt32, length(candidateindices))
    for (c, (i, j)) in enumerate(candidateindices)
        votecount = 0
        spdpath = spd[j].path
        neighborpath = neighborpaths[i]
        for sp in spdpath
            for np in neighborpath
                if abs(sp[1] - np[1]) < vectortolerance && abs(sp[2] - np[2]) < vectortolerance
                    votecount += 1
                end
            end
        end
        votes[c] += votecount
    end
    return votes
end

"""
    closest(X, a)

Return index of `X` that is closest distance to `a`.
"""
function closest(X::CoordinateVector{T}, a::SVector{2,T}) where T <: Real
    mindist = typemax(T)
    mindistidx = 0
    for (i, x) in enumerate(X)
        dist = norm(a - x)
        if dist < mindist
            mindist = dist
            mindistidx = i
        end
    end
    return mindistidx
end

"""
    match!(matches, winner, neighbors, neighborpath, vectortolerance)

Match `neighbors` to path indexes of `winner` and add them as `KnownStar`s to `matches`.
"""
function match!(matches, winner::SPDEntry, neighbors::AbstractVector{UnknownStar},
    neighborpath::CoordinateVector, vectortolerance=VECTOR_TOLERANCE) where T <: Real

    # The first neighbor _must_ be correct, otherwise we wouldn't have the correct solution
    push!(matches, KnownStar(winner.SPidx, neighbors[1].xy, neighbors[1].vec,
                             neighbors[1].uvec, neighbors[1].d))

    for i in eachindex(neighborpath)
        np = neighborpath[i]
        idx = closest(winner.path, np)
        @inbounds wp = winner.path[idx]

        if abs(np[1] - wp[1]) < vectortolerance && abs(np[2] - wp[2]) < vectortolerance
            push!(matches, KnownStar(winner.pathidxs[idx+1], neighbors[i+1].xy,
                                     neighbors[i+1].vec, neighbors[i+1].uvec,
                                     neighbors[i+1].d))
        end
    end

    return matches
end

"""
    generatespd(camera, catalog)

Generate a Star Pattern Database for `camera` and a `catalog` of `CatalogStar`s.
"""
function generatespd(camera::Camera, catalog::Vector{CatalogStar})
    # Common values
    halffov = camera.fov/2
    coshalffov = cosd(halffov)
    imagecenter = camera.img_center

    # Calculating star position in ECI is surprisingly expensive, so we do it once
    catalogeci = [radec2eci(s.ra, s.dec) for s in catalog]

    spd = PushVector{SPDEntry}(length(catalog)*N_NEAREST)
    neighbors = PushVector{KnownStar}()

    @showprogress 5 "Building SPD..." for (oi, starO) in enumerate(catalog)
        VeciO = catalogeci[oi]
        CO = cameraattitude(VeciO)
        xyO = imagecenter  # == camera2image(CO*VeciO, camera)

        # Identify neighboring stars in FOV
        for (si, s) in enumerate(catalog)
            @inbounds Veci = catalogeci[si]

            # Is star `s` within fov/2 of `starO`?
            if dot(Veci, VeciO) > coshalffov && Veci != VeciO
                # Star coordinates in camera frame
                Vcam = CO*Veci

                # Star coordinates on image plane
                xy = camera2image(Vcam, camera)

                # Vector from `starO` (center of image) to star `s`
                svec = xy - xyO

                # Distance from `starO`
                d = norm(svec)

                # Unit vector from `starO` to star `s`
                usvec = svec/d

                push!(neighbors, KnownStar(si, xy, svec, usvec, d))
            end
        end

        sort!(neighbors, by = x -> x.d)

        #==
        NOTE: Samirbhai removes the first N_NEAREST-1 entries of his boundary vector (path).

        They are included here because we tend to have few stars in our FOV and want to use
        as many of them as possible.

        By enforcing there be N_NEAREST+2 neighbors, we ensure the path is at least 2 long
        ==#
        length(neighbors) < (N_NEAREST+2) && continue

        @inbounds for i in 1:N_NEAREST
            V1 = neighbors[i].vec
            V1u = neighbors[i].uvec

            sn = @view neighbors[i:end]
            path = buildpath(sn, V1, V1u)
            pathidxs = [n.catalogidx for n in sn]

            push!(spd, SPDEntry(oi, neighbors[i].catalogidx, neighbors[i].d, path, pathidxs))
        end
        empty!(neighbors)
    end
    PushVectors.finish!(spd)
    return spd
end

"""
    solve(camera, imagestars, spd; distancetolerance=3, vectortolerance=VECTOR_TOLERANCE)

Return catalog matches from `spd` for `imagestars` extracted from an image taken with
`camera`. `distancetolerance` is the number of pixels within which the distance from the
nearest neighbor stars must be to the center star compared to the distances stored in the
`spd`. `vectortolerance` is the tolerance used in `vote` and `match` to check the `norm` of
the difference of path vectors.
"""
function solve(
    camera::Camera,
    imagestars::CoordinateVector{T},
    spd::AbstractVector{SPDEntry},
    distancetolerance=3,
    vectortolerance=VECTOR_TOLERANCE) where T <: Real

    if length(imagestars) <= N_NEAREST
        return nothing
    end

    imagecenter = camera.img_center

    # Determine star closest to center. This will be `starO`
    starOidx = closest(imagestars, imagecenter)
    starO = imagestars[starOidx]

    # Calculate vectors from each star to `starO`
    neighbors = Vector{UnknownStar}(undef, length(imagestars)-1)
    ni = 1
    for star in imagestars
        if star != starO
            svec = star - starO
            d = norm(svec)
            usvec = svec/d

            @inbounds neighbors[ni] = UnknownStar(star, svec, usvec, d)
            ni += 1
        # else starOid = i  # just fyi
        end
    end
    sort!(neighbors, by = x -> x.d)

    # TODO: Reuse `candidateindices` when we handle multiple images
    # Check for distance matches with SPD, approx that 1/20 are candidates
    candidateindices = PushVector{Tuple{Int,Int}}(floor(Int, length(spd)/20))
    for j in eachindex(spd)
        spddist = spd[j].d
        @inbounds for i in 1:N_NEAREST
            if norm(spddist - neighbors[i].d) < distancetolerance
                push!(candidateindices, (i, j))
            end
        end
    end
    PushVectors.finish!(candidateindices)

    neighborpaths = Vector{CoordinateVector{Float64}}(undef, N_NEAREST)
    @inbounds for i in 1:N_NEAREST
        V1 = neighbors[i].vec
        V1u = neighbors[i].uvec

        sn = @view neighbors[i:end]
        neighborpaths[i] = buildpath(sn, V1, V1u)
    end

    votecounts = vote(candidateindices, neighborpaths, spd, vectortolerance)

    winnerid = argmax(votecounts)
    winnervotecount = votecounts[winnerid]

    if winnervotecount < VOTE_THRESHOLD
        @warn "Winner vote count is below threshold! ($winnervotecount/$VOTE_THRESHOLD)"
        return nothing
    elseif count(votecounts .== winnervotecount) > 1
        @warn "Vote resulted in a tie!"
        return nothing
    end

    winningindices = candidateindices[winnerid]
    winner = spd[winningindices[2]]

    matches = PushVector{KnownStar}(length(neighbors))
    push!(matches, KnownStar(winner.starOidx, SVector{2,Float64}(starO), SVector(0.0, 0.0),
                             SVector(0.0, 0.0), 0.0))

    sn = @view neighbors[winningindices[1]:end]
    match!(matches, winner, sn, neighborpaths[winningindices[1]], vectortolerance)

    return matches
end

end # module


text:
# Some useful constants
const sixoverpi = 6.0 / π
const fourpi = 4.0 * π

"""
    ParticleSystem{T <: AbstractVector}

A handy type to manipulate the positions of a set of given particles.

# Fields
- `x::AbstractVector{T}`: stores the positions of the particles in the x
    directions, in rectangular coordinates
- `y::AbstractVector{T}`: stores the positions of the particles in the x
    directions, in rectangular coordinates
- `z::AbstractVector{T}`: stores the positions of the particles in the x
    directions, in rectangular coordinates
"""

mutable struct ParticleSystem{T<:AbstractArray}
    positions::T
    forces::T
end

function ParticleSystem(N::V, M::V, ::Type{T}) where {T<:AbstractFloat,V<:Integer}
    pos = zeros(T, N, M)
    forces = zeros(T, N, M)
    return ParticleSystem{AbstractArray{T}}(pos, forces)
end

"""
"""
mutable struct Parameters{U<:AbstractFloat}
    ϕ::U
    kT::U
    N::Integer
    τ::U
    seed::Integer
end

# Just define packing fraction, temperature and time step
function Parameters(
    ϕ::AbstractFloat,
    kT::AbstractFloat,
    τ::AbstractFloat,
    ::Type{U},
) where {U}
    return Parameters{U}(U(ϕ), U(kT), 512, U(τ), 0)
end

# Keep default seed, define total number of particles
function Parameters(
    ϕ::AbstractFloat,
    kT::AbstractFloat,
    N::Integer,
    τ::AbstractFloat,
    ::Type{U},
) where {U}
    return Parameters{U}(U(ϕ), U(kT), N, U(τ), 0)
end

# Define a different seed
function Parameters(
    ϕ::AbstractFloat,
    kT::AbstractFloat,
    N::Integer,
    τ::AbstractFloat,
    s::Integer,
    ::Type{U},
) where {U}
    return Parameters{U}(U(ϕ), U(kT), N, U(τ), s)
end

"""
    SimulationSystem{U<:AbstractFloat} <: PhysicalSystem

A type that stores the most important parameters to perform a simulation.
These values are computed through the `Parameters` type.

# Fields
- `ρ::U`: the packing fraction of the system
"""
mutable struct SimulationSystem{U<:AbstractFloat}
    params::Parameters
    ρ::U
    boxl::U
    rc::U
    system::ParticleSystem
    energy::U
    dims::Integer
end

function SimulationSystem(params::Parameters, dims::Integer, ::Type{T}) where {T}
    # Assign density for a 3D system
    ρ = params.ϕ * sixoverpi
    boxl = ∛(params.N / ρ)
    rc = boxl * 0.5
    # Build a system of particles with the total number of particles
    psys = ParticleSystem(params.N, dims, T)
    return SimulationSystem{T}(params, T(ρ), T(boxl), T(rc), psys, T(0.0), dims)
end

abstract type Structure end
mutable struct PairDistributionFunction <: Structure
    gofr::AbstractArray
    nm::Integer
    dr::AbstractFloat
    naverage::AbstractFloat
    norm_const::AbstractFloat
end

function PairDistributionFunction(s::SimulationSystem, nm::Integer)
    normalizing_constant = fourpi * s.ρ
    gofr = zeros(typeof(s.ρ), (nm, 2))
    dr = s.rc / nm
    return PairDistributionFunction(gofr, nm, dr, 0.0, normalizing_constant)
end

mutable struct StructureFactor{T<:AbstractFloat} <: Structure
    sq::AbstractArray
    dq::T
    nm::Integer
    naverage::T
    norm_const::T
end

function StructureFactor(s::SimulationSystem, rdf::PairDistributionFunction)
    compute_type = typeof(s.ρ)
    normalizing_constant = fourpi * s.ρ
    sq = zeros(compute_type, (rdf.nm, 2))
    dq = 0.1 * π / s.rc
    return StructureFactor{compute_type}(sq, dq, rdf.nm, rdf.naverage, normalizing_constant)
end

abstract type Thermodynamics end
mutable struct ZFactor{T<:AbstractFloat} <: Thermodynamics
    zval::T
    naverage::T
end

function compute_z(zfac::ZFactor, N::Real)
    return 1.0 - (zfac.zval / (3.0 * zfac.naverage * N))
end

abstract type Dynamics end
mutable struct MeanSquaredDisplacement{V<:Integer} <: Dynamics
    displacement::AbstractArray
    mt::V
    wt::AbstractArray
    naverage::V
    interval::V
end

function MeanSquaredDisplacement(s::SimulationSystem, mt::Integer, interval::Integer)
    displacement = zeros(mt, s.params.N, s.dims)
    wt = zeros(mt, 2)
    return MeanSquaredDisplacement(displacement, mt, wt, 0, interval)
end

mutable struct SelfScatteringFunction{V<:Integer} <: Dynamics
    mt::V
    dft::AbstractArray
    naverage::V
    interval::V
    κ::AbstractFloat
end

function SelfScatteringFunction(
    s::SimulationSystem,
    mt::Integer,
    interval::Integer,
    κ::AbstractFloat,
)
    dft = zeros(mt, 2)
    # Use the standard wave number value of 6.6
    return SelfScatteringFunction(mt, dft, 0, interval, κ)
end


text:
module SinkTest

using HPAT

@acc hpat function sink(file3)
    y = randn(10)
    DataSink(y,HDF5,"/y", file3)
end

end

SinkTest.sink("y.hdf5")
using MPI
if MPI.Comm_rank(MPI.COMM_WORLD)==0 rm("y.hdf5") end


text:
export RockPaperScissorsEnv

"""
    RockPaperScissorsEnv()

[Rock Paper Scissors](https://en.wikipedia.org/wiki/Rock_paper_scissors) is a
simultaneous, zero sum game.
"""
Base.@kwdef mutable struct RockPaperScissorsEnv <: AbstractEnv
    reward::Tuple{Int,Int} = (0, 0)
    is_done::Bool = false
end

RLBase.players(::RockPaperScissorsEnv) = (1, 2)

"""
Note that although this is a two player game. The current player is always a
dummy simultaneous player.
"""
RLBase.current_player(::RockPaperScissorsEnv) = SIMULTANEOUS_PLAYER

# Defining the `action_space` of each independent player can help to transform
# this SIMULTANEOUS environment into a SEQUENTIAL environment with
# [`simultaneous2sequential`](@ref).
RLBase.action_space(::RockPaperScissorsEnv, ::Int) = Space(('💎', '📃', '✂'))

RLBase.action_space(::RockPaperScissorsEnv, ::SimultaneousPlayer) =
    Space(Tuple((i, j) for i in ('💎', '📃', '✂') for j in ('💎', '📃', '✂')))

RLBase.action_space(env::RockPaperScissorsEnv) = action_space(env, SIMULTANEOUS_PLAYER)

RLBase.legal_action_space(env::RockPaperScissorsEnv, p) =
    is_terminated(env) ? () : action_space(env, p)

"Since it's a one-shot game, the state space doesn't have much meaning."
RLBase.state_space(::RockPaperScissorsEnv, ::Observation, p) = Space(OneTo(1))

"""
For multi-agent environments, we usually implement the most detailed one.
"""
RLBase.state(::RockPaperScissorsEnv, ::Observation, p) = 1

RLBase.reward(env::RockPaperScissorsEnv) = env.is_done ? env.reward : (0, 0)
RLBase.reward(env::RockPaperScissorsEnv, p::Int) = reward(env)[p]

RLBase.is_terminated(env::RockPaperScissorsEnv) = env.is_done
RLBase.reset!(env::RockPaperScissorsEnv) = env.is_done = false

function (env::RockPaperScissorsEnv)((x, y))
    if x == y
        env.reward = (0, 0)
    elseif x == '💎' && y == '✂' || x == '✂' && y == '📃' || x == '📃' && y == '💎'
        env.reward = (1, -1)
    else
        env.reward = (-1, 1)
    end
    env.is_done = true
end

RLBase.NumAgentStyle(::RockPaperScissorsEnv) = MultiAgent(2)
RLBase.DynamicStyle(::RockPaperScissorsEnv) = SIMULTANEOUS
RLBase.ActionStyle(::RockPaperScissorsEnv) = MINIMAL_ACTION_SET
RLBase.InformationStyle(::RockPaperScissorsEnv) = IMPERFECT_INFORMATION
RLBase.StateStyle(::RockPaperScissorsEnv) = Observation{Int}()
RLBase.RewardStyle(::RockPaperScissorsEnv) = TERMINAL_REWARD
RLBase.UtilityStyle(::RockPaperScissorsEnv) = ZERO_SUM
RLBase.ChanceStyle(::RockPaperScissorsEnv) = DETERMINISTIC


text:
function get_body_oriented_frame(obj::JPLCelestialBody)
    return jcall(obj, "getBodyOrientedFrame", Frame, ())
end

function get_gm(obj::JPLCelestialBody)
    return jcall(obj, "getGM", jdouble, ())
end

function get_inertially_oriented_frame(obj::JPLCelestialBody)
    return jcall(obj, "getInertiallyOrientedFrame", Frame, ())
end

function get_name(obj::JPLCelestialBody)
    return jcall(obj, "getName", JString, ())
end

function get_pv_coordinates(obj::JPLCelestialBody, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function get_pv_coordinates(obj::JPLCelestialBody, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build MPCBuilder
name = "ARB"
version = v"2.16.0"
sources = [
    "https://github.com/fredrik-johansson/arb/archive/$version.tar.gz" =>
    "77464be4d34a511bb004457f862fec857ff934b0ed58d56d6f52d76ebadd4daf",
]

# Bash recipe for building across all platforms

script = """cd arb-$version
"""*raw"""
UNAME=`uname`
if [[ ${UNAME} == MSYS_NT-6.3 ]]; then
    ln -sf $prefix/bin/libflint-13.dll $prefix/lib/libflint.dll.a
fi
./configure --prefix=$prefix --enable-shared --disable-static --with-gmp=$prefix --with-mpfr=$prefix --with-flint=$prefix
make -j
make install
if [[ ${UNAME} == MSYS_NT-6.3 ]]; then
    rm -f $prefix/lib/libflint.dll.a
fi
"""
# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Windows(:i686, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.Windows(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:i686, :glibc, :blank_abi)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libarb", :libarb)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl"
    "https://github.com/JuliaMath/MPFRBuilder/releases/download/v4.0.1/build.jl"
    "https://github.com/isuruf/FLINTBuilder/releases/download/v2.5.3%2Bdev0.ea9f102a0b/build_FLINT.v2.5.3+dev0.ea9f102a0b.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: isuruf/ARBBuilder

max_issues_repo_name: isuruf/ARBBuilder

max_forks_repo_name: isuruf/ARBBuilder


text:
module Types

export LLVMBool, LLVMPtr, ContextPtr, TypePtr, GlobalPtr, GlobalValuePtr, ValuePtr, ModulePtr,
       ConstPtr, UserPtr, GlobalAliasPtr, FunctionPtr, RawOStreamPtr, MemoryBufferPtr,
       SMDiagnosticPtr, BasicBlockPtr, BuilderPtr, MDNodePtr, InstructionPtr, ParamPtr,
       ExeEnginePtr, JITCompilerOpts, isnull

typealias LLVMBool  Cint
typealias CodeModel Cint

abstract LLVMPtr

isnull(h::LLVMPtr) = is(h.ptr, zero(Ptr{Void}))

Base.convert(::Type{Ptr{Void}}, ptr::LLVMPtr) = ptr.ptr

Base.(:(==))(p1::LLVMPtr, p2::LLVMPtr) = p1.ptr === p2.ptr
Base.isequal(p1::LLVMPtr, p2::LLVMPtr) = p1.ptr == p2.ptr

Base.hash(p::LLVMPtr) = hash(p.ptr)

Base.show(io::IO, ptr::LLVMPtr) = begin
    addr = "0x$(hex(unsigned(convert(Ptr{Void}, ptr)), WORD_SIZE>>2))"
    print(io, "LLVM.$(typeof(ptr))(@$addr)")
end

immutable SMDiagnosticPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable RawOStreamPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable MemoryBufferPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable BuilderPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable DataLayoutPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable TargetPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable TargetOptionsPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable TargetLibInfoPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable TargetMachinePtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable TargetLoweringPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable ExeEnginePtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable MemManagerPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable JITCompilerOpts
    optlevel::Cuint
    codemodel::CodeModel
    no_frame_ptr_elim::LLVMBool
    enable_fast_isel::LLVMBool
    memory_manager::MemManagerPtr

    JITCompilerOpts(optlevel,
                    codemodel,
                    no_frame_ptr_elim::Bool,
                    enable_fast_isel::Bool) = begin
        @assert optlevel >= 0
        memory_manager = MemManagerPtr(C_NULL)
        return new(optlevel,
                   codemodel,
                   no_frame_ptr_elim,
                   enable_fast_isel,
                   memory_manager)
    end
end

immutable MCJITOptsPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable PassManagerPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable PassBuilderPtr <: LLVMPtr
    ptr::Ptr{Void}
end

#---------------------------------#

immutable ContextPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable ModulePtr <: LLVMPtr
    ptr::Ptr{Void}
end

# Pointer Hierarchy

immutable ValuePtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable UserPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::UserPtr) = ValuePtr(ptr.ptr)

immutable ConstPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{UserPtr}, ptr::ConstPtr) = UserPtr(ptr.ptr)
Base.convert(::Type{ValuePtr}, ptr::ConstPtr) = ValuePtr(ptr.ptr)

abstract GlobalPtr <: LLVMPtr

immutable GlobalValuePtr <: GlobalPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::GlobalValuePtr) = ValuePtr(ptr.ptr)
Base.convert(::Type{UserPtr},  ptr::GlobalValuePtr) = UserPtr(ptr.ptr)
Base.convert(::Type{ConstPtr}, ptr::GlobalValuePtr) = ConstPtr(ptr.ptr)

immutable GlobalAliasPtr <: GlobalPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::GlobalAliasPtr) = ValuePtr(ptr.ptr)
Base.convert(::Type{UserPtr},  ptr::GlobalAliasPtr) = UserPtr(ptr.ptr)
Base.convert(::Type{ConstPtr}, ptr::GlobalAliasPtr) = ConstPtr(ptr.ptr)
Base.convert(::Type{GlobalValuePtr}, ptr::GlobalAliasPtr) = GlobalValuePtr(ptr.ptr)

immutable FunctionPtr <: GlobalPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::FunctionPtr) = ValuePtr(ptr.ptr)
Base.convert(::Type{UserPtr},  ptr::FunctionPtr) = UserPtr(ptr.ptr)
Base.convert(::Type{ConstPtr}, ptr::FunctionPtr) = ConstPtr(ptr.ptr)
Base.convert(::Type{GlobalValuePtr}, ptr::FunctionPtr) = GlobalValuePtr(ptr.ptr)

immutable BasicBlockPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::BasicBlockPtr) = ValuePtr(ptr.ptr)

immutable ParamPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::ParamPtr) = ValuePtr(ptr.ptr)

immutable InstructionPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::InstructionPtr) = ValuePtr(ptr.ptr)
Base.convert(::Type{UserPtr},  ptr::InstructionPtr)  = UserPtr(ptr.ptr)

#TODO: temp hack, we should have a constant instruction ptr
Base.convert(::Type{InstructionPtr}, ptr::ConstPtr) = InstructionPtr(ptr.ptr)

immutable BinaryOpPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::BinaryOpPtr) = ValuePtr(ptr.ptr)
Base.convert(::Type{UserPtr},  ptr::BinaryOpPtr) = UserPtr(ptr.ptr)
Base.convert(::Type{InstructionPtr}, ptr::BinaryOpPtr) = InstructionPtr(ptr.ptr)

immutable MDNodePtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::MDNodePtr) = ValuePtr(ptr.ptr)

immutable MDStringPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::MDStringPtr) = ValuePtr(ptr.ptr)

immutable NamedMetadataPtr <: LLVMPtr
    ptr::Ptr{Void}
end

immutable InlineAsmPtr <: LLVMPtr
    ptr::Ptr{Void}
end

Base.convert(::Type{ValuePtr}, ptr::InlineAsmPtr) = ValuePtr(ptr.ptr)

immutable TypePtr <: LLVMPtr
    ptr::Ptr{Void}
end

end


text:
struct Effect
    type::Int # What effect
    value::Int # How strong the effect is
end

struct Recipe
    required::Vector{Int}
    quantity::Vector{Int}
end

struct Item
    name::String
    id::Int
    edible::Bool
    fishing::Int
    exploit::Int
    effects::Vector{Effect}
    recipes::Vector{Recipe}
end

struct AllItems
    v::Vector{Item}
end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

"""
    OpticalSystem{T<:Real}

Abstract type for any optical system, must parameterized by the datatype of entities within the system `T`.
"""
abstract type OpticalSystem{T<:Real} end
export OpticalSystem

"""
    CSGOpticalSystem{T,D<:Real,S<:Surface{T},L<:LensAssembly{T}} <: OpticalSystem{T}

An optical system containing a lens assembly with all optical elements and a detector surface with associated image. The system can be at a specified temperature and pressure.

There are two number types in the type signature. The `T` type parameter is the numeric type for geometry in the optical system, the `D` type parameter is the numeric type of the pixels in the detector image.
This way you can have `Float64` geometry, where high precision is essential, but the pixels in the detector can be `Float32` since precision is much less critical for image data.

The detector can be any [`Surface`](@ref) which implements [`uv`](@ref), [`uvtopix`](@ref) and [`onsurface`](@ref), typically this is one of [`Rectangle`](@ref), [`Ellipse`](@ref) or [`SphericalCap`](@ref).

```julia
CSGOpticalSystem(assembly::LensAssembly, detector::Surface, detectorpixelsx = 1000, detectorpixelsy = 1000, ::Type{D} = Float32; temperature = OpticSim.GlassCat.TEMP_REF, pressure = OpticSim.GlassCat.PRESSURE_REF)
```
"""
struct CSGOpticalSystem{T,D<:Number,S<:Surface{T},L<:LensAssembly{T}} <: OpticalSystem{T}
    assembly::L
    detector::S
    detectorimage::HierarchicalImage{D}
    temperature::T
    pressure::T
    function CSGOpticalSystem(assembly::L, detector::S, detectorpixelsx::Int = 1000, detectorpixelsy::Int = 1000, ::Type{D} = Float32; temperature::Union{T,Unitful.Temperature} = T(OpticSim.GlassCat.TEMP_REF), pressure::T = T(OpticSim.GlassCat.PRESSURE_REF)) where {T<:Real,S<:Surface{T},L<:LensAssembly{T},D<:Number}
        @assert hasmethod(uv, (S, SVector{3,T})) "Detector must implement uv()"
        @assert hasmethod(uvtopix, (S, SVector{2,T}, Tuple{Int,Int})) "Detector must implement uvtopix()"
        @assert hasmethod(onsurface, (S, SVector{3,T})) "Detector must implement onsurface()"
        opticalinterface = interface(detector)
        @assert insidematerialid(opticalinterface) == outsidematerialid(opticalinterface) "Detector must have same material either side"
        @assert interface(detector) !== NullInterface(T) "Detector can't have null interface"
        image = HierarchicalImage{D}(detectorpixelsy, detectorpixelsx)
        if temperature isa Unitful.Temperature
            temperature = Unitful.ustrip(T, Unitful.u"°C", temperature)
        end
        return new{T,D,S,L}(assembly, detector, image, temperature, T(pressure))
    end
end
export CSGOpticalSystem

Base.copy(a::CSGOpticalSystem) = CSGOpticalSystem(a.assembly, a.detector, size(a.detectorimage)..., temperature = a.temperature * Unitful.u"°C", pressure = a.pressure)

#added this show method because the type of CSGOpticalSystem is gigantic and printing it in the REPL can crash the system.
Base.show(io::IO, a::CSGOpticalSystem{T}) where {T} = print(io, "CSGOpticalSystem{$T}($(temperature(a)), $(pressure(a)), $(assembly(a)), $(detector(a)))")

"""
    assembly(system::OpticalSystem{T}) -> LensAssembly{T}

Get the [`LensAssembly`](@ref) of `system`.
"""
assembly(system::CSGOpticalSystem{T}) where {T<:Real} = system.assembly
detector(system::CSGOpticalSystem) = system.detector
"""
    detectorimage(system::OpticalSystem{T}) -> HierarchicalImage{D}

Get the detector image of `system`.
`D` is the datatype of the detector image and is not necessarily the same as the datatype of the system `T`.
"""
detectorimage(system::CSGOpticalSystem) = system.detectorimage
detectorsize(system::CSGOpticalSystem) = size(system.detectorimage)
"""
    temperature(system::OpticalSystem{T}) -> T

Get the temperature of `system` in °C.
"""
temperature(system::CSGOpticalSystem{T}) where {T<:Real} = system.temperature
"""
    pressure(system::OpticalSystem{T}) -> T

Get the pressure of `system` in Atm.
"""
pressure(system::CSGOpticalSystem{T}) where {T<:Real} = system.pressure
"""
    resetdetector!(system::OpticalSystem{T})

Reset the deterctor image of `system` to zero.
"""
resetdetector!(system::CSGOpticalSystem{T}) where {T<:Real} = reset!(system.detectorimage)
export temperature, pressure, detectorimage, resetdetector!, assembly

Base.Float32(a::T) where {T<:ForwardDiff.Dual} = Float32(ForwardDiff.value(a))

"""
    trace(system::OpticalSystem{T}, ray::OpticalRay{T}; trackrays = nothing, test = false)

Traces `system` with `ray`, if `test` is enabled then fresnel reflections are disabled and the power distribution will not be correct.
Returns either a [`LensTrace`](@ref) if the ray hits the detector or `nothing` otherwise.

`trackrays` can be passed an empty vector to accumulate the `LensTrace` objects at each intersection of `ray` with a surface in the system.
"""
function trace(system::CSGOpticalSystem{T,D}, r::OpticalRay{T,N}; trackrays::Union{Nothing,Vector{LensTrace{T,N}}} = nothing, test::Bool = false) where {T<:Real,N,D<:Number}
    if power(r) < POWER_THRESHOLD
        return nothing
    end

    result = trace(system.assembly, r, temperature(system), pressure(system), trackrays = trackrays, test = test)

    if result === nothing || result === nopower
        emptyintervalpool!(T)
        return nothing
    else #ray intersected lens assembly so continue to see if ray intersects detector
        intsct = surfaceintersection(detector(system), ray(result))
        if intsct === nothing # no intersection of final ray with detector
            emptyintervalpool!(T)
            return nothing
        end

        detintsct = closestintersection(intsct)
        if detintsct === nothing
            emptyintervalpool!(T)
            return nothing
        else
            # need to modify power and path length accordingly for the intersection with the detector
            surfintsct = point(detintsct)
            nml = normal(detintsct)
            opticalinterface = interface(detintsct)::FresnelInterface{T}
            λ = wavelength(r)

            # Optical path length is measured in mm
            # in this case the result ray is exact so no correction for RAY_OFFSET is needed
            geometricpathlength = norm(surfintsct - origin(ray(result)))
            opticalpathlength = geometricpathlength
            pow = power(result)

            m = outsidematerialid(opticalinterface)
            # compute updated power based on absorption coefficient of material using Beer's law
            # this will almost always not apply as the detector will be in air, but it's possible that the detector is not in air, in which case this is necessary
            if !isair(m)
                mat = glassforid(m)::OpticSim.GlassCat.Glass
                nᵢ = index(mat, λ, temperature = temperature(system), pressure = pressure(system))::T
                α = absorption(mat, λ, temperature = temperature(system), pressure = pressure(system))::T
                if α > zero(T)
                    internal_trans = exp(-α * geometricpathlength)
                    if rand() >= internal_trans
                        return nothing
                    end
                    pow = pow * internal_trans
                end
                opticalpathlength = nᵢ * geometricpathlength
            end

            temp = LensTrace{T,N}(OpticalRay(ray(ray(result)), pow, wavelength(result), opl = pathlength(result) + opticalpathlength, nhits = nhits(result) + 1, sourcenum = sourcenum(r), sourcepower = sourcepower(r)), detintsct)
            if trackrays !== nothing
                push!(trackrays, temp)
            end

            # incerement the detector image
            pixu, pixv = uvtopix(detector(system), uv(detintsct), size(system.detectorimage))
            system.detectorimage[pixv, pixu] += D(sourcepower(r)) # TODO will need to handle different detector image types a bit better than this

            # should be okay to assume intersection will not be a DisjointUnion for all the types of detectors we will be using.
            emptyintervalpool!(T)
            return temp
        end
    end
end

######################################################################################################################

"""
    AxisymmetricOpticalSystem{T,C<:CSGOpticalSystem{T}} <: OpticalSystem{T}

Optical system which has lens elements and an image detector, created from a `DataFrame` containing prescription data.

These tags are supported for columns: `:Radius`, `:SemiDiameter`, `:Surface`, `:Thickness`, `:Conic`, `:Aspherics`, `:Reflectance`, `:Material`, `:OptimizeRadius`, `:OptimizeThickness`, `:OptimizeConic`.
These tags are supported for entries in a `:Surface` column: `:Object`, `:Image`, `:Stop`
Assumes the `:Image` row will be the last row in the `DataFrame`.

In practice a [`CSGOpticalSystem`](@ref) is generated automatically and stored within this system.

```julia
AxisymmetricOpticalSystem{T}(prescription::DataFrame, detectorpixelsx = 1000, detectorpixelsy:: = 1000, ::Type{D} = Float32; temperature = OpticSim.GlassCat.TEMP_REF, pressure = OpticSim.GlassCat.PRESSURE_REF)
```
"""
struct AxisymmetricOpticalSystem{T,C<:CSGOpticalSystem{T}} <: OpticalSystem{T}
    system::C
    prescription::DataFrame
    semidiameter::T

    function AxisymmetricOpticalSystem{T}(prescription::DataFrame, detectorpixelsx::Int = 1000, detectorpixelsy::Int = 1000, ::Type{D} = Float32; temperature::Union{T,Unitful.Temperature} = T(OpticSim.GlassCat.TEMP_REF), pressure::T = T(OpticSim.GlassCat.PRESSURE_REF)) where {T<:Real,D<:Number}
        pr = prescription
        elements = Array{Union{Surface{T},CSGTree{T}},1}(undef, 0)
        (rows, cols) = size(pr)
        systemsemidiameter = zero(T)

        let sumstart = 1
            if pr[1, :Surface] == :Object
                sumstart = 2
            end
            firstelement = true

            for i in 1:(rows - 1)
                if pr[i, :Surface] == :Stop
                    stop = CircularAperture(T(pr[i, :SemiDiameter]), SVector{3,T}(0.0, 0.0, 1.0), SVector{3,T}(0.0, 0.0, T(-sum(pr[sumstart:(i - 1), :Thickness]))))
                    push!(elements, stop)
                else
                    material = pr[i, :Material]
                    frontreflection = zero(T)
                    backreflection = zero(T)
                    frontconic = zero(T)
                    backconic = zero(T)
                    frontaspherics = nothing
                    backaspherics = nothing
                    if material != OpticSim.GlassCat.Air && !isequal(material, missing)
                        semidiameter::T = max(pr[i, :SemiDiameter], pr[i + 1, :SemiDiameter])
                        if firstelement
                            systemsemidiameter = semidiameter
                            firstelement = false
                        end

                        vertex::T = -sum(pr[sumstart:(i - 1), :Thickness])
                        frontradius::T = pr[i, :Radius]
                        backradius::T = pr[i + 1, :Radius]
                        if "Reflectance" in names(prescription)
                            temp = pr[i, :Reflectance]
                            if !ismissing(temp)
                                frontreflection = temp
                            end
                            temp = pr[i + 1, :Reflectance]
                            if !ismissing(temp)
                                backreflection = temp
                            end
                        end
                        if "Conic" in names(prescription)
                            temp = pr[i, :Conic]
                            if !ismissing(temp)
                                frontconic = temp
                            end
                            temp = pr[i + 1, :Conic]
                            if !ismissing(temp)
                                backconic = temp
                            end
                        end
                        if "Aspherics" in names(prescription)
                            temp = pr[i, :Aspherics]
                            if !ismissing(temp)
                                frontaspherics = temp
                            end
                            temp = pr[i + 1, :Aspherics]
                            if !ismissing(temp)
                                backaspherics = temp
                            end
                        end

                        lastmaterial = pr[i - 1, :Material]
                        nextmaterial = pr[i + 1, :Material]

                        thickness = T(pr[i, :Thickness])

                        if frontaspherics !== nothing || backaspherics !== nothing
                            newelt = AsphericLens(material, vertex, frontradius, frontconic, frontaspherics, backradius, backconic, backaspherics, thickness, semidiameter, lastmaterial = lastmaterial, nextmaterial = nextmaterial, frontsurfacereflectance = frontreflection, backsurfacereflectance = backreflection)
                        elseif frontconic != zero(T) || backconic != zero(T)
                            newelt = ConicLens(material, vertex, frontradius, frontconic, backradius, backconic, thickness, semidiameter, lastmaterial = lastmaterial, nextmaterial = nextmaterial, frontsurfacereflectance = frontreflection, backsurfacereflectance = backreflection)
                        else
                            newelt = SphericalLens(material, vertex, frontradius, backradius, thickness, semidiameter, lastmaterial = lastmaterial, nextmaterial = nextmaterial, frontsurfacereflectance = frontreflection, backsurfacereflectance = backreflection)
                        end

                        push!(elements, newelt())
                    end
                end
            end

            indexofimage = findfirst(isequal(:Image), prescription[!, :Surface])
            imagesize = prescription[indexofimage, :SemiDiameter]
            vertextoimage = T(-sum(pr[sumstart:(indexofimage - 1), :Thickness]))
            imagerad = prescription[indexofimage, :Radius]

            if imagerad != zero(T) && imagerad != typemax(T)
                det = SphericalCap(abs(imagerad), NaNsafeasin(imagesize / abs(imagerad)), imagerad < 0 ? SVector{3,T}(0, 0, 1) : SVector{3,T}(0, 0, -1), SVector{3,T}(0, 0, vertextoimage), interface = opaqueinterface(T))
            else
                det = Rectangle(T(imagesize), T(imagesize), SVector{3,T}(0, 0, 1), SVector{3,T}(0, 0, vertextoimage), interface = opaqueinterface(T))
            end

            system = CSGOpticalSystem(OpticSim.LensAssembly(elements...), det, detectorpixelsx, detectorpixelsy, D, temperature = temperature, pressure = pressure)

            return new{T,typeof(system)}(system, prescription, systemsemidiameter)
        end
    end

    AxisymmetricOpticalSystem(prescription::DataFrame) = AxisymmetricOpticalSystem{Float64}(prescription)
end
export AxisymmetricOpticalSystem

Base.show(io::IO, a::AxisymmetricOpticalSystem) = print(io, a.prescription)
Base.copy(a::AxisymmetricOpticalSystem) = AxisymmetricOpticalSystem(a.prescription, size(detectorimage(a))..., temperature = temperature(a), pressure = pressure(a))

trace(system::AxisymmetricOpticalSystem{T,C}, r::OpticalRay{T,N}; trackrays::Union{Nothing,Vector{LensTrace{T,N}}} = nothing, test::Bool = false) where {T<:Real,N,C<:CSGOpticalSystem{T}} = trace(system.system, r, trackrays = trackrays, test = test)

"""
    semidiameter(system::AxisymmetricOpticalSystem{T}) -> T

Get the semidiameter of `system`, that is the semidiameter of the entrance pupil (i.e. first surface) of the system.
"""
semidiameter(a::AxisymmetricOpticalSystem) = a.semidiameter
assembly(system::AxisymmetricOpticalSystem) = assembly(system.system)
detector(system::AxisymmetricOpticalSystem) = detector(system.system)
detectorimage(system::AxisymmetricOpticalSystem) = detectorimage(system.system)
detectorsize(system::AxisymmetricOpticalSystem) = detectorsize(system.system)
resetdetector!(system::AxisymmetricOpticalSystem) = resetdetector!(system.system)
temperature(system::AxisymmetricOpticalSystem) = temperature(system.system)
pressure(system::AxisymmetricOpticalSystem) = pressure(system.system)
export semidiameter

######################################################################################################################

trace(system::AxisymmetricOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false, outpath::Union{Nothing,String} = nothing) where {T<:Real} = trace(system.system, raygenerator, printprog = printprog, test = test, outpath = outpath)

"""
    trace(system::OpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog = true, test = false)

Traces `system` with rays generated by `raygenerator` on a single thread.
Optionally the progress can be printed to the REPL.
If `test` is enabled then fresnel reflections are disabled and the power distribution will not be correct.
If `outpath` is specified then the result will be saved to this path.

Returns the detector image of the system.
"""
function trace(system::CSGOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false, outpath::Union{Nothing,String} = nothing) where {T<:Real}
    start_time = time()
    update_timesteps = 1000
    total_traced = 0
    for (k, r) in enumerate(raygenerator)
        if k % update_timesteps == 0
            total_traced += update_timesteps
            dif = round(time() - start_time, digits = 1)
            left = round((time() - start_time) * (length(sources) / total_traced - 1), digits = 1)
            if printprog
                print("\rTraced: ~ $t / $(length(sources))        Elapsed: $(dif)s        Left: $(left)s           ")
            end
        end
        trace(system, r, test = test)
    end
    numrays = length(raygenerator)
    tracetime = round(time() - all_start_time)
    if printprog
        print("\rFinished tracing $numrays rays in $(tracetime)s")
        if tracetime != 0.0
            print(",  $(Int32(round(numrays/tracetime))) rays per second")
        end
        println()
    end

    det = detectorimage(system)

    if outpath !== nothing
        save(outpath, colorview(Gray, det ./ maximum(det)))
    end

    return det
end

function traceMT(system::AxisymmetricOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false, outpath::Union{Nothing,String} = nothing) where {T<:Real}
    traceMT(system.system, raygenerator, printprog = printprog, test = test, outpath = outpath)
end

"""
    traceMT(system::OpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog = true, test = false)

Traces `system` with rays generated by `raygenerator` using as many threads as possible.
Optionally the progress can be printed to the REPL.
If `test` is enabled then fresnel reflections are disabled and the power distribution will not be correct.
If `outpath` is specified then the result will be saved to this path.

Returns the accumulated detector image from all threads.
"""
function traceMT(system::CSGOpticalSystem{T,S}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false, outpath::Union{Nothing,String} = nothing) where {T<:Real,S<:Number}
    if printprog
        println("Initialising...")
    end

    all_start_time = time()
    N = Threads.nthreads()
    copies = Vector{Tuple{CSGOpticalSystem{T},OpticalRayGenerator{T}}}(undef, N)
    for i in eachindex(copies)
        copies[i] = (copy(system), copy(raygenerator))
    end

    if printprog
        println("Tracing...")
    end
    total_traced = Threads.Atomic{Int}(0)
    update_timesteps = 1000
    Threads.@threads for sys in copies
        (system, sources) = sys
        i = Threads.threadid()
        len, rem = divrem(length(sources), N)
        # not enough iterations for all the threads
        if len == 0
            if i > rem
                return
            end
            len, rem = 1, 0
        end
        # compute this thread's iterations
        f = firstindex(sources) + ((i - 1) * len)
        l = f + len - 1
        # distribute remaining iterations evenly
        if rem > 0
            if i <= rem
                f = f + (i - 1)
                l = l + i
            else
                f = f + rem
                l = l + rem
            end
        end
        if i == 1
            start_time = time()
            for k in f:l
                if k % update_timesteps == 0
                    Threads.atomic_add!(total_traced, update_timesteps)
                    dif = round(time() - start_time, digits = 1)
                    t = total_traced[]
                    left = round((time() - start_time) * (length(sources) / t - 1), digits = 1)
                    if printprog
                        print("\rTraced: ~ $t / $(length(sources))        Elapsed: $(dif)s        Left: $(left)s           ")
                    end
                end
                trace(system, generateray(sources, k), test = test)
            end
        else
            for k in f:l
                if k % update_timesteps == 0
                    Threads.atomic_add!(total_traced, update_timesteps)
                end
                trace(system, generateray(sources, k), test = test)
            end
        end
    end

    if printprog
        println("\rAccumulating images...                                                                       ")
    end

    det = OpticSim.detectorimage(copies[1][1])
    # sum all the copies
    @inbounds for i in 2:N
        OpticSim.sum!(det, OpticSim.detectorimage(copies[i][1]))
    end

    numrays = length(raygenerator)
    tracetime = round(time() - all_start_time)
    if printprog
        print("\rFinished tracing $numrays rays in $(tracetime)s,  ")
        if tracetime != 0.0
            print("$(Int32(round(numrays/tracetime))) rays per second")
        end
        println()
    end

    if outpath !== nothing
        save(outpath, colorview(Gray, det ./ maximum(det)))
    end

    return det
end

function tracehitsMT(system::AxisymmetricOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false) where {T<:Real}
    tracehitsMT(system.system, raygenerator, printprog = printprog, test = test)
end

"""
    tracehitsMT(system::OpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog = true, test = false)

Traces `system` with rays generated by `raygenerator` using as many threads as possible.
Optionally the progress can be printed to the REPL.
If `test` is enabled then fresnel reflections are disabled and the power distribution will not be correct.

Returns a list of [`LensTrace`](@ref)s which hit the detector, accumulated from all threads.
"""
function tracehitsMT(system::CSGOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false) where {T<:Real}
    if printprog
        println("Initialising...")
    end

    all_start_time = time()
    N = Threads.nthreads()
    copies = Vector{Tuple{CSGOpticalSystem{T},OpticalRayGenerator{T}}}(undef, N)
    results = Vector{Vector{LensTrace{T,3}}}(undef, N)
    for i in eachindex(copies)
        copies[i] = (copy(system), copy(raygenerator))
        results[i] = Vector{LensTrace{T,3}}(undef, 0)
    end

    if printprog
        println("Tracing...")
    end
    total_traced = Threads.Atomic{Int}(0)
    update_timesteps = 1000
    Threads.@threads for sys in copies
        (system, sources) = sys
        i = Threads.threadid()
        len, rem = divrem(length(sources), N)
        # not enough iterations for all the threads
        if len == 0
            if i > rem
                return
            end
            len, rem = 1, 0
        end
        # compute this thread's iterations
        f = firstindex(sources) + ((i - 1) * len)
        l = f + len - 1
        # distribute remaining iterations evenly
        if rem > 0
            if i <= rem
                f = f + (i - 1)
                l = l + i
            else
                f = f + rem
                l = l + rem
            end
        end
        if i == 1
            start_time = time()
            for k in f:l
                if k % update_timesteps == 0
                    Threads.atomic_add!(total_traced, update_timesteps)
                    dif = round(time() - start_time, digits = 1)
                    t = total_traced[]
                    left = round((time() - start_time) * (length(sources) / t - 1), digits = 1)
                    if printprog
                        print("\rTraced: ~ $t / $(length(sources))        Elapsed: $(dif)s        Left: $(left)s           ")
                    end
                end
                lt = trace(system, generateray(sources, k), test = test)
                if lt !== nothing
                    push!(results[i], lt)
                end
            end
        else
            for k in f:l
                if k % update_timesteps == 0
                    Threads.atomic_add!(total_traced, update_timesteps)
                end
                lt = trace(system, generateray(sources, k), test = test)
                if lt !== nothing
                    push!(results[i], lt)
                end
            end
        end
    end

    if printprog
        println("\rAccumulating results...                                                                       ")
    end

    accumres = results[1]
    # sum all the copies
    @inbounds for i in 2:N
        accumres = vcat(accumres, results[i])
    end

    numrays = length(raygenerator)
    tracetime = round(time() - all_start_time)
    if printprog
        print("\rFinished tracing $numrays rays in $(tracetime)s")
        if tracetime != 0.0
            print(",  $(Int32(round(numrays/tracetime))) rays per second")
        end
        println()
    end

    return accumres
end

tracehits(system::AxisymmetricOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false) where {T<:Real} = tracehits(system.system, raygenerator, printprog = printprog, test = test)

"""
    tracehits(system::OpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog = true, test = false)

Traces `system` with rays generated by `raygenerator` on a single thread.
Optionally the progress can be printed to the REPL.
If `test` is enabled then fresnel reflections are disabled and the power distribution will not be correct.

Returns a list of [`LensTrace`](@ref)s which hit the detector.
"""
function tracehits(system::CSGOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; printprog::Bool = true, test::Bool = false) where {T<:Real}
    start_time = time()
    update_timesteps = 1000
    total_traced = 0
    res = Vector{LensTrace{T,3}}(undef, 0)
    for (k, r) in enumerate(raygenerator)
        if k % update_timesteps == 0
            total_traced += update_timesteps
            dif = round(time() - start_time, digits = 1)
            left = round((time() - start_time) * (length(sources) / total_traced - 1), digits = 1)
            if printprog
                print("\rTraced: ~ $t / $(length(sources))        Elapsed: $(dif)s        Left: $(left)s           ")
            end
        end
        lt = trace(system, r, test = test)
        if lt !== nothing
            push!(res, lt)
        end
    end
    numrays = length(raygenerator)
    tracetime = round(time() - all_start_time)
    if printprog
        print("\rFinished tracing $numrays rays in $(tracetime)s")
        if tracetime != 0.0
            print(",  $(Int32(round(numrays/tracetime))) rays per second")
        end
        println()
    end

    return res
end

export trace, traceMT, tracehits, tracehitsMT


text:
# ================================
# Flowpipe
# ================================

"""
    Flowpipe{N, RT<:AbstractReachSet{N}, VRT<:AbstractVector{RT}} <: AbstractFlowpipe

Type that wraps a flowpipe, which is an iterable collection of reach-sets that
behaves like their set union.

### Fields

- `Xk`  -- array of reach-sets
- `ext` -- extension dictionary; field used by extensions

### Notes

The dimension of the flowpipe corresponds to the dimension of the underlying
reach-sets; in this type, it is is assumed that the dimension is the same for
the different reach-sets.
"""
struct Flowpipe{N, RT<:AbstractReachSet{N}, VRT<:AbstractVector{RT}} <: AbstractFlowpipe
    Xk::VRT
    ext::Dict{Symbol, Any}
end

# getter functions
@inline array(fp::Flowpipe) = fp.Xk
@inline flowpipe(fp::Flowpipe) = fp

# constructor from empty extension dictionary
function Flowpipe(Xk::AbstractVector{RT}) where {N, RT<:AbstractReachSet{N}}
    return Flowpipe(Xk, Dict{Symbol, Any}())
end

# undef initializer given a set type
function Flowpipe(::UndefInitializer, ST::Type{<:LazySet{N}}, k::Int) where {N}
    return Flowpipe(Vector{ReachSet{N, ST}}(undef, k), Dict{Symbol, Any}())
end

# undef initializer given a reach-set type
function Flowpipe(::UndefInitializer, RT::Type{<:AbstractReachSet}, k::Int)
    return Flowpipe(Vector{RT}(undef, k), Dict{Symbol, Any}())
end

# undef initializer given a continuous post-operator
function Flowpipe(::UndefInitializer, cpost::AbstractContinuousPost, k::Int)
    RT = rsetrep(cpost)
    return Flowpipe(Vector{RT}(undef, k), Dict{Symbol, Any}())
end

# constructor from a single reach-set
Flowpipe(R::AbstractReachSet) = Flowpipe([R])

function Base.similar(fp::Flowpipe{N, RT, VRT}) where {N, RT, VRT}
   return Flowpipe(VRT())
end

Base.IndexStyle(::Type{<:Flowpipe}) = IndexLinear()
Base.eltype(::Flowpipe{N, RT}) where {N, RT} = RT
Base.size(fp::Flowpipe) = (length(fp.Xk),)
Base.view(fp::Flowpipe, args...) = view(fp.Xk, args...)
Base.push!(fp::Flowpipe, args...) = push!(fp.Xk, args...)
Base.keys(fp::Flowpipe) = keys(fp.Xk)
Base.append!(fp::Flowpipe{N, RT, VT}, v::VT) where {N, RT, VT} = append!(fp.Xk, v)

numtype(::Flowpipe{N}) where {N} = N
setrep(fp::Flowpipe{N, RT}) where {N, RT} = setrep(RT)
setrep(::Type{<:Flowpipe{N, RT}}) where {N, RT} = setrep(RT)
rsetrep(fp::Flowpipe{N, RT}) where {N, RT} = RT
rsetrep(::Type{<:Flowpipe{N, RT}}) where {N, RT} = RT
numrsets(fp::Flowpipe) = length(fp)

# getter functions for hybrid systems

"""
    location(F::Flowpipe)

Return the location of a flowpipe within a hybrid system, or `missing` if it is
not defined.

### Input

- `F` -- flowpipe

### Output

The `:loc_id` value of stored in the flowpipe's extension field.
"""
function location(fp::Flowpipe)
    return get(fp.ext, :loc_id, missing)
end

# evaluate a flowpipe at a given time point: gives a reach set
# here it would be useful to layout the times contiguously in a vector
# (see again array of struct vs struct of array)
function (fp::Flowpipe)(t::Number)
    Xk = array(fp)
    @inbounds for (i, X) in enumerate(Xk)
        if t ∈ tspan(X) # exit on the first occurrence
            if i < length(Xk) && t ∈ tspan(Xk[i+1])
                return view(Xk, i:i+1)
            else
                return X
            end
        end
    end
    throw(ArgumentError("time $t does not belong to the time span, " *
                        "$(tspan(fp)), of the given flowpipe"))
end

# evaluate a flowpipe at a given time interval: gives possibly more than one reach set
# i.e. first and last sets and those in between them
function (fp::Flowpipe)(dt::TimeInterval)
    idx = Vector{Int}()
    α = inf(dt)
    β = sup(dt)
    Xk = array(fp)
    for (i, X) in enumerate(Xk)
        if !isempty(tspan(X) ∩ dt)
            push!(idx, i)
        end
    end

    if isempty(idx)
        throw(ArgumentError("the time interval $dt does not intersect the time span, " *
                            "$(tspan(fp)), of the given flowpipe"))
    end
    return view(Xk, idx)
end

# concrete projection of a flowpipe along variables `vars`
function project(fp::Flowpipe, vars::NTuple{D, T}) where {D, T<:Integer}
    Xk = array(fp)
    πfp = map(X -> project(X, vars), Xk)
    return Flowpipe(πfp, fp.ext)
end

project(fp::Flowpipe, vars::Int) = project(fp, (vars,))
project(fp::Flowpipe, vars::AbstractVector{<:Int}) = project(fp, Tuple(vars))
project(fp::Flowpipe; vars) = project(fp, Tuple(vars))
project(fp::Flowpipe, i::Int, vars) = project(fp[i], vars)

# concrete linear map of a flowpipe for a given matrix
function linear_map(M::AbstractMatrix, fp::Flowpipe)
    Xk = array(fp)
    πfp = map(X -> linear_map(M, X), Xk)
    return Flowpipe(πfp, fp.ext)
end

"""
    shift(fp::Flowpipe{N, <:AbstractReachSet}, t0::Number) where {N}

Return the time-shifted flowpipe by the given number.

### Input

- `fp` -- flowpipe
- `t0` -- time shift

### Output

A new flowpipe such that the time-span of each constituent reach-set has been
shifted by `t0`.

### Notes

See also `Shift` for the lazy counterpart.
"""
function shift(fp::Flowpipe{N, <:AbstractReachSet}, t0::Number) where {N}
    return Flowpipe([shift(X, t0) for X in array(fp)], fp.ext)
end

"""
    convexify(fp::Flowpipe{N, <:AbstractLazyReachSet}) where {N}

Return a reach-set representing the convex hull array of the flowpipe.

### Input

- `fp` -- flowpipe

### Output

A reach-set that contains the convex hull array, `ConvexHullArray`, of the given
flowpipe.

### Notes

The time span of this reach-set is the same as the time-span of the flowpipe.

This function allocates an array to store the sets of the flowpipe.
"""
function convexify(fp::Flowpipe{N, <:AbstractLazyReachSet}) where {N}
    Y = ConvexHullArray([set(X) for X in array(fp)])
    return ReachSet(Y, tspan(fp))
end

"""
    convexify(fp::AbstractVector{<:AbstractLazyReachSet{N}}) where {N}

Return a reach-set representing the convex hull array of the array of the array of
reach-sets.

### Input

- `fp` -- array of reach-sets

### Output

A reach-set that contains the convex hull array, `ConvexHullArray`, of the given
flowpipe.

### Notes

The time span of this reach-set corresponds to the minimum (resp. maximum) of the
time span of each reach-set in `fp`.

This function allocates an array to store the sets of the flowpipe.

The function doesn't assume that the reach-sets are time ordered.
"""
function convexify(fp::AbstractVector{<:AbstractLazyReachSet{N}}) where {N}
    Y = ConvexHullArray([set(X) for X in fp])
    ti = minimum(tstart, fp)
    tf = maximum(tend, fp)
    return ReachSet(Y, TimeInterval(ti, tf))
end

# the dimension of sparse flowpipes is known in the type
function LazySets.dim(::Flowpipe{N, SparseReachSet{N, ST, D}}) where {N, ST, D}
    return D
end

tstart(F::Flowpipe, arr::UnitRange) = tstart(view(array(F), arr), contiguous=true)
tstart(F::Flowpipe, arr::AbstractVector) = tstart(view(array(F), arr))
tend(F::Flowpipe, arr::UnitRange) = tend(view(array(F), arr), contiguous=true)
tend(F::Flowpipe, arr::AbstractVector) = tend(view(array(F), arr))
tspan(F::Flowpipe, arr::UnitRange) = tspan(view(array(F), arr), contiguous=true)
tspan(F::Flowpipe, arr::AbstractVector) = tspan(view(array(F), arr))

# further setops
LazySets.isdisjoint(F::Flowpipe{N, <:AbstractLazyReachSet}, Y::LazySet) where {N} = all(X -> _is_intersection_empty(X, Y), array(F))
Base.:⊆(F::Flowpipe, X::LazySet) = all(R ⊆ X for R in F)
Base.:⊆(F::Flowpipe, Y::AbstractLazyReachSet) = all(R ⊆ set(Y) for R in F)

# lazy projection of a flowpipe
function Projection(F::Flowpipe, vars::NTuple{D, T}) where {D, T<:Integer}
    Xk = array(F)
    out = map(X -> Projection(X, vars), Xk)
    return Flowpipe(out)
end
Projection(F::Flowpipe; vars) = Projection(F, Tuple(vars))
Projection(F::Flowpipe, vars::AbstractVector{M}) where {M<:Integer} = Projection(F, Tuple(vars))

# membership test
function ∈(x::AbstractVector{N}, fp::Flowpipe{N, <:AbstractLazyReachSet{N}}) where {N}
    return any(R -> x ∈ set(R), array(fp))
end

function ∈(x::AbstractVector{N}, fp::VT) where {N, RT<:AbstractLazyReachSet{N}, VT<:AbstractVector{RT}}
    return any(R -> x ∈ set(R), fp)
end

function LazySets.linear_map(M, fp::Flowpipe)
    out = [linear_map(M, R) for R in fp]
    return Flowpipe(out, fp.ext)
end

function LazySets.affine_map(M, b, fp::Flowpipe)
    out = [reconstruct(R, affine_map(M, set(R), b)) for R in fp]
    return Flowpipe(out, fp.ext)
end

function LazySets.translate(fp::Flowpipe, v::AbstractVector)
    out = [reconstruct(R, translate(set(R), v)) for R in fp]
    return Flowpipe(out, fp.ext)
end

# lazy intersection
# use the flag `filter` (optional, default: `true`) to remove empty sets
@commutative function Base.:∩(fp::Flowpipe, X::LazySet; filter=true)
    out = [reconstruct(R, set(R) ∩ X) for R in fp]

    if filter
        # keep the sets which are not empty
        filter!(!isempty, out)
    end
    return Flowpipe(out, fp.ext)
end

# concrete intersection
# use the flag `filter` (optional, default: `true`) to remove empty sets
@commutative function LazySets.intersection(fp::Flowpipe, X::LazySet; filter=true)
    out = [reconstruct(R, intersection(set(R), X)) for R in fp]

    if filter
        # keep the sets which are not empty
        filter!(!isempty, out)
    end
    return Flowpipe(identity.(out), fp.ext)
end

# --------------------------------------------
# Specialized methods for template flowpipes
# --------------------------------------------

# assumes that the first reach-set is representative
function support_function_matrix(fp::Flowpipe{N, <:TemplateReachSet}) where {N}
    return support_function_matrix(first(fp))
end

# it is assumed that rows = (idx_pos_dir, idx_neg_dir) is such that each integer
# idx_pos_dir and idx_neg_dir refers to the row of the support function matrix
# correponding to a positive (resp. negative) direction
function flatten(fp::Flowpipe{N, <:TemplateReachSet}, rows=(1, 2)) where {N}
    # get the matrix of support function evaluations
    mat = support_function_matrix(fp)

    @assert size(mat, 1) ≥ 2 || throw(ArgumentError("the number of rows of the support function matrix should be at least 2, got $(size(mat, 1))"))
    @assert length(rows) == 2 || throw(ArgumentError("expected the number of rows of the support function matrix to be 2, got $(length(rows))"))
    idx_pos_dir = rows[1]
    idx_neg_dir = rows[2]

    RT = ReachSet{N, Interval{N, IA.Interval{N}}}
    out = Vector{RT}(undef, length(fp))

    @inbounds for (k, Rk) in enumerate(fp)
        Ik = Interval(-mat[idx_neg_dir, k], mat[idx_pos_dir, k])
        dtk = tspan(Rk)
        out[k] = ReachSet(Ik, dtk)
    end
    return Flowpipe(out, fp.ext)
end


text:
# Modelo de torre de refrigeracao
type Torre
	Torre()=begin
		new(
			CalorEspecifico(),
			VazaoMassica(),
			Potencia(Dict{Symbol,Any}(
				:Symbol=>"_{in1}",
				:PosX=>0.1,
				:PosY=>0.5
			)),
			Dif_Temp(),
			Dif_Temp(),
			Dif_Temp(),
			Temperatura(),
			Temperatura(),
			Temperatura(),
			Temperatura(),
			CoefGlobal_area(),
			[
				:(DTar = Tar_h - Tar_c),
				:(DTh = Th - Tar_h),
				:(DTc = Tc - Tar_c),
				:(F * cpa * (Th - Tc) = Qin),
				:(Uat * (DTh - DTc) = Qin * ln(abs(DTh/DTc))),
			],
			[
				"","","","","",
			],
			[:cpa,],
			[:F,:Qin,:DTh,:DTc,:DTar,:Th,:Tc,:Tar_c,:Tar_h,:Uat,]
		)
	end
	cpa::CalorEspecifico
	F::VazaoMassica
	Qin::Potencia
	DTh::Dif_Temp
	DTc::Dif_Temp
	DTar::Dif_Temp
	Th::Temperatura
	Tc::Temperatura
	Tar_c::Temperatura
	Tar_h::Temperatura
	Uat::CoefGlobal_area
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Torre
function setEquationFlow(in::Torre)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	#	Uat * 0.5 * (DTh + DTc) = Qin;

end
function atributes(in::Torre,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/torreresf"
	drive!(fields,_)
	return fields
end
Torre(_::Dict{Symbol,Any})=begin
	newModel=Torre()
	newModel.attributes=atributes(newModel,_)
	newModel
end