text:
const Indexable = Union{AbstractVector,Tuple}

# Point
geomtype(::AbstractVector{T}) where {T <: Real} = Point()
geomtype(::Tuple{T,U}) where {T,U <: Real} = Point()
geomtype(::Tuple{T,U,V}) where {T,U,V <: Real} = Point()
ncoord(::AbstractPoint, geom::Indexable) = length(geom)
getcoord(::AbstractPoint, geom::Indexable, i::Integer) = geom[i]

# LineString
ncoord(::AbstractLineString, geom::Indexable) =
    ncoord(Point, getpoint(geom, 1))
npoint(::AbstractLineString, geom::Indexable) = length(geom)
getpoint(::AbstractLineString, geom::Indexable, i::Integer) = geom[i]

# Polygon
ncoord(::AbstractPolygon, geom::Indexable) =
    ncoord(LineString, getexterior(geom))
getexterior(::AbstractPolygon, geom::Indexable) = geom[1]
nhole(::AbstractPolygon, geom::Indexable) = length(geom) - 1
gethole(::AbstractPolygon, geom::Indexable, i::Integer) = geom[i + 1]

# MultiPoint
ncoord(::AbstractMultiPoint, geom::Indexable) =
    ncoord(Point, getpoint(geom, 1))
npoint(::AbstractMultiPoint, geom::Indexable) = length(geom)
getpoint(::AbstractMultiPoint, geom::Indexable, i::Integer) = geom[i]

# MultiLineString
ncoord(::AbstractMultiLineString, geom::Indexable) =
    ncoord(LineString, getlinestring(geom, 1))
nlinestring(::AbstractMultiLineString, geom::Indexable) = length(geom)
getlinestring(::AbstractMultiLineString, geom::Indexable, i::Integer) =
    geom[i]

# MultiPolygon
ncoord(::AbstractMultiPolygon, geom::Indexable) =
    ncoord(Polygon, getpolygon(geom, 1))
npolygon(::AbstractMultiPolygon, geom::Indexable) = length(geom)
getpolygon(::AbstractMultiPolygon, geom::Indexable, i::Integer) = geom[i]

# GeometryCollection
ncoord(::AbstractGeometryCollection, geom::Indexable) =
    ncoord(geomtype(geom), getgeom(geom, 1))
ngeom(::AbstractGeometryCollection, collection::Indexable) =
    length(collection)
getgeom(::AbstractGeometryCollection, collection::Indexable, i::Integer) =
    collection[i]


text:
substitute(σ) = Base.Fix2(replace, σ)
replace(t::AbstractTerm, σ) = map(substitute(σ), t)
replace(t::Union{Variable,AbstractTerm}) = Base.Fix1(replace, t)

function many_matchers(ps::Vector, V)
    matchers = similar(ps, Union{AbstractMatcher,Variable})
    for (i, p) ∈ enumerate(ps)
        matchers[i], V::Set{Variable} = matcher(p, V)
    end
    return matchers, V
end

function compatible(p::Dict, q::Dict)
    length(p) < length(q) && ((p, q) = (q, p))
    for (x, t) ∈ q
        haskey(p, x) && (p[x] == t || return false)
    end
    return true
end

@inline _aiterate(p) = (p, ())
@inline function _aiterate(p, iter1, rest...)
    restres = _aiterate(p, rest...)
    restres === nothing && return nothing
    (P, states) = restres

    thisres = iterate(Matches(P, iter1))
    while thisres === nothing
        restres = _aiterate1(p, rest, states)
        restres === nothing && return nothing
        (P, states) = restres

        thisres = iterate(Matches(P, iter1))
    end
    (P′, state) = thisres

    return (P′, ((P, state), states...))
end

@inline _aiterate1(p, ::Tuple{}, ::Tuple{}) = nothing
@inline function _aiterate1(p, iters, states)
    iter1 = first(iters)
    (P, state1) = first(states)
    next = iterate(Matches(P, iter1), state1)
    reststates = Base.tail(states)

    if next === nothing
        restnext = _aiterate1(p, Base.tail(iters), reststates)
        restnext === nothing && return nothing
        (P, reststates) = restnext

        next = iterate(Matches(P, iter1))
        next === nothing && return nothing
    end

    (P′, state) = next
    return (P′, ((P, state), reststates...))
end

struct LazyMap{T,F}
    f::F
    iter::T
end
(Base.IteratorSize(::Type{<:LazyMap{T}}) where T) = Base.IteratorSize(T)
Base.length(iter::LazyMap) = length(iter.iter)
function Base.iterate(iter::LazyMap)
    next = iterate(iter.iter)
    next === nothing && return nothing
    (result, state) = next
    return (iter.f(result), state)
end
function Base.iterate(iter::LazyMap, state)
    next = iterate(iter.iter, state)
    next === nothing && return nothing
    (result, state) = next
    return (iter.f(result), state)
end

struct LazyFlatten{T}
    iter::T
end
(Base.IteratorSize(::Type{<:LazyFlatten{T}}) where T) = Base.IteratorSize(T)
Base.length(iter::LazyFlatten) = length(iter.iter)
function _iterate_lazyflatten(next)
    while next !== nothing
        (inner, state) = next

        inext = iterate(inner)
        if inext === nothing
            (inner, state) = iterate(iter.iter, state)
            continue
        end

        (result, istate) = inext
        return result, (state, (inner, istate))
    end

    return nothing
end
function Base.iterate(iter::LazyFlatten)
    next = iterate(iter.iter)
    return _iterate_lazyflatten(next)
end
function Base.iterate(iter::LazyFlatten, (state, (inner, istate)))
    inext = iterate(inner, istate)
    if inext !== nothing
        (result, istate) = inext
        return result, (state, (inner, istate))
    end

    next = iterate(iter.iter, state)
    return _iterate_lazyflatten(next)
end

ext: jl


text:
ExodusConstant = Int64

# TODO: cleanup up below
#
const MAX_LINE_LENGTH = 80
const MAX_STR_LENGTH = 32

const cpu_word_size = Ref{Int64}(sizeof(Float64))
const IO_word_size = Ref{Int64}(8)

# TODO: make this be read in from the OS or something like that
#
const version_number_1 = 8
const version_number_2 = 15
const version_number = 8.15

# TODO: add the entirity from the exoudsII header

const EX_WRITE         = 0x0001          #< ex_open(): open existing file for appending. */
const EX_READ          = 0x0002          #< ex_open(): open file for reading (default) */

const EX_NOCLOBBER     = 0x0004          #< Don't overwrite existing database, default */
const EX_CLOBBER       = 0x0008          #< Overwrite existing database if it exists */
const EX_NORMAL_MODEL  = 0x0010          #< disable mods that permit storage of larger models */
const EX_64BIT_OFFSET  = 0x0020          #< enable mods that permit storage of larger models */
const EX_LARGE_MODEL   = EX_64BIT_OFFSET #< enable mods that permit storage of larger models */
const EX_64BIT_DATA    = 0x400000        #< CDF-5 format: classic model but 64 bit dimensions and sizes */
const EX_NETCDF4       = 0x0040          #< use the hdf5-based netcdf4 output */
const EX_NOSHARE       = 0x0080          #< Do not open netcdf file in "share" mode */
const EX_SHARE         = 0x0100          #< Do open netcdf file in "share" mode */
const EX_NOCLASSIC     = 0x0200          #< Do not force netcdf to classic mode in netcdf4 mode */

const EX_DISKLESS      = 0x100000        #< Experimental */
const EX_MMAP          = 0x200000        #< Experimental */

const EX_MAPS_INT64_DB = 0x0400          #< All maps (id, order, ...) store int64_t values */
const EX_IDS_INT64_DB  = 0x0800          #< All entity ids (sets, blocks, maps) are int64_t values */
const EX_BULK_INT64_DB = 0x1000          #< All integer bulk data (local indices, counts, maps); not ids                         \*/

# TODO: figure out how to define EX_ALL_INT64_DB

const EX_MAPS_INT64_API = 0x2000         #< All maps (id, order, ...) store int64_t values */
const EX_IDS_INT64_API  = 0x4000         #< All entity ids (sets, blocks, maps) are int64_t values */
const EX_BULK_INT64_API = 0x8000         #< All integer bulk data (local indices, counts, maps); not ids */
const EX_INQ_INT64_API  = 0x10000        #< Integers passed to/from ex_inquire() are int64_t */

# TODO: figure out how to define EX_ALL_INT64_API

# Parallel IO mode flags
#
const EX_MPIIO          = 0x20000
const EX_MPIPOSIX       = 0x40000        #< \deprecated As of libhdf5 1.8.13. */
const EX_PNETCDF        = 0x80000

# ex_inquire() stuff
#
@enum ExInquiry begin
    EX_INQ_FILE_TYPE                  = 1  #< EXODUS file type (deprecated) */
    EX_INQ_API_VERS                   = 2  #< API version number (float) */
    EX_INQ_DB_VERS                    = 3  #< database version number (float) */
    EX_INQ_TITLE                      = 4  #< database title. MAX_LINE_LENGTH+1 char* size */
    EX_INQ_DIM                        = 5  #< number of dimensions */
    EX_INQ_NODES                      = 6  #< number of nodes    */
    EX_INQ_ELEM                       = 7  #< number of elements */
    EX_INQ_ELEM_BLK                   = 8  #< number of element blocks */
    EX_INQ_NODE_SETS                  = 9  #< number of node sets*/
    EX_INQ_NS_NODE_LEN                = 10 #< length of node set node list */
    EX_INQ_SIDE_SETS                  = 11 #< number of side sets*/
    EX_INQ_SS_NODE_LEN                = 12 #< length of side set node list */
    EX_INQ_SS_ELEM_LEN                = 13 #< length of side set element list */
    EX_INQ_QA                         = 14 #< number of QA records */
    EX_INQ_INFO                       = 15 #< number of info records */
    EX_INQ_TIME                       = 16 #< number of time steps in the database */
    EX_INQ_EB_PROP                    = 17 #< number of element block properties */
    EX_INQ_NS_PROP                    = 18 #< number of node set properties */
    EX_INQ_SS_PROP                    = 19 #< number of side set properties */
    EX_INQ_NS_DF_LEN                  = 20 #< length of node set distribution factor list*/
    EX_INQ_SS_DF_LEN                  = 21 #< length of side set distribution factor list*/
    EX_INQ_LIB_VERS                   = 22 #< API Lib vers number (float) */
    EX_INQ_EM_PROP                    = 23 #< number of element map properties */
    EX_INQ_NM_PROP                    = 24 #< number of node map properties */
    EX_INQ_ELEM_MAP                   = 25 #< number of element maps */
    EX_INQ_NODE_MAP                   = 26 #< number of node maps*/
    EX_INQ_EDGE                       = 27 #< number of edges    */
    EX_INQ_EDGE_BLK                   = 28 #< number of edge blocks */
    EX_INQ_EDGE_SETS                  = 29 #< number of edge sets   */
    EX_INQ_ES_LEN                     = 30 #< length of concat edge set edge list       */
    EX_INQ_ES_DF_LEN                  = 31 #< length of concat edge set dist factor list*/
    EX_INQ_EDGE_PROP                  = 32 #< number of properties stored per edge block    */
    EX_INQ_ES_PROP                    = 33 #< number of properties stored per edge set      */
    EX_INQ_FACE                       = 34 #< number of faces */
    EX_INQ_FACE_BLK                   = 35 #< number of face blocks */
    EX_INQ_FACE_SETS                  = 36 #< number of face sets */
    EX_INQ_FS_LEN                     = 37 #< length of concat face set face list */
    EX_INQ_FS_DF_LEN                  = 38 #< length of concat face set dist factor list*/
    EX_INQ_FACE_PROP                  = 39 #< number of properties stored per face block */
    EX_INQ_FS_PROP                    = 40 #< number of properties stored per face set */
    EX_INQ_ELEM_SETS                  = 41 #< number of element sets */
    EX_INQ_ELS_LEN                    = 42 #< length of concat element set element list       */
    EX_INQ_ELS_DF_LEN                 = 43 #< length of concat element set dist factor list*/
    EX_INQ_ELS_PROP                   = 44 #< number of properties stored per elem set      */
    EX_INQ_EDGE_MAP                   = 45 #< number of edge maps                     */
    EX_INQ_FACE_MAP                   = 46 #< number of face maps                     */
    EX_INQ_COORD_FRAMES               = 47 #< number of coordinate frames */
    EX_INQ_DB_MAX_ALLOWED_NAME_LENGTH = 48 #< size of MAX_NAME_LENGTH dimension on database */
    EX_INQ_DB_MAX_USED_NAME_LENGTH    = 49 #< size of MAX_NAME_LENGTH dimension on database */
    EX_INQ_MAX_READ_NAME_LENGTH       = 50 #< client-specified max size of returned names */
    EX_INQ_DB_FLOAT_SIZE              = 51 #< size of floating-point values stored on database */
    EX_INQ_NUM_CHILD_GROUPS           = 52 #< number of groups contained in this (exoid) group */
    EX_INQ_GROUP_PARENT               = 53 #< id of parent of this (exoid) group; returns exoid if at root */
    EX_INQ_GROUP_ROOT                 = 54 #< id of root group "/" of this (exoid) group; returns exoid if at root */
    EX_INQ_GROUP_NAME_LEN             = 55 #< length of name of group exoid */
    EX_INQ_GROUP_NAME                 = 56 #< name of group exoid. "/" returned for root group (char* GROUP_NAME_LEN+1 size) */
    EX_INQ_FULL_GROUP_NAME_LEN        = 57 #< length of full path name of this (exoid) group */
    EX_INQ_FULL_GROUP_NAME            = 58 #< full "/"-separated path name of this (exoid) group */
    EX_INQ_THREADSAFE                 = 59 #< Returns 1 if library is thread-safe; 0 otherwise */
    EX_INQ_ASSEMBLY                   = 60 #< number of assemblies */
    EX_INQ_BLOB                       = 61 #< number of blobs */
    EX_INQ_NUM_NODE_VAR               = 62 #< number of nodal variables */
    EX_INQ_NUM_EDGE_BLOCK_VAR         = 63 #< number of edge block variables */
    EX_INQ_NUM_FACE_BLOCK_VAR         = 64 #< number of face block variables */
    EX_INQ_NUM_ELEM_BLOCK_VAR         = 65 #< number of element block variables */
    EX_INQ_NUM_NODE_SET_VAR           = 66 #< number of node set variables */
    EX_INQ_NUM_EDGE_SET_VAR           = 67 #< number of edge set variables */
    EX_INQ_NUM_FACE_SET_VAR           = 68 #< number of face set variables */
    EX_INQ_NUM_ELEM_SET_VAR           = 69 #< number of element set variables */
    EX_INQ_NUM_SIDE_SET_VAR           = 70 #< number of sideset variables */
    EX_INQ_NUM_GLOBAL_VAR             = 71 #< number of global variables */
    EX_INQ_INVALID                    = -1
end

# TODO: finish below
# options type
#
# @enum ExOptionType begin

# end

# entities
#
@enum ExEntityType begin
    EX_NODAL      = 14 #< nodal "block" for variables*/
    # TODO: figure out how to enable below enum
    # TODO: julia enum won't let you define something
    # TODO: with the same value
    # EX_NODE_BLOCK = 14 #< alias for EX_NODAL         */
    EX_NODE_SET   = 2  #< node set property code     */
    EX_EDGE_BLOCK = 6  #< edge block property code   */
    EX_EDGE_SET   = 7  #< edge set property code     */
    EX_FACE_BLOCK = 8  #< face block property code   */
    EX_FACE_SET   = 9  #< face set property code     */
    EX_ELEM_BLOCK = 1  #< element block property code*/
    EX_ELEM_SET   = 10 #< face set property code     */

    EX_SIDE_SET   = 3  #< side set property code     */

    EX_ELEM_MAP   = 4  #< element map property code  */
    EX_NODE_MAP   = 5  #< node map property code     */
    EX_EDGE_MAP   = 11 #< edge map property code     */
    EX_FACE_MAP   = 12 #< face map property code     */

    EX_GLOBAL     = 13 #< global "block" for variables*/
    EX_COORDINATE = 15 #< kluge so some internal wrapper functions work */
    EX_ASSEMBLY   = 16 #< assembly property code */
    EX_BLOB       = 17 #< blob property code */
    EX_INVALID    = -1
end


text:
"""
Runs Impute 5 on chromosome 20 over 4 chunks.
Each chunk imputes 20 million SNPs using 2 or 3 threads.
Each script runs a different chunk and are initialized simultaneously.
"""
function impute5_HRC_chr20(
    reffile::AbstractString,
    tgtfile::AbstractString,
    )
    chr = 20                # chromosome 20
    allale_start = 1        # first SNP is 60309
    allele_end = 70000000   # last SNP is 69960292
    chunk_size = 20000000   # each chunk imputes 20cM
    nchunks = 4             # chunk 1 imputes first 20cM region, chunk 2 imputes next 20cM ...etc
    cthreads = [3, 3, 2, 2] # chunk 1 uses 3 threads, chunk 2 uses 3 threads...etc

    # start 4 processes, each imputing different regions, using 2-3 threads
    Threads.@threads for i in 1:nchunks
        cur_start = (i - 1) * chunk_size + 1
        cur_end = i == nchunks ? allele_end : chunk_size * i
        cur_threads = cthreads[i]
        run(`nohup /usr/bin/time /home/biona001/impute5_v1.1.4/impute5_1.1.4_static
            --h $reffile
            --g $tgtfile
            --r $chr:$cur_start-$cur_end
            --o impute5.chr$chr.result.chunk$i.vcf.gz
            --threads $cur_threads`)
    end
end

"""
Runs Impute 5 on chromosome 10 over 14 chunks.
Each chunk imputes 10 million SNPs using 1 thread.
Up to 8 chunks can be ran simultaneously (10 concurrent threads throw memory error).

Note using 20cM chunks (as recommended) would result in 7 chunks,
"""
function impute5_HRC_chr10(
    reffile::AbstractString,
    tgtfile::AbstractString,
    )
    chr = 10                # chromosome 10
    allale_start = 1        # first SNP is 94316
    allele_end = 140000000  # last SNP is 135523061
    chunk_size = 10000000   # each chunk imputes 10cM (running 20cM chunks with 7 chunks throws Our Of Memory Error)
    nchunks = 14            # chunk 1 imputes first 10cM region, chunk 2 imputes next 10cM ...etc
    cthreads = ones(Int, nchunks)  # each chunk runs 1 threads

    # run each chunk with 1 thread
    Threads.@threads for i in 1:nchunks
        cur_start = (i - 1) * chunk_size + 1
        cur_end = i == nchunks ? allele_end : chunk_size * i
        cur_threads = cthreads[i]
        run(`/home/biona001/impute5_v1.1.4/impute5_1.1.4_static
            --h $reffile
            --g $tgtfile
            --r $chr:$cur_start-$cur_end
            --o impute5.chr$chr.result.chunk$i.vcf.gz
            --threads $cur_threads`)
    end
end

chr = parse(Int, ARGS[1])
reffile = ARGS[2]
tgtfile = ARGS[3]

if Threads.nthreads() == 1
    error("start Julia with >1 thread for this script to work!")
end

if chr == 20
    _, t, _, _, _ = @timed impute5_HRC_chr20(reffile, tgtfile)
elseif chr == 10
    _, t, _, _, _ = @timed impute5_HRC_chr10(reffile, tgtfile)
end

println("Total elapsed time = $t seconds")

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using Test
using BlockSystems
using ModelingToolkit
using ModelingToolkit: get_iv, get_eqs, get_states, value
using Graphs

@info "Tests of BlockSystems.jl"

@testset "BlockSystems.jl" begin
    @testset "namespaced accessors" begin
        @parameters t i1(t) i2(t) a
        @variables x1(t) x2(t) o1(t) o2(t)
        D = Differential(t)
        eqs  = [D(x1) ~ a*i1, o1~i1, D(x2) ~ i2, o2~i2]
        iob = IOBlock(eqs, [i1, i2], [o1, o2], name=:ns)
        @test Set(BlockSystems.namespace_inputs(iob)) == Set([iob.i1, iob.i2])
        @test Set(BlockSystems.namespace_outputs(iob)) == Set([iob.o1, iob.o2])
        @test Set(BlockSystems.namespace_istates(iob)) == Set([iob.x1, iob.x2])
        @test Set(BlockSystems.namespace_iparams(iob)) == Set([iob.a])
    end

    @testset "creation of IOBlocks" begin
        @parameters t i1(t) i2(t) a b
        @variables x1(t) x2(t) o1(t) o2(t)
        D = Differential(t)
        eqs = [D(x1) ~ a*i1,
               D(x2) ~ b*i2,
               o1 ~ a*x1,
               o2 ~ b*x2]

        iob = IOBlock(eqs, [i1, i2], [o1, o2], name=:iob)
        @test Set(iob.inputs) == Set([i1, i2])
        @test Set(iob.iparams) == Set([a, b])
        @test Set(iob.istates) == Set([x1, x2])
        @test Set(iob.outputs) == Set([o1, o2])
        @test Set(iob.removed_states) == Set()

        # x1 is a state not an input. Assign it as input will elad to warn + name conflict between inputs and states
        @test_throws ArgumentError IOBlock(eqs, [x1], [o1,o2], warn=false)

        @test_throws ArgumentError IOBlock(eqs, [i1,i2], [i1,o1,o2], warn=false)

        @parameters i a(t)
        @variables x(t) o(t)
        sys = ODESystem( [D(x) ~ a * i], name=:foo)
        aeq = [i ~ 2*a + i1]
        # name conflict betwen ios and odes
        @test_throws ArgumentError IOBlock(:name, [i.val], [a.val], [], [x.val], sys, aeq; warn=true)

        @parameters t a(t) p
        @variables x(t) y(t)
        D = Differential(t)
        IOBlock([x ~ 2 + a], [], [x])

        # implecit output, warn
        @test_logs (:warn, ) IOBlock([D(y) ~ 2 + a], [], [x])

        @test_throws ArgumentError IOBlock([a ~ 2 + x], [], [x])
        @test_throws ArgumentError IOBlock([p ~ 2 + x], [p], [x], iv=t)
        @test_throws ArgumentError IOBlock([a ~ 2 + x], [a], [x], iv=t)
        @test_logs (:warn, ) IOBlock([a ~ 2 + x + 2*a], [a], [x], iv=t)
    end

    @testset "test of create_namespace_promotions" begin
        using BlockSystems: create_namespace_promotions
        @parameters t
        Aa, Ba, Ab, Bc = value.(@parameters A₊a B₊a(t) A₊b B₊c)
        Ax, Bx, Ay, Bz = value.(@variables A₊x(t) B₊x(t) A₊y(t) B₊z(t))
        b, c = value.(@parameters b c)
        y, z = value.(@variables y(t) z(t))

        prom = create_namespace_promotions([Aa, Ba, Ab, Bc, Ax, Bx, Ay, Bz], [], true)
        @test Set(values(prom)) == Set([Aa, Ba, b, c, Ax, Bx, y, z])

        prom = create_namespace_promotions([Aa, Ba, Ab, Bc, Ax, Bx, Ay, Bz], [b, y], true)
        @test Set(values(prom)) == Set([Aa, Ba, Ab, c, Ax, Bx, Ay, z])

        prom = create_namespace_promotions([Aa, Ba, Ab, Bc, Ax, Bx, Ay, Bz], [], false)
        @test Set(values(prom)) == Set([Aa, Ba, Ab, Bc, Ax, Bx, Ay, Bz])

        prom = create_namespace_promotions([Aa, Ba, Ab, Bc, Ax, Bx, Ay, Bz], [b, y], false)
        @test Set(values(prom)) == Set([Aa, Ba, Ab, Bc, Ax, Bx, Ay, Bz])
    end

    @testset "IOBlock from other IOBlock" begin
        @parameters t i1(t) i2(t) a b
        @variables x1(t) x2(t) o1(t) o2(t)
        D = Differential(t)
        eqs = [D(x1) ~ a*i1,
               D(x2) ~ b*i2,
               o1 ~ a*x1,
               o2 ~ b*x2]
        iob1 = IOBlock(eqs, [i1, i2], [o1, o2], name=:iob1)
        iob2 = IOBlock(iob1, name=:iob2)
        iob3 = IOBlock(iob1)
        @test iob1.name != iob2.name != iob3.name
        @test Set(iob1.inputs) == Set(iob2.inputs) == Set(iob3.inputs)
        @test Set(iob1.iparams) == Set(iob2.iparams) == Set(iob3.iparams)
        @test Set(iob1.istates) == Set(iob2.istates) == Set(iob3.istates)
        @test Set(iob1.outputs) == Set(iob2.outputs) == Set(iob3.outputs)
        @test get_eqs(iob1.system) == get_eqs(iob2.system) == get_eqs(iob3.system)
        @test iob1.system.name == iob1.name == :iob1
        @test iob2.system.name == iob2.name == :iob2
        @test iob3.system.name == iob3.name
    end

    @testset "test creation of namespace map" begin
        @parameters t i(t) a b
        @variables x(t) o(t)
        D = Differential(t)
        eqs  = [D(x) ~ a*i, o ~ b*i]

        @parameters i2(t) a
        @variables x2(t) o(t)
        eqs2  = [D(x) ~ a*i2, D(x2) ~ i2, o~i2]

        iob1 = IOBlock(eqs, [i], [o], name=:iob1)
        iob2 = IOBlock(eqs2, [i2], [o], name=:iob2)

        @test Set(iob1.inputs) == Set(i)
        @test Set(iob2.inputs) == Set(i2)
        @test Set(iob1.outputs) == Set(o)
        @test Set(iob2.outputs) == Set(o)
        @test Set(iob1.istates) == Set(x)
        @test Set(iob2.istates) == Set([x, x2])
        @test Set(iob1.iparams) == Set([a, b])
        @test Set(iob2.iparams) == Set(a)

        sys = IOSystem([], [iob1, iob2])
        @test Set(sys.inputs) == Set([i, i2])
        @test Set(sys.istates) == Set([iob1.x, iob2.x, x2])
        @test Set(sys.iparams) == Set([iob1.a, iob2.a, b])
        @test Set(sys.outputs) == Set([iob1.o, iob2.o])
    end

    @testset "iosystem asserts" begin
        @parameters t i(t)
        @variables o(t)
        iob1 = IOBlock([o ~ i],[i],[o],name=:name)
        iob2 = IOBlock([o ~ i],[i],[o],name=:name)

        # namespace collision
        @test_throws ArgumentError IOSystem([iob1.o => iob2.i], [iob1, iob2])

        iob2 = IOBlock([o ~ i],[i],[o])
        # multiple connections to same input
        @test_throws ArgumentError IOSystem([iob1.o => iob1.i, iob2.o => iob1.i], [iob1, iob2])
        # make sure that alle of form input => output
        @test_throws ArgumentError IOSystem([iob1.o => iob1.o], [iob1, iob2])
        @test_throws ArgumentError IOSystem([iob1.i => iob1.i], [iob1, iob2])

        # assert that input maps refer to open inputs
        @test_throws ArgumentError IOSystem([iob1.o => iob2.i], [iob1, iob2],
                                            namespace_map = [iob2.i => i])
        # assert that rhs of input map is unique
        iob3 = IOBlock([o ~ i],[i],[o])
        @test_throws ArgumentError IOSystem([iob1.o => iob2.i],
                                             [iob1, iob2, iob3],
                                             namespace_map = [iob1.i => i, iob3.i => i])

        # test assertions for iparams and istats map
        @parameters t a i(t) b c
        @variables x(t) o(t) y(t)
        D = Differential(t)
        iob1 = IOBlock([D(x)~ i, o~a*x], [i], [o], name=:iob1)
        iob2 = IOBlock([D(x)~ i, o~a*x], [i], [o], name=:iob2)
        IOSystem([iob1.o => iob2.i], [iob1, iob2])
        # rhs unique
        @test_throws ArgumentError IOSystem([iob1.o => iob2.i], [iob1, iob2],
                                            namespace_map = [iob1.a=>b, iob2.a=>b])
        @test_throws ArgumentError IOSystem([iob1.o => iob2.i], [iob1, iob2],
                                            namespace_map = [iob1.x=>y, iob2.x=>y])

        # rhs unique
        @test_throws ArgumentError IOSystem([iob1.o => iob2.i], [iob1, iob2],
                                            namespace_map = [iob1.o=>y, iob2.o=>y])
    end

    function test_complete_namespace_promotions(ios)
        eqs = vcat([ModelingToolkit.namespace_equations(iob.system) for iob in ios.systems]...)
        allvars = [(get_variables(eq.lhs) ∪ get_variables(eq.rhs)) for eq in eqs]
        allvars = union(allvars...) |> unique
        allvars = setdiff(allvars, [BlockSystems.get_iv(ios)])
        allkeys = keys(ios.namespace_map)
        # closed inputs should not appear in allkeys
        @test isempty(Set(allkeys) ∩ Set(last.(ios.connections)))
        # but there should be no namespace collision with them either..
        allkeys = Set(allkeys) ∪ Set(last.(ios.connections))
        @test allunique(allkeys)
        @test Set(allkeys) == Set(allvars)
    end

    @testset "test creation of systems" begin
        #=
                    +------------+
        in1 --> i1 -|    iob1    |
        in2 --> i2 -|(x1, x2)(a) |-o--+     +-----+
                    +------------+    +-ina-| add |- add ---> out
                    +------------+    +-inb-|     |
        in3 --> i1 -|    iob2    |-o--+     +-----+
        in4 --> i2 -|(x1, x2)(b) |
                    +------------+
        =#
        @parameters t i1(t) i2(t) a b ina(t) inb(t)
        @variables x1(t) x2(t) o(t) add(t)
        D = Differential(t)
        eqs1  = [D(x1) ~ a*i1, D(x2)~i2, o~x1+x2]
        iob1 = IOBlock(eqs1, [i1, i2], [o], name=:iob1)

        eqs2  = [D(x1) ~ b*i1, D(x2)~i2, o~x1+x2]
        iob2 = IOBlock(eqs2, [i1, i2], [o], name=:iob2)

        ioadd = IOBlock([add ~ ina + inb], [ina, inb], [add], name=:add)

        # try with auto namespacing
        sys = IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                       [iob1, iob2, ioadd],
                       name=:sys)
        @test Set(sys.inputs) == Set([iob1.i1, iob1.i2, iob2.i1, iob2.i2])
        @test Set(sys.iparams) == Set([a, b])
        @test Set(sys.istates) == Set([iob1.x1, iob1.x2, iob2.x1, iob2.x2])
        @test Set(sys.outputs) == Set([iob1.o, iob2.o, add])
        test_complete_namespace_promotions(sys)

        # provide maps
        @parameters in1(t) in2(t) in3(t) in4(t) p1 p2
        @variables out(t) y1(t) y2(t)
        sys1 = IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                       [iob1, iob2, ioadd],
                       namespace_map = Dict(iob1.i1 => in1,
                                            iob1.i2 => in2,
                                            iob2.i1 => in3,
                                            iob2.i2 => in4,
                                            iob1.a => p1,
                                            iob2.b => p2,
                                            iob1.x1 => y1,
                                            iob1.x2 => y2,
                                            iob2.x1 => x1,
                                            iob2.x2 => x2,
                                            ioadd.add => out),
                       outputs = [ioadd.add],
                       name=:sys)
        sys2 = IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                       [iob1, iob2, ioadd],
                       namespace_map = Dict(iob1.i1 => :in1,
                                            iob1.i2 => :in2,
                                            iob2.i1 => :in3,
                                            iob2.i2 => :in4,
                                            iob1.a => :p1,
                                            iob2.b => :p2,
                                            iob1.x1 => :y1,
                                            iob1.x2 => :y2,
                                            iob2.x1 => :x1,
                                            iob2.x2 => :x2,
                                            ioadd.add => :out),
                       outputs = [:out],
                       name=:sys)
        @test Set(sys1.inputs) == Set([in1, in2, in3, in4])
        @test Set(sys1.iparams) == Set([p1, p2])
        @test Set(sys1.istates) == Set([y1, y2, x1, x2, iob1.o, iob2.o])
        @test Set(sys1.outputs) == Set([out])
        test_complete_namespace_promotions(sys1)

        @test Set(sys2.inputs) == Set([in1, in2, in3, in4])
        @test Set(sys2.iparams) == Set([p1, p2])
        @test Set(sys2.istates) == Set([y1, y2, x1, x2, iob1.o, iob2.o])
        @test Set(sys2.outputs) == Set([out])
        test_complete_namespace_promotions(sys2)

        # provide partial maps
        sys = IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                       [iob1, iob2, ioadd],
                       namespace_map = Dict(iob1.i1 => in1,
                                            iob1.i2 => in2,
                                            iob1.a => p1,
                                            iob1.x1 => y1,
                                            iob1.x2 => y2,
                                            ioadd.add => out),
                       outputs = [out], # provide outputs as namespaced variable
                       name=:sys)
        @test Set(sys.inputs) == Set([in1, in2, i1, i2])
        @test Set(sys.iparams) == Set([p1, b])
        @test Set(sys.istates) == Set([y1, y2, x1, x2, iob1.o, iob2.o])
        @test Set(sys.outputs) == Set([out])
        test_complete_namespace_promotions(sys)

        # check for argument error for bad outputs argument
        sys1 = IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                       [iob1, iob2, ioadd],
                       namespace_map = Dict(ioadd.add => out),
                       outputs = [ioadd.add], # provide outputs as namespaced variable
                       name=:sys)
        sys2 = IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                       [iob1, iob2, ioadd],
                       namespace_map = Dict(ioadd.add => out),
                       outputs = [out], # provide outputs as namespaced variable
                       name=:sys)
        @test Set(sys1.istates) == Set(sys2.istates)
        @test Set(sys1.outputs) == Set(sys2.outputs)
        @test_throws ArgumentError IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                                            [iob1, iob2, ioadd],
                                            namespace_map = Dict(ioadd.add => out),
                                            outputs = [p1], # provide outputs as namespaced variable
                                            name=:sys)
    end

    @testset "algebraic connections" begin
        @parameters t p1 p2 i(t)
        @variables o1(t) o2(t) o3(t)
        dt = Differential(t)

        @named b1 = IOBlock([o1 ~ p1], [], [o1])
        @named b2 = IOBlock([o2 ~ p2], [], [o2])
        @named b3 = IOBlock([dt(o3) ~ i], [i], [o3])

        sys = IOSystem([b1.o1 + b2.o2 => b3.i],
                       [b1, b2, b3], outputs=[b3.o3])
        con = connect_system(sys, verbose=false)
        @test equations(con) == [dt(o3) ~ p1 + p2]
    end

    @testset "test BlockSpec" begin
        # test the constructors
        @parameters t
        @parameters uᵢ(t) uᵣ(t)
        @variables x(t) iᵢ(t) iᵣ(t)
        bs1 = BlockSpec([:uᵣ, :uᵢ], [:iᵣ, :iᵢ])
        bs2 = BlockSpec(value.([uᵣ, uᵢ]), value.([iᵣ, iᵢ]))
        bs3 = BlockSpec([uᵣ, uᵢ], [iᵣ, iᵢ])
        @test bs1.inputs == bs2.inputs == bs3.inputs
        @test bs1.outputs == bs2.outputs == bs3.outputs

        D = Differential(t)
        eqs  = [D(x) ~ x, iᵢ ~ uᵢ, iᵣ ~ uᵣ]
        iob1 = IOBlock(eqs, [uᵢ, uᵣ], [iᵢ, iᵣ])
        iob2 = IOBlock(eqs, [uᵢ], [iᵢ, iᵣ])
        iob3 = IOBlock(eqs, [uᵢ, uᵣ], [iᵢ])

        @test fulfills(iob1, bs1) == bs1(iob1) == true
        @test fulfills(iob2, bs1) == bs1(iob2) == false
        @test fulfills(iob3, bs1) == bs1(iob3) == false

        sys1 = IOSystem([], [iob1, iob2])
        sys2 = IOSystem([], [iob1, iob2], outputs=[iᵢ, iᵣ],
                        namespace_map = [iob1.iᵢ => iᵢ,
                                         iob1.iᵣ => iᵣ,
                                         iob1.uᵢ => uᵢ,
                                         iob1.uᵣ => uᵣ])

        @test BlockSpec([:uᵣ, :uᵢ], [:iᵣ, :iᵢ])(sys1) == false # no name promotion
        @test BlockSpec([:uᵣ, :uᵢ], [:iᵣ, :iᵢ])(sys2) == false # to much inputs

        # test the strict stuff
        @parameters t
        @parameters i1(t) i2(t)
        @variables o1(t) o2(t)
        iob = IOBlock([o1 ~ i1, o2 ~ i2], [i1, i2], [o1, o2], iv=t)

        @test BlockSpec([:i1, :i2], [:o1])(iob) == true
        @test BlockSpec([:i1, :i2], [:o1], in_strict=1, out_strict=1)(iob) == false
        @test BlockSpec([:i1], [:o1], in_strict=1, out_strict=1)(iob) == false
        @test BlockSpec([:i1], [:o1], in_strict=0, out_strict=0)(iob) == true

        @test BlockSpec([:i1], [:o1, :o2])(iob) == false
        @test BlockSpec([:i1], [:o1, :o2], in_strict=0, out_strict=0)(iob) == true
        @test BlockSpec([:i1], [:o1, :o2], in_strict=0, out_strict=1)(iob) == true
    end
end


text:
# NOTE:
# tests in this file may not be so robust against the future changes in JET.jl or
# even those in julia itself

@testset "print toplevel errors" begin
    let
        io = IOBuffer()
        s = """
            a = begin
                b =
            end
            """

        res = analyze_text(s, @__FILE__)
        print_reports(io, res.toplevel_error_reports)
        let s = String(take!(io))
            @test occursin("1 toplevel error found", s)
            @test occursin(Regex("@ $(@__FILE__):\\d"), s)
            @test occursin("syntax: unexpected \"end\"", s)
        end

        res = analyze_text(s, "foo")
        print_reports(io, res.toplevel_error_reports)
        let s = String(take!(io))
            @test occursin("1 toplevel error found", s)
            @test occursin(r"@ foo:\d", s)
            @test occursin("syntax: unexpected \"end\"", s)
        end
    end
end

@testset "print inference errors" begin
    let
        res = @analyze_toplevel begin
            s = "julia"
            sum(s)
        end

        io = IOBuffer()
        @test print_reports(io, res.inference_error_reports)
        let s = String(take!(io))
            @test occursin("2 possible errors found", s)
            @test occursin(Regex("@ $(escape_string(@__FILE__)):$((@__LINE__)-7)"), s) # toplevel call signature
        end
    end

    @testset "special case splat call signature" begin
        let
            res = @analyze_toplevel begin
                foo(args...) = sum(args)
                foo(rand(Char, 1000000000)...)
            end

            io = IOBuffer()
            @test print_reports(io, res.inference_error_reports, JET.gen_postprocess(res.actual2virtual))
            let s = String(take!(io))
                @test occursin("1 possible error found", s)
                @test occursin(Regex("@ $(escape_string(@__FILE__)):$((@__LINE__)-8)"), s) # toplevel call signature
                @test occursin("foo(rand(Char, 1000000000)...)", s)
            end
        end
    end
end


text:
# TestNet: https://en.bitcoin.it/wiki/Testnet
# - Faucet: http://faucet.xeno-genesis.com/
# - Explorer: http://tbtc.blockr.io/
# - Peerlist: nslookup testnet-seed.bitcoin.petertodd.org

Crypto.init()

const LISTEN_PORT = 18333 # 8333 for mainnet
const RPC_PORT = 18332 # 8332 for mainnet

priv_key = "7e68e472a6b41f165c0d13b57b6ac88440367c25361b5d6d407bb5fe2cd05a12"
# moZRBjoq3ELstSp6CPeBXir5FSpjjHyQMp
pub_key = get_pub_key(priv_key, network_id = 0x6f)

connect(ip"107.170.32.58", 18333)


text:
"""
    abstract type RotationGenerator{N,T} <: StaticMatrix{N,N,T}

An abstract type representing `N`-dimensional rotation generator. More abstractly, they represent
skew-symmetric real `N`×`N` matrices.
"""
abstract type RotationGenerator{N,T} <: StaticMatrix{N,N,T} end

Base.@pure StaticArrays.Size(::Type{RotationGenerator{N}}) where {N} = Size(N,N)
Base.@pure StaticArrays.Size(::Type{RotationGenerator{N,T}}) where {N,T} = Size(N,N)
Base.@pure StaticArrays.Size(::Type{R}) where {R<:RotationGenerator} = Size(supertype(R))
Base.adjoint(r::RotationGenerator) = -r
Base.transpose(r::RotationGenerator{N,T}) where {N,T<:Real} = -r
Base.real(R::Type{<:RotationGenerator}) = R

# Generate identity-matrix with SMatrix
# Note that ones(RotationGenerator3,dims...) is not Array{<:RotationGenerator} but Array{<:StaticMatrix{3,3}}
Base.one(::RotationGenerator{N,T}) where {N,T} = one(SMatrix{N,N,T})
Base.one(::Type{RotationGenerator}) = error("The dimension of rotation is not specified.")
Base.one(::Type{<:RotationGenerator{N}}) where N = one(SMatrix{N,N})
Base.one(::Type{<:RotationGenerator{N,T}}) where {N,T} = one(SMatrix{N,N,T})
Base.ones(::Type{R}) where {R<:RotationGenerator} = ones(R, ()) # avoid StaticArray constructor
Base.ones(::Type{R}, dims::Base.DimOrInd...) where {R<:RotationGenerator} = ones(typeof(one(R)),dims...)
Base.ones(::Type{R}, dims::NTuple{N, Integer}) where {R<:RotationGenerator, N} = ones(typeof(one(R)),dims)

# Generate zero rotation matrix
Base.zero(r::RotationGenerator) = zero(typeof(r))

# difference
Base.:-(r1::RotationGenerator, r2::RotationGenerator) = r1 + (-r2)

# `convert` goes through the constructors, similar to e.g. `Number`
Base.convert(::Type{R}, rot::RotationGenerator{N}) where {N,R<:RotationGenerator{N}} = R(rot)

# RotationGenerator matrices should be skew-symmetric. Only the operations we define,
# like addition, will stay as RotationGenerators, otherwise users will get an
# SMatrix{3,3} (e.g. rot1 * rot2 -> SMatrix)
Base.@pure StaticArrays.similar_type(::Union{R,Type{R}}) where {R <: RotationGenerator} = SMatrix{size(R)..., eltype(R), prod(size(R))}
Base.@pure StaticArrays.similar_type(::Union{R,Type{R}}, ::Type{T}) where {R <: RotationGenerator, T} = SMatrix{size(R)..., T, prod(size(R))}

################################################################################
################################################################################
"""
    struct RotMatrixGenerator{N,T} <: RotationGenerator{N,T}

A statically-sized, N×N skew-symmetric matrix.

Note: the skew-symmetricity of the input matrix is *not* checked by the constructor.
"""
struct RotMatrixGenerator{N,T,L} <: RotationGenerator{N,T} # which is <: AbstractMatrix{T}
    mat::SMatrix{N, N, T, L} # The final parameter to SMatrix is the "length" of the matrix, 3 × 3 = 9
    RotMatrixGenerator{N,T,L}(x::AbstractArray) where {N,T,L} = new{N,T,L}(convert(SMatrix{N,N,T,L}, x))
    # fixes #49 ambiguity introduced in StaticArrays 0.6.5
    RotMatrixGenerator{N,T,L}(x::StaticArray) where {N,T,L} = new{N,T,L}(convert(SMatrix{N,N,T,L}, x))
end
RotMatrixGenerator(x::SMatrix{N,N,T,L}) where {N,T,L} = RotMatrixGenerator{N,T,L}(x)

Base.zero(::Type{RotMatrixGenerator}) = error("The dimension of rotation is not specified.")

# These functions (plus size) are enough to satisfy the entire StaticArrays interface:
for N in 2:3
    L = N*N
    RotMatrixGeneratorN = Symbol(:RotMatrixGenerator, N)
    @eval begin
        @inline RotMatrixGenerator(t::NTuple{$L})  = RotMatrixGenerator(SMatrix{$N,$N}(t))
        @inline (::Type{RotMatrixGenerator{$N}})(t::NTuple{$L}) = RotMatrixGenerator(SMatrix{$N,$N}(t))
        @inline RotMatrixGenerator{$N,T}(t::NTuple{$L}) where {T} = RotMatrixGenerator(SMatrix{$N,$N,T}(t))
        @inline RotMatrixGenerator{$N,T,$L}(t::NTuple{$L}) where {T} = RotMatrixGenerator(SMatrix{$N,$N,T}(t))
        const $RotMatrixGeneratorN{T} = RotMatrixGenerator{$N, T, $L}
    end
end
Base.@propagate_inbounds Base.getindex(r::RotMatrixGenerator, i::Int) = r.mat[i]
@inline Base.Tuple(r::RotMatrixGenerator) = Tuple(r.mat)

@inline RotMatrixGenerator(v::Number) = RotMatrixGenerator{2}(v)
@inline function (::Type{RotMatrixGenerator{2}})(v::Number)
    RotMatrixGenerator(@SMatrix [zero(v) -v; v zero(v)])
end
@inline function RotMatrixGenerator{2,T}(v::Number) where T
    RotMatrixGenerator(@SMatrix T[zero(v) -v; v zero(v)])
end

# A rotation is more-or-less defined as being an orthogonal (or unitary) matrix
Base.:-(r::RotMatrixGenerator) = RotMatrixGenerator(r.mat')

# By default, composition of rotations will go through RotMatrixGenerator, unless overridden
@inline Base.:+(r1::RotationGenerator, r2::RotationGenerator) = RotMatrixGenerator(r1) + RotMatrixGenerator(r2)
@inline Base.:+(r1::RotMatrixGenerator, r2::RotationGenerator) = r1 + RotMatrixGenerator(r2)
@inline Base.:+(r1::RotationGenerator, r2::RotMatrixGenerator) = RotMatrixGenerator(r1) + r2
@inline Base.:+(r1::RotMatrixGenerator, r2::RotMatrixGenerator) = RotMatrixGenerator(r1.mat + r2.mat)

@inline Base.:*(t::Number, r::RotMatrixGenerator) = RotMatrixGenerator(t*r.mat)
@inline Base.:*(r::RotMatrixGenerator, t::Number) = t*r
@inline Base.:/(r::RotMatrixGenerator, t::Number) = RotMatrixGenerator(r.mat/t)

# Special case multiplication of 2×2 rotation matrices: speedup using skew-symmetricity
@inline function Base.:+(r1::RotMatrixGenerator{2}, r2::RotMatrixGenerator{2})
    v = r1[2,1]+r2[2,1]
    s = @SMatrix [0  -v
                  v   0]
    return RotMatrixGenerator(s)
end

# Special case multiplication of 3×3 rotation matrices: speedup using skew-symmetricity
@inline function Base.:+(r1::RotMatrixGenerator{3}, r2::RotMatrixGenerator{3})
    x = r1[6]+r2[6]
    y = r1[7]+r2[7]
    z = r1[2]+r2[2]
    s = @SMatrix [ 0 -z  y
                   z  0 -x
                  -y  x  0]

    return RotMatrixGenerator(s)
end

"""
    struct Angle2dGenerator{T} <: RotationGenerator{2,T}
        v::T
    end

A 2×2 rotation generator matrix (i.e. skew-symmetric matrix).
[ 0 -v
  v  0 ]
"""
struct Angle2dGenerator{T} <: RotationGenerator{2,T}
    v::T
    Angle2dGenerator{T}(r) where T = new{T}(r)
end

@inline function Angle2dGenerator(r::T) where T <: Number
    Angle2dGenerator{T}(r)
end

@inline Angle2dGenerator(r::RotationGenerator{2}) = Angle2dGenerator(r[2,1])
@inline Angle2dGenerator{T}(r::RotationGenerator{2}) where {T} = Angle2dGenerator{T}(r[2,1])

@inline Base.zero(::Type{A}) where {A<: Angle2dGenerator} = A(0.0)

@inline Base.:+(r1::Angle2dGenerator, r2::Angle2dGenerator) = Angle2dGenerator(r1.v + r2.v)
@inline Base.:-(r1::Angle2dGenerator, r2::Angle2dGenerator) = Angle2dGenerator(r1.v - r2.v)
@inline Base.:*(t::Number, r::Angle2dGenerator) = Angle2dGenerator(t*r.v)
@inline Base.:*(r::Angle2dGenerator, t::Number) = t*r
@inline Base.:/(r::Angle2dGenerator, t::Number) = Angle2dGenerator(r.v/t)
@inline Base.:-(r::Angle2dGenerator) = Angle2dGenerator(-r.v)

@inline function Base.getindex(r::Angle2dGenerator{T}, i::Int) where T
    if i == 1
        zero(T)
    elseif i == 2
        r.v
    elseif i == 3
        -r.v
    elseif i == 4
        zero(T)
    else
        throw(BoundsError(r,i))
    end
end

@inline params(r::Angle2dGenerator) = SVector{1}(r.v)

"""
    struct RotationVecGenerator{T} <: RotationGenerator{2,T}
        x::T
        y::T
        z::T
    end

A 3×3 rotation generator matrix (i.e. skew-symmetric matrix).
[ 0 -z  y
  z  0 -x
 -y  x  0]
"""
struct RotationVecGenerator{T} <: RotationGenerator{3,T}
    x::T
    y::T
    z::T
    RotationVecGenerator{T}(x,y,z) where T = new{T}(x,y,z)
end

@inline function RotationVecGenerator(x::X,y::Y,z::Z) where {X<:Number, Y<:Number, Z<:Number}
    RotationVecGenerator{promote_type(X,Y,Z)}(x, y, z)
end

@inline RotationVecGenerator(r::RotationGenerator{3}) = RotationVecGenerator(r[6], r[7], r[2])
@inline RotationVecGenerator{T}(r::RotationGenerator{3}) where {T} = RotationVecGenerator{T}(r[6], r[7], r[2])

@inline Base.zero(::Type{R}) where {R<: RotationVecGenerator} = R(0.0,0.0,0.0)

@inline Base.:+(r1::RotationVecGenerator, r2::RotationVecGenerator) = RotationVecGenerator(r1.x + r2.x, r1.y + r2.y, r1.z + r2.z)
@inline Base.:-(r1::RotationVecGenerator, r2::RotationVecGenerator) = RotationVecGenerator(r1.x - r2.x, r1.y - r2.y, r1.z - r2.z)
@inline Base.:*(t::Number, r::RotationVecGenerator) = RotationVecGenerator(t*r.x, t*r.y, t*r.z)
@inline Base.:*(r::RotationVecGenerator, t::Number) = t*r
@inline Base.:/(r::RotationVecGenerator, t::Number) = RotationVecGenerator(r.x/t, r.y/t, r.z/t)
@inline Base.:-(r::RotationVecGenerator) = RotationVecGenerator(-r.x, -r.y, -r.z)

@inline function Base.getindex(r::RotationVecGenerator{T}, i::Int) where T
    if i == 1
        zero(T)
    elseif i == 2
        r.z
    elseif i == 3
        -r.y
    elseif i == 4
        -r.z
    elseif i == 5
        zero(T)
    elseif i == 6
        r.x
    elseif i == 7
        r.y
    elseif i == 8
        -r.x
    elseif i == 9
        zero(T)
    else
        throw(BoundsError(r,i))
    end
end

@inline params(r::RotationVecGenerator) = SVector{3}(r.x,r.y,r.z)

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

# A simplification and specialization of the Base.show function for AbstractArray makes
# everything sensible at the REPL.
function Base.show(io::IO, ::MIME"text/plain", X::RotationGenerator)
    if !haskey(io, :compact)
        io = IOContext(io, :compact => true)
    end
    summary(io, X)
    if !isa(X, RotMatrixGenerator)
        n_fields = length(fieldnames(typeof(X)))
        print(io, "(")
        for i in 1:n_fields
            print(io, getfield(X, i))
            if i < n_fields
                print(io, ", ")
            end
        end
        print(io, ")")
    end
    print(io, ":")
    println(io)
    io = IOContext(io, :typeinfo => eltype(X))
    Base.print_array(io, X)
end

"""
    isrotationgenerator(r)

Check whether `r` is a rotation generator matrix (skew-symmetric matrix).
"""
isrotationgenerator

function isrotationgenerator(r::StaticMatrix{N,N,T}) where {N,T<:Real}
    # This method can be removed, but is a little faster than isrotationgenerator(r::AbstractMatrix{T}).
    m = SMatrix(r)
    return m == -m'
end

function isrotationgenerator(r::AbstractMatrix{T}) where T<:Real
    indsm, indsn = axes(r)
    if indsm != indsn
        return false
    end
    for i in indsn, j in i:last(indsn)
        if r[i,j] != -r[j,i]
            return false
        end
    end
    return true
end

function isrotationgenerator(r::AbstractMatrix{T}) where T<:Number
    if !isreal(r)
        return false
    end
    return isrotationgenerator(real(r))
end


text:
using AxisKeys
using PrettyTables
export PsisLoo, PsisLooMethod, Psis, BayesCV

const POINTWISE_LABELS = (:cv_est, :naive_est, :overfit, :ess, :pareto_k)
const CV_DESC = """
# Fields

  - `estimates::KeyedArray`: A KeyedArray with columns `:total, :se_total, :mean, :se_mean`,
    and rows `:cv_est, :naive_est, :overfit`. See `# Extended help` for more.
      - `:cv_est` contains estimates for the out-of-sample prediction error, as
        predicted using the jackknife (LOO-CV).
      - `:naive_est` contains estimates of the in-sample prediction error.
      - `:overfit` is the difference between the previous two estimators, and estimates
        the amount of overfitting. When using the log probability score, it is equal to
        the effective number of parameters -- a model with an overfit of 2 is "about as
        overfit" as a model with 2 independent parameters that have a flat prior.
  - `pointwise::KeyedArray`: A `KeyedArray` of pointwise estimates with 5 columns --
      - `:cv_est` contains the estimated out-of-sample error for this point, as measured
        using leave-one-out cross validation.
      - `:naive_est` contains the in-sample estimate of error for this point.
      - `:overfit` is the difference in the two previous estimates.
      - `:ess` is the effective sample size, which measures the simulation error caused by
        using Monte Carlo estimates. It is *not* related to the actual sample size, and it
        does not measure how accurate your predictions are.
    - `:pareto_k` is the estimated value for the parameter `ξ` of the generalized Pareto
      distribution. Values above .7 indicate that PSIS has failed to approximate the true
      distribution.
  - `psis_object::Psis`: A `Psis` object containing the results of Pareto-smoothed
    importance sampling.

# Extended help

The total score depends on the sample size, and summarizes the weight of evidence for or
against a model. Total scores are on an interval scale, meaning that only differences of
scores are meaningful. *It is not possible to interpret a total score by looking at it.*
The total score is not a relative goodness-of-fit statistic (for this, see the average
score).

The overfit is equal to the difference between the in-sample and out-of-sample predictive
accuracy. When using the log probability score, it is equal to the "effective number of
parameters" -- a model with an overfit of 2 is "about as overfit" as a model with 2
free parameters and flat priors.

The average score is the total score, divided by the sample size. It estimates the expected
log score, i.e. the expectation of the log probability density of observing the next point.
The average score is a relative goodness-of-fit statistic which does not depend on sample
size.

Unlike for chi-square goodness of fit tests, models do not have to be nested for model
comparison using cross-validation methods.
"""

###########################
### IMPORTANCE SAMPLING ###
###########################

"""
    Psis{V<:AbstractVector{F},I<:Integer} where {F<:Real}

A struct containing the results of Pareto-smoothed importance sampling.

# Fields
  - `weights`: A vector of smoothed, truncated, and normalized importance sampling weights.
  - `pareto_k`: Estimates of the shape parameter `k` of the generalized Pareto distribution.
  - `ess`: Estimated effective sample size for each LOO evaluation.
  - `tail_len`: Vector indicating how large the "tail" is for each observation.
  - `dims`: Named tuple of length 2 containing `s` (posterior sample size) and `n` (number
    of observations).
"""
struct Psis{
    F<:Real,
    AF<:AbstractArray{F,3},
    VF<:AbstractVector{F},
    I<:Integer,
    VI<:AbstractVector{I},
}
    weights::AF
    pareto_k::VF
    ess::VF
    r_eff::VF
    tail_len::VI
    posterior_sample_size::I
    data_size::I
end

function _throw_pareto_k_warning(ξ)
    if any(ξ .≥ .7)
        @warn "Some Pareto k values are very high (>0.7), indicating that PSIS has " *
        "failed to approximate the true distribution."
    elseif any(ξ .≥ .5)
        @info "Some Pareto k values are slightly high (>0.5); some pointwise estimates " *
        "may be slow to converge or have high variance."
    end
end

function Base.show(io::IO, ::MIME"text/plain", psis_object::Psis)
    table = hcat(psis_object.pareto_k, psis_object.ess)
    post_samples = psis_object.posterior_sample_size
    data_size = psis_object.data_size
    println("Results of PSIS with $post_samples Monte Carlo samples and " *
    "$data_size data points.")
    _throw_pareto_k_warning(psis_object.pareto_k)
    return pretty_table(
        table;
        compact_printing=false,
        header=[:pareto_k, :ess],
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end

##########################
#### CROSS VALIDATION ####
##########################

"""
    AbstractCV
An abstract type used in cross-validation.
"""
abstract type AbstractCV end

"""
    AbstractCVMethod
An abstract type used to dispatch the correct method for cross validation.
"""
abstract type AbstractCVMethod end

##########################
######## PSIS-LOO ########
##########################

"""
    PsisLooMethod

Use Pareto-smoothed importance sampling together with leave-one-out cross validation to
estimate the out-of-sample predictive accuracy.
"""
struct PsisLooMethod <: AbstractCVMethod end

"""
    PsisLoo{
        F <: Real,
        AF <: AbstractArray{F},
        VF <: AbstractVector{F},
        I <: Integer,
        VI <: AbstractVector{I},
    } <: AbstractCV

A struct containing the results of jackknife (leave-one-out) cross validation using Pareto
smoothed importance sampling.

$CV_DESC

See also: [`loo`]@ref, [`bayes_cv`]@ref, [`psis_loo`]@ref, [`Psis`]@ref
"""
struct PsisLoo{
    F <: Real,
    AF <: AbstractArray{F},
    VF <: AbstractVector{F},
    I <: Integer,
    VI <: AbstractVector{I},
} <: AbstractCV
    estimates::KeyedArray
    pointwise::KeyedArray
    psis_object::Psis{F, AF, VF, I, VI}
end

function Base.show(io::IO, ::MIME"text/plain", loo_object::PsisLoo)
    table = loo_object.estimates
    _throw_pareto_k_warning(loo_object.pointwise(:pareto_k))
    post_samples = loo_object.psis_object.posterior_sample_size
    data_size = loo_object.psis_object.data_size
    println("Results of PSIS-LOO-CV with $post_samples Monte Carlo samples and " *
    "$data_size data points.")
    return pretty_table(
        table;
        compact_printing=false,
        header=table.statistic,
        row_names=table.criterion,
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end

##########################
### BAYESIAN BOOTSTRAP ###
##########################

"""
    BayesCV{
        F <: Real,
        AF <: AbstractArray{F},
        VF <: AbstractVector{F},
        I <: Integer,
        VI <: AbstractVector{I},
    } <: AbstractCV

A struct containing the results of cross-validation using the Bayesian bootstrap.

$CV_DESC

See also: [`bayes_cv`]@ref, [`psis_loo`]@ref, [`psis`]@ref, [`Psis`]@ref
"""
struct BayesCV{
    F <: Real,
    AF <: AbstractArray{F},
    VF <: AbstractVector{F},
    I <: Integer,
    VI <: AbstractVector{I},
} <: AbstractCV
    estimates::KeyedArray
    posteriors::KeyedArray
    psis_object::Psis{F, AF, VF, I, VI}
end

function Base.show(io::IO, ::MIME"text/plain", cv_object::BayesCV)
    table = cv_object.estimates
    post_samples = cv_object.psis_object.posterior_sample_size
    data_size = cv_object.psis_object.data_size
    _throw_pareto_k_warning(cv_object.resamples(:pareto_k))
    println("Results of Bayesian bootstrap CV with $post_samples Monte Carlo samples and " *
    "$data_size data points.")
    return pretty_table(
        table;
        compact_printing=false,
        header=table.statistic,
        row_names=table.criterion,
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end


text:
function excess_gibbs_free_energy(model::ActivityModel,p,T,z)
    γ = activity_coefficient(model,p,T,z)
    return sum(z[i]*R̄*T*log(γ[i]) for i ∈ @comps)
end

function eos(model::ActivityModel,V,T,z)
    Σz = sum(z)
    lnΣz = log(Σz)
    g_E = excess_gibbs_free_energy(model,V,T,z)
    pures = model.puremodel
    p = sum(z[i]*pressure(pures[i],V,T) for i ∈ @comps)/Σz
    g_ideal = sum(z[i]*R̄*T*(log(z[i])-lnΣz) for i ∈ @comps)
    g_pure = sum(z[i]*VT_gibbs_free_energy(pures[i],V,T) for i ∈ @comps)
    return g_E+g_ideal+g_pure-p*V
end

function eos_res(model::ActivityModel,V,T,z)
    Σz = sum(z)
    g_E = excess_gibbs_free_energy(model,V,T,z)
    pures = model.puremodel
    g_pure_res = sum(VT_chemical_potential_res(pures[i],V,T)[1]*z[i] for i ∈ @comps)
    p = sum(z[i]*pressure(pures[i],V,T) for i ∈ @comps)/Σz
    p_res = p - Σz*R̄*T/V
    return g_E+g_pure_res-p_res*V
end

function bubble_pressure(model::ActivityModel,T,x)
    sat = saturation_pressure.(model.puremodel,T)
    p_sat = [tup[1] for tup in sat]
    γ     = activity_coefficient(model,1e-4,T,x)
    p     = sum(x.*γ.*p_sat)
    y     = x.*γ.*p_sat ./ p
    return (p,y)
end

function bubble_temperature(model::ActivityModel,p,x;T0=nothing)
    f(z) = Obj_bubble_temperature(model,z,p,x)
    if T0===nothing
        pure = model.puremodel
        sat = saturation_temperature.(pure,p)
        Ti   = zero(x)
        for i ∈ 1:length(x)
            if isnan(sat[i][1])
                Tc,pc,vc = crit_pure(pure[i])
                g(x) = p-pressure(pure[i],vc,x,[1.])
                Ti[i] = Roots.find_zero(g,(Tc))
            else
                Ti[i] = sat[i][1]
            end
        end
        T = Roots.find_zero(f,(minimum(Ti)*0.9,maximum(Ti)*1.1))
    else
        T = Roots.find_zero(f,T0)
    end
    p,y = bubble_pressure(model,T,x)
    return (T,y)
end

function Obj_bubble_temperature(model::ActivityModel,T,p,x)
    sat = saturation_pressure.(model.puremodel,T)
    p_sat = [tup[1] for tup in sat]
    γ     = activity_coefficient(model,1e-4,T,x)
    y     = x.*γ.*p_sat ./ p
    return sum(y)-1
end

function mixing(model::ActivityModel,p,T,z,::typeof(enthalpy))
    f(x) = excess_gibbs_free_energy(model,p,x,z)/x
    df(x) = Solvers.derivative(f,x)
    return -df(T)*T^2
end

function mixing(model::ActivityModel,p,T,z,::typeof(gibbs_free_energy))
    x = z./sum(z)
    return excess_gibbs_free_energy(model,p,T,z)+dot(z,log.(x))*R̄*T
end

function mixing(model::ActivityModel,p,T,z,::typeof(entropy))
    f(x) = excess_gibbs_free_energy(model,p,x,z)/x
    g,dg = Solvers.f∂f(f,T)
    return -dg*T-g
end

function lb_volume(model::ActivityModel,z = SA[1.0])
    b = maximum([model.puremodel[i].params.b.values[1,1] for i ∈ @comps])
    return b
end

function T_scale(model::ActivityModel,z=SA[1.0])
    n = sum(z)
    invn2 = one(n)/(n*n)
    Ωa,Ωb = ab_consts(model.puremodel[1])
    _a = [model.puremodel[i].params.a.values[1,1] for i in @comps]
    _b = [model.puremodel[i].params.b.values[1,1] for i in @comps]
    a = dot(z, _a)*invn2/Ωa
    b = dot(z, _b)*invn2/Ωb
    return a/b/R̄
end

function p_scale(model::ActivityModel,z=SA[1.0])
    n = sum(z)
    invn2 = (1/n)^2
    Ωa,Ωb = ab_consts(model.puremodel[1])
    _a = [model.puremodel[i].params.a.values[1,1] for i in @comps]
    _b = [model.puremodel[i].params.b.values[1,1] for i in @comps]
    a = dot(z, _a)*invn2/Ωa
    b = dot(z, _b)*invn2/Ωb
    return a/ (b^2) # Pc mean
end


text:
using ArgParse
using Random
using JLD
using NPZ
using UUIDs

include("./MAXENT.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--temperature", "-T"
            help = "Temperature of system."
            arg_type = Float64
            default = 0.0
        "--number_of_iterations", "-N"
            help = "Number of generations before algorithm quits."
            arg_type = Int64
            default = 1000
        "--regularization_constant"
            help = "Regularization constant."
            arg_type = Float64
            default = 0.0
        "--save_file_dir"
            help = "Directory to save results in."
            arg_type = String
            default = "./maxentresults"
        "--number_of_blas_threads"
            help = "Number of BLAS threads to utilize if > 0 otherwise uses Julia default value."
            arg_type = Int64
            default = 0
        "--allow_f_increases"
            help = "Allow fitness function increases to avoid local minima."
            action = :store_true
        "qmc_data"
            help = "*.npz or *.jld file containing quantum Monte Carlo data with columns: IMAGINARY_TIME, INTERMEDIATE_SCATTERING_FUNCTION, ERROR"
            arg_type = String
            required = true
        "initial_dsf"
            help = "*.npz or *.jld file containing initial guess for dynamic structure factor and frequency with columns: FREQUENCY, DSF"
            arg_type = String
            required = true
        "default_dsf"
            help = "*.npz or *.jld file containing default dynamic structure factor with columns: FREQUENCY,DSF"
            arg_type = String
            required = true
    end

    return parse_args(s)
end

function main()
    start = time();
    parsed_args = parse_commandline()

    # create data directory
    try
        mkpath(parsed_args["save_file_dir"]);
    catch
        nothing
    end
    save_dir = parsed_args["save_file_dir"];

    _ext = splitext(parsed_args["qmc_data"])[2];
    if _ext == ".npz"
        qmcdata=NPZ.npzread(parsed_args["qmc_data"]);
    elseif _ext == ".jld"
        qmcdata = load(parsed_args["qmc_data"]);
    else
        throw(AssertionError("qmc_data must be *.jld or *.npz"));
    end

    imaginary_time = qmcdata["tau"];
    isf = qmcdata["isf"];
    isf_error = qmcdata["error"];

    _ext = splitext(parsed_args["initial_dsf"])[2];
    if _ext == ".npz"
        dsfdata=NPZ.npzread(parsed_args["initial_dsf"]);
    elseif _ext == ".jld"
        dsfdata = load(parsed_args["initial_dsf"]);
    else
        throw(AssertionError("initial_dsf must be *.jld or *.npz"));
    end
    dsf = dsfdata["dsf"];
    frequency_bins = dsfdata["frequency"];

    _ext = splitext(parsed_args["default_dsf"])[2];
    if _ext == ".npz"
        defaultdsfdata=NPZ.npzread(parsed_args["default_dsf"]);
    elseif _ext == ".jld"
        defaultdsfdata = load(parsed_args["default_dsf"]);
    else
        throw(AssertionError("default_dsf must be *.jld or *.npz"));
    end
    default_dsf = defaultdsfdata["dsf"];

    u4 = uuid4();
    results, quality_of_fit = MAXENT.maxent(dsf,default_dsf,isf,isf_error,
                           frequency_bins,imaginary_time,
                           temperature = parsed_args["temperature"],
                           regularization_constant = parsed_args["regularization_constant"],
                           number_of_iterations = parsed_args["number_of_iterations"],
                           allow_f_increases = parsed_args["allow_f_increases"],
                           number_of_blas_threads = parsed_args["number_of_blas_threads"])
    elapsed = time() - start;
    regularization_constant=parsed_args["regularization_constant"];
    filename = "$(save_dir)/maxent_results_$(regularization_constant)_$u4.jld";
    println("Saving results to $filename");
    save(filename,
         "u4",u4,
         "dsf",results.minimizer,
         "frequency",frequency_bins,
         "elapsed_time",elapsed,
         "quality_of_fit",quality_of_fit,
         "regularization_constant",regularization_constant,
         "iterations",results.iterations);
    filename = "$(save_dir)/maxent_params_$(regularization_constant)_$u4.jld";
    println("Saving parameters to $filename");
    save(filename,
         "u4",u4,
         "imaginary_time",imaginary_time,
         "isf",isf,
         "isf_error",isf_error,
         "frequency",frequency_bins,
         "initial_dsf",dsf,
         "default_dsf",default_dsf,
         "number_of_iterations",parsed_args["number_of_iterations"],
         "temperature",parsed_args["temperature"],
         "regularization_constant",parsed_args["regularization_constant"])
    nothing
end

main()


text:
abstract RealInterval     <: Real

abstract AbstractBoundary
abstract BoundingState    <: AbstractBoundary
abstract ExactState          <: BoundingState
abstract ClosedState         <: BoundingState
abstract OpenState           <: BoundingState

# singleton types
type Exact  <: ExactState  end; ExactBound  = Exact()
type Closed <: ClosedState end; ClosedBound = Closed()
type Open   <: OpenState   end; OpenBound   = Open()

abstract RealInterval     <: Real

abstract AbstractBoundary
abstract BoundingState    <: AbstractBoundary
abstract ExactState          <: BoundingState
abstract ClosedState         <: BoundingState
abstract OpenState           <: BoundingState

# singleton types
type Exact  <: ExactState  end; ExactBound  = Exact()
type Closed <: ClosedState end; ClosedBound = Closed()
type Open   <: OpenState   end; OpenBound   = Open()

type FlexibleInterval{LoState,HiState} <: Real
    lo::Real
    hi::Real
end

# fi = FlexibleIterval{LoState,HiState}(1.0,2.0)
# FlexibleInterval{Exact,Closed}(1.0, 2.0)


text:
# Toy problem to test the solver

using luteos

P = P1() # Polynomial order of solution

mesh   = Mesh2D( "square", P, N = 11)
master = Master2D( P )

mat = Material( E = 1, ν = 0.0, dim = mesh.dim )

# source function
function funcS( p::Array{Float64} )
  return fill(1.0, size(p,1), 2)
end

function funcB( p::Array{Float64} )
  return fill(0.0, size(p,1), 2)
end

prob = Elas( "Example 1", mat, funcS, [1,1,1,1], true, [funcB, funcB, funcB, funcB] )

(uhath, uh, σh ) = hdgSolve( master, mesh, prob )

# write solution
# writeTecplot( "bla.dat", prob, mesh, uh, σh )


text:
"""
    synchronous(A::AbstractMatrix)

Calculate a synchronous two-dimensional correlation spectrum.

`A` is a dataset with multiple spectra in rows.
They don't need to be at equally spaced measurements of the perturbation.
Normally, the data is given as the difference from a reference spectrum (see
below).

# Examples
```jldoctest
julia> A = rand(50, 1000);

julia> Φ = synchronous(A .- mean(A, dims=1));

julia> size(Φ)
(1000, 1000)
```
"""
function synchronous(A::AbstractMatrix)
    return A'A / (size(A, 1) - 1)
end