text:
##
# layer 3 access
# read data as records which are named tuple representations of the schema

##
# Column cursor iterates through all values of the column, including null values.
# Each iteration returns the value (as a Union{T,Nothing}), definition level, and repetition level for each value.
# Row can be deduced from repetition level.
mutable struct ColCursor{T}
    par::Parquet.File
    colname::Vector{String}                                 # column name (full path in schema)
    colnamesym::Vector{Symbol}                              # column name converted to symbols
    required_at::Vector{Bool}                               # at what level in the schema the column is required
    repeated_at::Vector{Bool}                               # at what level in the schema the column is repeated
    logical_converter_fn::Function                          # column converter function as per schema or options (identity if none)
    maxdefn::Int32                                          # maximum definition level of the column as per schema

    row::Int                                                # current row
    rows::UnitRange{Int}                                    # row range to limit to
    rowgroups::Vector{RowGroup}                             # list of row groups
    row_positions::Vector{Int64}                            # starting positions if each rowgroup
    rgidx::Int                                              # current rowgroup
    colchunks::Union{Vector{ColumnChunk},Nothing}           # list of column chunk in the current rowgroup (metadata only)
    ccidx::Int                                              # index of column chunk that contains the current position

    pageiter::Union{Nothing,ColumnChunkPageValues{T}}       # iterator for the current page
    pagerange::Union{UnitRange{Int64},Nothing}              # range of positions (rows) that the current page contains
    pageiternext::Int64
    pagevals::OutputState{T}
    page_defn::OutputState{Int32}
    page_repn::OutputState{Int32}

    valpos::Int64                                           # current position within values of the current page
    levelpos::Int64                                         # current position within levels of the current page
    levelend::Int64

    function ColCursor{T}(par::Parquet.File, row_positions::Vector{Int64}, colname::Vector{String}, rows::UnitRange) where T
        colnamesym = [Symbol(name) for name in colname]
        sch = schema(par)

        len_colname_parts = length(colname)
        required_at = Array{Bool}(undef, len_colname_parts)
        repeated_at = Array{Bool}(undef, len_colname_parts)
        for idx in 1:len_colname_parts
            partname = colname[1:idx]
            required_at[idx] = isrequired(sch, partname)         # determine whether field is optional and repeated
            repeated_at[idx] = isrepeated(sch, partname)
        end

        logical_converter_fn = logical_converter(sch, colname)
        maxdefn = max_definition_level(sch, colname)
        new{T}(par, colname, colnamesym, required_at, repeated_at, logical_converter_fn, maxdefn, 0, rows, rowgroups(par), row_positions, 0, nothing, 0, nothing, nothing, 0)
    end
end

function ColCursor(par::Parquet.File, colname::Vector{String}; rows::UnitRange=1:nrows(par), row::Signed=first(rows))
    row_positions = rowgroup_row_positions(par)
    @assert last(rows) <= nrows(par)
    @assert first(rows) >= 1
    T = elemtype(schema(par), colname)
    cursor = ColCursor{T}(par, row_positions, colname, rows)
    setrow(cursor, Int64(row))
    cursor
end

eltype(::Type{ColCursor{T}}) where {T} = NamedTuple{(:value, :defn_level, :repn_level),Tuple{Union{Nothing,T},Int64,Int64}}
length(cursor::ColCursor) = length(cursor.rows)

function setrow(cursor::ColCursor{T}, row::Int64) where {T}
    # check if cursor is done
    if ((cursor.row > 0) && !(cursor.row in cursor.rows)) || isempty(cursor.rowgroups)
        cursor.colchunks = nothing
        cursor.ccidx = 0
        cursor.valpos = cursor.levelpos = 0
        cursor.levelend = -1
        return
    end

    # check if we are already on that row
    (cursor.row == row) && return

    par = cursor.par
    prevrow = cursor.row
    cursor.row = row

    # advance row group if needed
    if (cursor.rgidx == 0) || (cursor.row_positions[cursor.rgidx+1] <= row)
        rgidx = findfirst(x->x>row, cursor.row_positions)
        if rgidx === nothing
            cursor.rgidx = length(cursor.rowgroups) + 1
            return
        else
            cursor.rgidx = rgidx - 1
            rg = cursor.rowgroups[cursor.rgidx]
            cursor.colchunks = columns(par, rg, cursor.colname)
            cursor.pagerange = nothing
            cursor.pageiternext = 0
            cursor.pageiter = nothing
            cursor.ccidx = 0
        end
    end

    # advance column chunk and page within column  chunk if needed
    while (cursor.pagerange === nothing) || !(row in cursor.pagerange)
        cursor.valpos = cursor.levelpos = 0
        cursor.levelend = -1
        startrow = (cursor.pagerange === nothing) ? cursor.row_positions[cursor.rgidx] : (last(cursor.pagerange) + 1)
        if cursor.pageiter === nothing
            # need to start a page cursor for a new column chunk
            cursor.ccidx += 1
            cursor.pageiter = ColumnChunkPageValues(par, cursor.colchunks[cursor.ccidx], T, cursor.logical_converter_fn)
            cursor.pageiternext = 0
        end
        page_iter_result = (cursor.pageiternext > 0) ? iterate(cursor.pageiter, cursor.pageiternext) : iterate(cursor.pageiter)
        if page_iter_result === nothing
            # need to start on a new column chunk
            cursor.pageiter = nothing
        else
            resultdata,cursor.pageiternext = page_iter_result
            cursor.pagevals = resultdata.value
            cursor.page_repn = resultdata.repn_level
            cursor.page_defn = resultdata.defn_level
            nrowspage = 0
            if cursor.pageiter.has_repn_levels
                for i in 1:cursor.page_repn.offset
                    (cursor.page_repn.data[i] !== Int32(0)) && (nrowspage += 1)
                end
            else
                if cursor.pageiter.has_defn_levels
                    nrowspage = cursor.page_defn.offset
                else
                    nrowspage = cursor.pagevals.offset # number of values is number of rows
                end
            end
            if nrowspage > 0
                cursor.pagerange = startrow:(startrow+nrowspage-1)
            end
        end
    end

    # advance value and level positions within the page
    if cursor.pageiter.has_defn_levels
        if cursor.levelpos == 0
            for idx in 1:(row - first(cursor.pagerange) + 1)
                cursor.levelpos += 1
                (cursor.page_defn.data[cursor.levelpos] === cursor.maxdefn) && (cursor.valpos += 1)
            end
        elseif cursor.page_defn.data[cursor.levelpos] === cursor.maxdefn
            cursor.valpos += 1
        end
    else
        # all entries are required, so there must be a corresponding value
        cursor.valpos = cursor.levelpos = row - first(cursor.pagerange) + 1
    end

    if cursor.pageiter.has_repn_levels
        # multiple entries may constitute one row
        if cursor.levelend == -1
            cursor.levelend = Int64(something(findnext(x->x===Int32(0), cursor.page_repn.data, cursor.levelpos+1), cursor.page_repn.offset+1) - 1)
        end
    else
        # no repetitions, so each entry corresponds to one full row
        cursor.levelend = cursor.levelpos
    end
    nothing
end

function _start(cursor::ColCursor)
    setrow(cursor, Int64(first(cursor.rows)))
    cursor.row, cursor.levelpos
end
function _done(cursor::ColCursor, rowandlevel::Tuple{Int64,Int64})
    row, levelpos = rowandlevel
    (levelpos > cursor.levelend) || !(row in cursor.rows)
end
function _next(cursor::ColCursor{T}, rowandlevel::Tuple{Int64,Int64}) where {T}
    # find values for current row and level in row
    row, levelpos = rowandlevel
    (levelpos == cursor.levelpos) || throw(InvalidStateException("Invalid column cursor state", :levelpos))

    maxdefn = cursor.maxdefn
    defn_level = cursor.pageiter.has_defn_levels ? cursor.page_defn.data[cursor.levelpos] : maxdefn
    repn_level = cursor.pageiter.has_repn_levels ? cursor.page_repn.data[cursor.levelpos] : Int32(0)
    if defn_level == maxdefn
        val = (cursor.pagevals.data[cursor.valpos])::T
    else
        val = nothing
    end

    # advance row
    cursor.levelpos += 1
    if cursor.levelpos > cursor.levelend
        row += 1
        cursor.levelend = -1
        setrow(cursor, Int64(row))
    end

    NamedTuple{(:value, :defn_level, :repn_level),Tuple{Union{Nothing,T},Int64,Int64}}((val, defn_level, repn_level)), (row, cursor.levelpos)
end

function Base.iterate(cursor::ColCursor{T}, state) where {T}
    _done(cursor, state) && return nothing
    return _next(cursor, state)
end

function Base.iterate(cursor::ColCursor)
    r = iterate(cursor, _start(cursor))
    return r
end

##

mutable struct BatchedColumnsCursor{T}
    par::Parquet.File
    colnames::Vector{Vector{String}}
    colcursors::Vector{ColCursor}
    colstates::Vector{Tuple{Int64,Int64}}
    colbuffs::Vector{Union{Nothing,Vector}}
    batchid::Int
    rows::UnitRange{Int64}
    row::Int64
    batchsize::Int64
    nbatches::Int
    reusebuffer::Bool
    use_threads::Bool
end

"""
Create cursor to iterate over batches of column values. Each iteration returns a named tuple of column names with batch of column values. Files with nested schemas can not be read with this cursor.

```julia
BatchedColumnsCursor(par::Parquet.File; kwargs...)
```

Cursor options:
- `rows`: the row range to iterate through, all rows by default.
- `batchsize`: maximum number of rows to read in each batch (default: row count of first row group).
- `reusebuffer`: boolean to indicate whether to reuse the buffers with every iteration; if each iteration processes the batch and does not need to refer to the same data buffer again, then setting this to `true` reduces GC pressure and can help significantly while processing large files.
- `use_threads`: whether to use threads while reading the file; applicable only for Julia v1.3 and later and switched on by default if julia processes is started with multiple threads.
"""
function BatchedColumnsCursor(par::Parquet.File;
        rows::UnitRange=1:nrows(par),
        batchsize::Signed=length(rows) > 0 ? min(length(rows), first(rowgroups(par)).num_rows) : 0,
        reusebuffer::Bool=false,
        use_threads::Bool=(nthreads() > 1))

    sch = schema(par)

    # supports only non nested columns as of now
    if !all(num_children(schemaelem) == 0 for schemaelem in sch.schema[2:end])
        error("nested schemas are not supported with BatchedColumnsCursor yet")
    end

    colcursors = [ColCursor(par, colname; rows=rows) for colname in colnames(par)]
    rectype = ntcolstype(sch, sch.schema[1])
    nbatches = batchsize > 0 ? ceil(Int, length(rows)/batchsize) : 0
    colbuffs = Union{Nothing,Vector}[nothing for idx in 1:length(colcursors)]

    BatchedColumnsCursor{rectype}(par, colnames(par), colcursors, Array{Tuple{Int64,Int64}}(undef, length(colcursors)), colbuffs, 1, rows, first(rows), batchsize, nbatches, reusebuffer, (VERSION < v"1.3") ? false : use_threads)
end

eltype(::Type{BatchedColumnsCursor{T}}) where {T} = T
length(cursor::BatchedColumnsCursor) = cursor.nbatches

function colcursor_advance(colcursor::ColCursor, rows_by::Int64, vals_by::Int64=rows_by)
    if (colcursor.row + rows_by) > last(colcursor.pagerange)
        setrow(colcursor, Int64(colcursor.row+rows_by))
    else
        colcursor.valpos += vals_by
        colcursor.row += rows_by
        colcursor.levelpos += rows_by
    end
    nothing
end

function colcursor_values(colcursor::ColCursor{T}, batchsize::Int64, ::Type{Vector{Union{Missing,T}}}, cache) where {T}
    row = colcursor.row
    batchsize = min(batchsize, last(colcursor.rows)-row+1)
    vals = (cache === nothing) ? Array{Union{Missing,T}}(undef, batchsize) : resize!(cache::Vector{Union{Missing,T}}, batchsize)

    fillpos = 1
    while fillpos <= batchsize
        pagevals = colcursor.pagevals
        defn_levels = colcursor.page_defn
        val_idx = (colcursor.valpos == 0) ? 0 : (colcursor.valpos-1)
        nvals_from_page = min(batchsize - fillpos + 1, defn_levels.offset - colcursor.levelpos + 1)
        @inbounds for idx in 1:nvals_from_page
            if defn_levels.data[colcursor.levelpos+idx-1] === Int32(1)
                vals[fillpos+idx-1] = pagevals.data[val_idx+=1]
            else
                vals[fillpos+idx-1] = missing
            end
        end
        fillpos += nvals_from_page
        valposincr = (colcursor.valpos == 0) ? val_idx : (val_idx - colcursor.valpos + 1)
        colcursor_advance(colcursor, nvals_from_page, Int64(valposincr))
    end
    vals
end

function colcursor_values(colcursor::ColCursor{T}, batchsize::Int64, ::Type{Vector{T}}, cache) where {T}
    row = colcursor.row
    batchsize = min(batchsize, last(colcursor.rows)-row+1)
    vals = (cache === nothing) ? Array{T}(undef, batchsize) : resize!(cache::Vector{T}, batchsize)

    fillpos = 1
    while fillpos <= batchsize
        pagevals = colcursor.pagevals
        nvals_from_page = min(batchsize - fillpos + 1, pagevals.offset - colcursor.valpos + 1)
        @inbounds for idx in 1:nvals_from_page
            vals[fillpos+idx-1] = pagevals.data[pagevals.offset+idx]
        end
        fillpos += nvals_from_page
        colcursor_advance(colcursor, nvals_from_page)
    end
    vals
end

function Base.iterate(cursor::BatchedColumnsCursor{T}, batchid) where {T}
    (batchid > length(cursor)) && (return nothing)

    colcursors = cursor.colcursors
    coltypes = T.types
    if cursor.use_threads
        L = length(colcursors)
        colvals = Array{Any}(undef, L)
        @threads for idx in 1:L
            colbuff = cursor.reusebuffer ? cursor.colbuffs[idx] : nothing
            colvals[idx] = colcursor_values(colcursors[idx], cursor.batchsize, coltypes[idx], colbuff)
        end
        cursor.reusebuffer && (cursor.colbuffs = colvals)
    else
        if cursor.reusebuffer
            colvals = cursor.colbuffs = [colcursor_values(colcursor,cursor.batchsize,coltype,colbuff) for (colcursor,coltype,colbuff) in zip(colcursors,coltypes,cursor.colbuffs)]
        else
            colvals = [colcursor_values(colcursor,cursor.batchsize,coltype,nothing) for (colcursor,coltype) in zip(colcursors,coltypes)]
        end
    end

    cursor.row += cursor.batchsize
    cursor.batchid += 1
    T(colvals), cursor.batchid
end

function Base.iterate(cursor::BatchedColumnsCursor{T}) where {T}
    cursor.row = first(cursor.rows)
    for colcursor in cursor.colcursors
        setrow(colcursor, Int64(cursor.row))
    end
    iterate(cursor, cursor.batchid)
end

##

mutable struct RecordCursor{T}
    par::Parquet.File
    colnames::Vector{Vector{String}}
    colcursors::Vector{ColCursor}
    colstates::Vector{Tuple{Int64,Int64}}
    rows::UnitRange{Int64}                      # rows to scan over
    row::Int64                                  # current row
end

"""
Create cursor to iterate over records. In parallel mode, multiple remote cursors can be created and iterated on in parallel.

```julia
RecordCursor(par::Parquet.File; kwargs...)
```

Cursor options:
- `rows`: the row range to iterate through, all rows by default.
- `colnames`: the column names to retrieve; all by default
"""
function RecordCursor(par::Parquet.File; rows::UnitRange=1:nrows(par), colnames::Vector{Vector{String}}=colnames(par))
    colcursors = [ColCursor(par, colname; rows=rows, row=first(rows)) for colname in colnames]
    sch = schema(par)
    rectype = ntelemtype(sch, sch.schema[1])
    RecordCursor{rectype}(par, colnames, colcursors, Array{Tuple{Int64,Int64}}(undef, length(colcursors)), rows, first(rows))
end

eltype(::Type{RecordCursor{T}}) where {T} = T
length(cursor::RecordCursor) = length(cursor.rows)

function Base.iterate(cursor::RecordCursor{T}, row) where {T}
    (row > last(cursor.rows)) && (return nothing)

    states = cursor.colstates
    cursors = cursor.colcursors

    colvals = Dict{Symbol,Any}()
    col_repeat_state = Dict{Tuple{Int,Int},Int}()
    for colid in 1:length(states)  # for each column
        colcursor = cursors[colid]
        colstate = states[colid]
        states[colid] = update_record(cursor.par, colvals, colid, colcursor, colstate, col_repeat_state)
    end
    cursor.row += 1
    _nt(colvals, T), cursor.row
end

function Base.iterate(cursor::RecordCursor{T}) where {T}
    cursor.row = first(cursor.rows)
    cursor.colstates = [_start(colcursor) for colcursor in cursor.colcursors]
    iterate(cursor, cursor.row)
end

function _val_or_missing(dict::Dict{Symbol,Any}, k::Symbol, ::Type{T}) where {T}
    v = get(dict, k, missing)
    if isa(v, Vector)
        elt = eltype(v)
        if Dict{Symbol,Any} <: elt
            nonmissing_elt = nonmissingtype(eltype(T))
            v = [el === missing ? el : _nt(el,nonmissing_elt) for el in v]
        end
    end
    (isa(v, Dict{Symbol,Any}) ? _nt(v, nonmissingtype(T)) : v)::T
end

@generated function _nt(dict::Dict{Symbol,Any}, ::Type{T}) where {T}
    names = fieldnames(T)
    strnames = ["$n" for n in names]
    quote
        return T(($([:(_val_or_missing(dict,Symbol($(strnames[i])),$(fieldtype(T,i)))) for i in 1:length(names)]...),))
    end
end

default_init(::Type{Vector{T}}) where {T} = Vector{T}()
default_init(::Type{Dict{Symbol,Any}}) = Dict{Symbol,Any}()
default_init(::Type{T}) where {T} = ccall(:jl_new_struct_uninit, Any, (Any,), T)::T

function update_record(par::Parquet.File, row::Dict{Symbol,Any}, colid::Int, colcursor::ColCursor, colcursor_state::Tuple{Int64,Int64}, col_repeat_state::Dict{Tuple{Int,Int},Int})
    colpos = colcursor.row
    # iterate all repeated values from the column cursor (until it advances to the next row)
    while !_done(colcursor, colcursor_state)
        colval, colcursor_state = _next(colcursor, colcursor_state)                                                        # for each value, defn level, repn level in column
        update_record(par, row, colid, colcursor, colval.value, colval.defn_level, colval.repn_level, col_repeat_state)    # update record
        (colcursor.row > colpos) && break
    end
    colcursor_state # return new colcursor state
end

function update_record(par::Parquet.File, row::Dict{Symbol,Any}, colid::Int, colcursor::ColCursor, val, defn_level::Int64, repn_level::Int64, col_repeat_state::Dict{Tuple{Int,Int},Int})
    nameparts = colcursor.colname
    symnameparts = colcursor.colnamesym
    required_at = colcursor.required_at
    repeated_at = colcursor.repeated_at

    lparts = length(nameparts)
    sch = par.schema
    F = row  # the current field corresponding to the level in nameparts
    Fdefn = 0
    Frepn = 0

    # for each name part of colname (a field)
    for idx in 1:lparts
        colname = view(nameparts, 1:idx)
        #@debug("updating part $colname of $nameparts isnull:$(val === nothing), def:$(defn_level), rep:$(repn_level)")
        leaf = nameparts[idx]
        symleaf = symnameparts[idx]

        required = required_at[idx]                 # determine whether field is optional and repeated
        repeated = repeated_at[idx]
        required || (Fdefn += 1)                    # if field is optional, increment defn level
        repeated && (Frepn += 1)                    # if field can repeat, increment repn level

        defined = ((val === nothing) || (idx < lparts)) ? haskey(F, symleaf) : false
        mustdefine = defn_level >= Fdefn
        mustrepeat = repeated && (repn_level <= Frepn)
        repidx = 0
        if mustrepeat
            repkey = (colid, idx)
            repidx = get(col_repeat_state, repkey, 0)
            repidx += 1
            col_repeat_state[repkey] = repidx
        end
        nreps = (defined && isa(F[symleaf], Vector)) ? length(F[symleaf]) : 0

        #@debug("repeat:$mustrepeat, nreps:$nreps, repidx:$repidx, defined:$defined, mustdefine:$mustdefine")
        if mustrepeat && (nreps < repidx)
            if !defined && mustdefine
                Vtyp = elemtype(sch, colname)
                Vrep = F[symleaf] = default_init(Vtyp)
            else
                Vrep = F[symleaf]
            end
            if length(Vrep) < repidx
                resize!(Vrep, repidx)
                if !isbits(eltype(Vrep))
                    Vrep[repidx] = default_init(eltype(Vrep))
                end
            end
            F = Vrep[repidx]
        elseif !defined && mustdefine
            if idx == length(nameparts)
                V = val
            else
                Vtyp = elemtype(sch, colname)
                V = default_init(Vtyp)
            end
            F[symleaf] = V
            F = V
        else
            F = get(F, symleaf, nothing)
        end
    end
    nothing
end


text:
immutable Gumbel <: ContinuousUnivariateDistribution
    μ::Float64   # location
    β::Float64 # scale

    function Gumbel(μ::Real, β::Real)
        β > zero(β) || error("The scale of Gumbel must be positive")
        @compat new(Float64(μ), Float64(β))
    end

    Gumbel(μ::Real) = Gumbel(μ, 1.0)
    Gumbel() = new(0.0, 1.0)
end

@distr_support Gumbel -Inf Inf

const DoubleExponential = Gumbel

#### Parameters

location(d::Gumbel) = d.μ
scale(d::Gumbel) = d.β
params(d::Gumbel) = (d.μ, d.β)

#### Statistics

mean(d::Gumbel) = d.μ + d.β * 0.57721566490153286

median(d::Gumbel) = d.μ + d.β * 0.366512920581664327

mode(d::Gumbel) = d.μ

var(d::Gumbel) = 1.6449340668482264 * d.β^2

skewness(d::Gumbel) = 1.13954709940464866

kurtosis(d::Gumbel) = 2.4

entropy(d::Gumbel) = 1.57721566490153286 + log(d.β)

#### Evaluation

zval(d::Gumbel, x::Float64) = (x - d.μ) / d.β
xval(d::Gumbel, z::Float64) = x * d.β + d.μ

function pdf(d::Gumbel, x::Float64)
    z = zval(d, x)
    exp(-z - exp(-z)) / d.β
end

function logpdf(d::Gumbel, x::Float64)
    z = zval(d, x)
    - (z + exp(-z) + log(d.β))
end

cdf(d::Gumbel, x::Float64) = exp(-exp(-zval(d, x)))
logcdf(d::Gumbel, x::Float64) = -exp(-zval(d, x))

quantile(d::Gumbel, p::Float64) = d.μ - d.β * log(-log(p))

gradlogpdf(d::Gumbel, x::Float64) = - (1.0 + exp((d.μ - x) / d.β)) / d.β

#### Sampling

rand(d::Gumbel) = quantile(d, rand())

ext: jl


text:
"""
    generate_mixture_normal(n, m, p)

generate `n` observations from mixture of univariate normals
each with standard deviation 1 and mean parameters `m` and proportion `p`.
"""
function generate_mixture_normal(n::Int64, m::Array{T,1}=[-20., 0., 20.], p::Array{T,1}=ones(length(m))/length(m)) where T <: AbstractFloat
    mixed_normal = Array{Float64}(undef, 0)  # placeholder for the resulting mixture
    n_mixtures = Int64.(floor.(n * p))  # sample sizes for each univariate normals

    # leftover n's is assigned to the first distribution > sum(n_mixtures)
    n_mixtures[1] += (n - sum(n_mixtures))

    # sample from Normal(0, 1)
    for i=1:length(m)
        mixed_normal = [mixed_normal; randn(n_mixtures[i]) .+ m[i]]
    end

    mixed_normal
end


text:
function koopman(g,prob,u0,p,args...;kwargs...)
    g(solve(remake(prob,u0=u0,p=p),args...;kwargs...))
end

function koopman_expectation(g,u0s,ps,prob,args...;maxiters=0,
                      batch = 0,
                      quadalg = HCubatureJL(),
                      ireltol = 1e-2, iabstol=1e-2,kwargs...)
  n = length(u0s)
  if batch == 0
    _f = function (x,p)
      u0 = x[1:n]
      p  = x[n+1:end]
      k = koopman(g,prob,u0,p,args...;kwargs...)
      w = prod(pdf(a,b) for (a,b) in zip(u0s,u0))*
          prod(pdf(a,b) for (a,b) in zip(ps,p))
      k*w
    end
  else
    _f = function (dx,x,p)
      trajectories = size(x,2)
      prob_func = (prob,i,repeat) -> remake(prob,u0=@view(x[1:n,i]),
                                                 p=@view(x[n+1:end,i]))
      output_func = function (sol,i)
        k = g(sol)
        u0= @view(x[1:n,i])
        p = @view(x[n+1:end,i])
        w = prod(pdf(a,b) for (a,b) in zip(u0s,u0))*
            prod(pdf(a,b) for (a,b) in zip(ps,p))
        k*w,false
      end

      ensembleprob = EnsembleProblem(prob,prob_func=prob_func,
                                     output_func = output_func)
      sol = solve(ensembleprob,args...;trajectories=trajectories,kwargs...)
      dx .= vec(sol.u)
      nothing
    end
  end
  xs = [u0s;ps]
  intprob = QuadratureProblem(_f,minimum.(xs),maximum.(xs),batch=batch)
  sol = solve(intprob,quadalg,reltol=ireltol,
              abstol=iabstol,maxiters = maxiters)
end

function montecarlo_expectation(g,u0s,ps,prob,args...;trajectories,kwargs...)
  prob_func = function (prob,i,repeat)
    remake(prob,u0=rand.(u0s),p=rand.(ps))
  end
  output_func = (sol,i) -> (g(sol),false)
  monte_prob = EnsembleProblem(prob;
                                 output_func = output_func,
                                 prob_func = prob_func)
  mean(solve(monte_prob,args...;trajectories=trajectories,kwargs...).u)
end


text:
using Distributed
addprocs(2)

function main()
    if length(ARGS) != 2
        println("Usage: ping_pong <# iterations> <vector_length>")
        exit(1)
    end

    # Parse parameters
    argv = map(x -> parse(Int, x), ARGS)
    iterations, payload = argv

    # Setup the synchronous channels
    channels = [
        RemoteChannel(()->Channel{Vector{Float64}}(0), proc)
        for proc in workers()
    ]

    @sync for proc in workers()
        @spawnat proc work(iterations, payload, channels)
    end
end

@everywhere function work(iterations, payload, channels)
    partner_rank = if (myid() == 2) 3 else 2 end
    local_channel_ref = channels[myid()-1]
    other_channel_ref = channels[partner_rank-1]

    vector = Vector{Float64}(undef, payload)

    local t0, t1
    for k in 1 : iterations * 2
        if k == iterations + 1
            t0 = time_ns()
        end
        if k % 2 == myid() - 2
            j = k % payload + 1
            i = if (j != 0) j else payload end
            vector[i] += Float64(k)
            put!(other_channel_ref, vector)
        else
            vector = take!(local_channel_ref)
        end
    end
    t1 = time_ns()

    throughput = iterations * payload * 8.0 * 1e-6 / ((t1-t0)*1e-9)

    if myid() == 2
        println("Rate (MB/s): $throughput")
    end
end

main()


text:
# A quadratic irrational : (√a + b)/c
type QuadraticIrrational
  a_square::Integer
  b::Integer
  c::Integer
end

qi_sqrt(a::Integer) = QuadraticIrrational(a, 0, 1)

eval(q::QuadraticIrrational) = (sqrt(q.a_square) + q.b) / q.c

function -(q::QuadraticIrrational, i::Integer)
  return QuadraticIrrational(q.a_square, q.b - q.c * i, q.c)
end

function +(q::QuadraticIrrational, i::Integer)
  return QuadraticIrrational(q.a_square, q.b + q.c * i, q.c)
end

function /(i::Integer, q::QuadraticIrrational)
  # (i*c)/(√a + b) = ((i*c)*(√a + b))/(a - b^2)
  den = q.a_square - q.b^2
  num_f = i * q.c
  common_div = gcd(num_f, den)
  den /= common_div
  num_f /= common_div
  return QuadraticIrrational(num_f^2 * q.a_square, -num_f * q.b, den)
end

function show(io::IO, q::QuadraticIrrational)
  print(io, "(√$(q.a_square) + $(q.b)) / $(q.c)")
end


text:
module BaytesCore

############################################################################################
#Import External packages
import Base: Base, split, resize!
import Random: Random, shuffle!
import LogExpFunctions: LogExpFunctions, logsumexp, logaddexp

using DocStringExtensions:
    DocStringExtensions, TYPEDEF, TYPEDFIELDS, FIELDS, SIGNATURES, FUNCTIONNAME
using ArgCheck: ArgCheck, @argcheck
using UnPack: @unpack, @pack!
using Random: Random, AbstractRNG, GLOBAL_RNG
using Statistics: Statistics, middle

############################################################################################
#Import
include("Algorithm/Algorithm.jl")
include("Core/Core.jl")
include("Tempering/Tempering.jl")
include("Buffer/Buffer.jl")
include("Sampling/Sampling.jl")

############################################################################################
#export

end


text:
# using GLGM06 + template hull intersection
Dict([
    "opts" => Dict(:T=>3.0, :alg=>GLGM06(δ=1e-2, max_order=10),
                   :intersection_method=>TemplateHullIntersection(BoxDirections(2)),
                   :clustering_method=>ZonotopeClustering(),
                   :intersect_source_invariant=>false),
    "X0"   => Hyperrectangle(low=[10.0, 0.0], high=[10.2, 0.0]),
    "info" => Dict(
        "name" => "bouncing_ball",
        "dim" => 2,
        "linear" => true,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# Autogenerated wrapper script for Ncurses_jll for aarch64-apple-darwin
export libform, libmenu, libncurses, libpanel

JLLWrappers.@generate_wrapper_header("Ncurses")
JLLWrappers.@declare_library_product(libform, "@rpath/libformw.6.dylib")
JLLWrappers.@declare_library_product(libmenu, "@rpath/libmenuw.6.dylib")
JLLWrappers.@declare_library_product(libncurses, "@rpath/libncursesw.6.dylib")
JLLWrappers.@declare_library_product(libpanel, "@rpath/libpanelw.6.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libform,
        "lib/libform.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmenu,
        "lib/libmenu.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libncurses,
        "lib/libncurses.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpanel,
        "lib/libpanel.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ITensors,
  ITensors.NDTensors,
  ITensorGPU,
  LinearAlgebra, # For tr()
  Combinatorics, # For permutations()
  CUDA,
  Test

@testset "cuITensor $T Contractions" for T in (Float64, ComplexF64)
  mi, mj, mk, ml, ma = 2, 3, 4, 5, 6, 7
  i = Index(mi, "i")
  j = Index(mj, "j")
  k = Index(mk, "k")
  l = Index(ml, "l")
  a = Index(ma, "a")
  @testset "Test contract cuITensors" begin
    Aij = cuITensor(randomITensor(T, i, j))
    Aji = cuITensor(randomITensor(T, j, i))
    Bij = cuITensor(randomITensor(T, i, j))
    Aik = cuITensor(randomITensor(T, i, k))
    Ajk = cuITensor(randomITensor(T, j, k))
    Ajl = cuITensor(randomITensor(T, j, l))
    Akl = cuITensor(randomITensor(T, k, l))
    Dv = rand(T, mi)
    D = itensor(ITensors.tensor(NDTensors.Diag(CuVector(Dv)), IndexSet(i, i')))
    Ev = rand(T, mi)
    E = itensor(ITensors.tensor(NDTensors.Diag(CuVector(Ev)), IndexSet(i, i'')))
    @testset "Test contract cuITensors (Matrix*Diag -> Matrix)" begin
      C = Aij * D
      @test collect(CuArray(C)) ≈ collect(CuMatrix(Aij, j, i)) * diagm(0 => Dv)
    end
    @testset "Test contract cuDiagITensors (Diag*Diag -> Diag)" begin
      C = E * D
      cC = CuArray(C)
      @test collect(cC) ≈ diagm(0 => Ev) * diagm(0 => Dv)
    end
    @testset "Test contract cuDiagITensors (UniformDiag*Diag -> Diag)" begin
      scal = itensor(ITensors.tensor(NDTensors.Diag(2.0), IndexSet(i, i'')))
      C = scal * D
      @test collect(CuArray(C)) ≈ 2.0 .* diagm(0 => Dv)
      C = D * scal
      @test collect(CuArray(C)) ≈ 2.0 .* diagm(0 => Dv)
    end
    @testset "Test contract cuITensors (Matrix*UniformDiag -> Matrix)" begin
      scal = itensor(ITensors.tensor(NDTensors.Diag(T(2.0)), IndexSet(i, i')))
      C = scal * Aij
      @test cpu(C) ≈ 2.0 * cpu(replaceind(Aij, i, i')) atol = 1e-4
      C = Aij * scal
      @test_broken cpu(C) ≈ 2.0 * cpu(replaceind(permute(Aij, j, i), i, i')) atol = 1e-4
    end
  end # End contraction testset
end

@testset "cuITensor $T1, $T2 Contractions" for T1 in (Float64, ComplexF64),
  T2 in (Float64, ComplexF64)

  mi, mj, mk, ml, ma = 2, 3, 4, 5, 6, 7
  i = Index(mi, "i")
  j = Index(mj, "j")
  k = Index(mk, "k")
  l = Index(ml, "l")
  a = Index(ma, "a")
  @testset "Test contract cuITensors" begin
    Aij = cuITensor(randomITensor(T1, i, j))
    Aji = cuITensor(randomITensor(T1, j, i))
    Bij = cuITensor(randomITensor(T1, i, j))
    Dv = rand(T2, mi)
    D = itensor(ITensors.tensor(NDTensors.Diag(CuVector(Dv)), IndexSet(i, i')))
    Ev = rand(T2, mi)
    E = itensor(ITensors.tensor(NDTensors.Diag(CuVector(Ev)), IndexSet(i, i'')))
    @testset "Test contract cuITensors (Matrix*Diag -> Matrix)" begin
      C = Aij * D
      cC = CuArray(C)
      @test collect(cC) ≈ collect(CuMatrix(Aij, j, i)) * diagm(0 => Dv)
    end
    @testset "Test contract cuDiagITensors (Diag*Diag -> Diag)" begin
      C = E * D
      cC = CuArray(C)
      @test collect(cC) ≈ diagm(0 => Ev) * diagm(0 => Dv)
    end
    @testset "Test contract cuDiagITensors (UniformDiag*Diag -> Diag)" begin
      scal = itensor(ITensors.tensor(NDTensors.Diag(T2(2.0)), IndexSet(i, i'')))
      C = scal * D
      cC = CuArray(C)
      @test collect(cC) ≈ 2.0 .* diagm(0 => Dv)
      C = D * scal
      cC = CuArray(C)
      @test collect(cC) ≈ 2.0 .* diagm(0 => Dv)
    end
    @testset "Test contract cuITensors (Matrix*UniformDiag -> Matrix)" begin
      scal = itensor(ITensors.tensor(NDTensors.Diag(T2(2.0)), IndexSet(i, i')))
      C = scal * Aij
      cC = CuArray(C)
      @test collect(cC) ≈ array(2.0 * cpu(replaceind(Aij, i, i'))) atol = 1e-4
      C = Aij * scal
      cC = CuArray(C)
      @test_broken collect(cC) ≈ array(2.0 * cpu(replaceind(permute(Aij, j, i), i, i'))) atol =
        1e-4
    end
  end # End contraction testset
end


text:
mutable struct ContinuousSystem
    jacobian::SparseMatrixCSC{Float64,Int64}
    jacobianTranspose::SparseMatrixCSC{Float64,Int64}
    observation::Array{Float64,1}
    variance::Array{Float64,1}
    data::String
end

mutable struct ContinuousGraph
    Nvariable::Int64
    Nfactor::Int64
    Nindirect::Int64
    Nlink::Int64
    virtualMean::Float64
    virtualVariance::Float64
    meanDirect::Array{Float64,1}
    weightDirect::Array{Float64,1}
    meanIndirect::Array{Float64,1}
    varianceIndirect::Array{Float64,1}
    coefficient::Array{Float64,1}
    toVariable::SparseMatrixCSC{Int64, Int64}
    toFactor::SparseMatrixCSC{Int64, Int64}
    rowptr::Vector{Vector{Int64}}
    colptr::Vector{Vector{Int64}}
    colptrMarginal::Vector{Vector{Int64}}
    alphaNew::Array{Float64,1}
    alphaOld::Array{Float64,1}
    iterateFactor::Array{Int64,1}
    iterateVariable::Array{Int64,1}
    iterateMarginal::Array{Int64,1}
    dynamic::Array{Int64,1}
end

mutable struct ContinuousTreeGraph
    Nvariable::Int64
    Nfactor::Int64
    Nlink::Int64
    root::Int64
    virtualMean::Float64
    virtualVariance::Float64
    meanDirect::Array{Float64,1}
    weightDirect::Array{Float64,1}
    rowForward::Vector{Vector{Int64}}
    rowBackward::Vector{Vector{Int64}}
    colForward::Vector{Vector{Int64}}
    colBackward::Vector{Vector{Int64}}
    incomingToFactor::Vector{Vector{Int64}}
    incomingToVariable::Vector{Vector{Int64}}
    iterateFactor::Array{Int64,1}
    iterateVariable::Array{Int64,1}
    passFactorVariable::Int64
    passVariableFactor::Int64
    forward::Bool
    backward::Bool
end

mutable struct ContinuousInference
    fromFactor::Array{Int64,1}
    toVariable::Array{Int64,1}
    meanFactorVariable::Array{Float64,1}
    varianceFactorVariable::Array{Float64,1}
    fromVariable::Array{Int64,1}
    toFactor::Array{Int64,1}
    meanVariableFactor::Array{Float64,1}
    varianceVariableFactor::Array{Float64,1}
    mean::Array{Float64,1}
    variance::Array{Float64,1}
end

struct ContinuousModel
    graph::ContinuousGraph
    inference::ContinuousInference
    system::ContinuousSystem
end

struct ContinuousTreeModel
    graph::ContinuousTreeGraph
    inference::ContinuousInference
    system::ContinuousSystem
end

########## Form a graph with continuous variables and initialize messages ##########
function continuousModel(
    args...;
    prob::Float64 = 0.6,
    alpha::Float64 = 0.4,
    mean::Float64 = 0.0,
    variance::Float64 = 1e10)

    checkKeywords(prob, alpha, variance)
    if checkFileOrArguments(args)
        system = readContinuousFile(args)
    else
        system = readContinuousArguments(args)
    end

    graph, inference = makeContinuousGraph(system, mean, variance, prob, alpha)

    return ContinuousModel(graph, inference, system)
end

########## Form a tree factor graph and initialize messages ##########
function continuousTreeModel(
    args...;
    mean::Float64 = 0.0,
    variance::Float64 = 1e10,
    root::Int64 = 1)

    if checkFileOrArguments(args)
        system = readContinuousFile(args)
    else
        system = readContinuousArguments(args)
    end

    graph, inference = makeContinuousTreeGraph(system, mean, variance, root)

    return ContinuousTreeModel(graph, inference, system)
end

########## Load from HDF5 or XLSX files ##########
function readContinuousFile(args)
    fullpath, extension, dataname = checkImportFile(args)

    #### Read from HDF5 or XLSX file
    if extension == ".h5"
        list = h5read(fullpath, "/jacobian")::Array{Float64,2}
        jacobian = sparse(list[:,1], list[:,2], list[:,3])::SparseMatrixCSC{Float64,Int64}
        jacobianTranspose = sparse(list[:,2], list[:,1], list[:,3])::SparseMatrixCSC{Float64,Int64}

        observation = h5read(fullpath, "/observation")::Array{Float64,1}
        variance = h5read(fullpath, "/variance")::Array{Float64,1}
    elseif extension == ".xlsx"
        xf = XLSX.openxlsx(fullpath, mode = "r")
        if "jacobian" in XLSX.sheetnames(xf)
            start = startxlsx(xf["jacobian"])
            list = xf["jacobian"][:][start:end, :]
            jacobian = sparse(list[:, 1], list[:, 2], list[:, 3])
            jacobianTranspose = sparse(list[:, 2], list[:, 1], list[:, 3])
        else
            throw(ErrorException("error opening sheet jacobian"))
        end
        if "observation" in XLSX.sheetnames(xf)
            start = startxlsx(xf["observation"])
            observation = xf["observation"][:][start:end]
        else
            throw(ErrorException("error opening sheet observation"))
        end
        if "variance" in XLSX.sheetnames(xf)
            start = startxlsx(xf["variance"])
            variance = xf["variance"][:][start:end]
        else
            throw(ErrorException("error opening sheet variance"))
        end
    else
        error("the input data is not a valid format")
    end

    return ContinuousSystem(jacobian, jacobianTranspose, observation, variance, dataname)
end

########## Read in-Julia system model ##########
function readContinuousArguments(args)
    if typeof(args[1]) == Array{Float64, 2}
        jacobian = sparse(args[1])
    else
        jacobian = args[1]
    end

    jacobianTranspose = copy(transpose(jacobian))
    observation = args[2]
    variance = args[3]

    return ContinuousSystem(jacobian, jacobianTranspose, observation, variance, "noname")
end

########## Produce the graphical model ##########
function makeContinuousGraph(system, meanVirtual, varianceVirtual, dampProbability, dampAlpha)
    ### Number of factor and variable nodes
    Nfactor, Nvariable = size(system.jacobian)

    ### Find graph numbers, set the direct mean and variance and internal factor numeration
    Ndirect = 0; Nlink = 0; Nindirect = 0
    meanDirect = fill(meanVirtual / varianceVirtual, Nvariable)
    weightDirect = fill(1 / varianceVirtual, Nvariable)
    dynamic = fill(0, Nfactor)

    @inbounds for i = 1:Nfactor
        NvariableInRow = system.jacobianTranspose.colptr[i + 1] - system.jacobianTranspose.colptr[i]
        if NvariableInRow == 1
            Ndirect += NvariableInRow
            variable = system.jacobianTranspose.rowval[system.jacobianTranspose.colptr[i]]
            meanDirect[variable] = 0.0
            weightDirect[variable] = 0.0
        else
            Nlink += NvariableInRow
            Nindirect += 1
            dynamic[i] = Nindirect
        end
    end

    ### Pass through the columns
    colptr = [Int[] for i = 1:Nvariable]
    colptrMarginal = [Int[] for i = 1:Nvariable]
    toVariable = fill(0, Nlink)
    fromFactor = similar(toVariable)
    idxi = 1
    @inbounds for col = 1:Nvariable
        for i = system.jacobian.colptr[col]:(system.jacobian.colptr[col + 1] - 1)
            row = system.jacobian.rowval[i]
            NvariableInRow = system.jacobianTranspose.colptr[row + 1] - system.jacobianTranspose.colptr[row]
            if NvariableInRow == 1
                meanDirect[col] += system.observation[row] * system.jacobian[row, col] / system.variance[row]
                weightDirect[col] += system.jacobian[row, col]^2 / system.variance[row]
            else
                push!(colptr[col], idxi)
                push!(colptrMarginal[col], idxi)
                toVariable[idxi] = col
                fromFactor[idxi] = row
                idxi += 1
            end
        end
    end

    ### Pass through the rows and send messages from singly-connected factor nodes to all indirect links
    coefficient = fill(0.0, Nlink)
    meanIndirect = similar(coefficient)
    varianceIndirect = similar(coefficient)
    meanFactorVariable = similar(coefficient)
    varianceFactorVariable = similar(coefficient)
    meanVariableFactor = similar(coefficient)
    varianceVariableFactor = similar(coefficient)

    fromVariable = fill(0, Nlink)
    toFactor = similar(fromVariable)
    rowptr = [Int[] for i = 1:Nindirect]
    idxi = 1
    @inbounds for (col, val) in enumerate(dynamic)
        if val != 0
            for i = system.jacobianTranspose.colptr[col]:(system.jacobianTranspose.colptr[col + 1] - 1)
                row = system.jacobianTranspose.rowval[i]

                coefficient[idxi] = system.jacobianTranspose[row, col]
                meanIndirect[idxi] = system.observation[col]
                varianceIndirect[idxi] = system.variance[col]

                toFactor[idxi] = col
                fromVariable[idxi] = row
                push!(rowptr[dynamic[col]], idxi)

                varianceVariableFactor[idxi] = 1 / weightDirect[row]
                meanVariableFactor[idxi] = meanDirect[row] * varianceVariableFactor[idxi]
                idxi += 1
            end
        end
    end

    ### Message send indices
    links = collect(1:idxi - 1)
    sendToFactor = sparse(toFactor, fromVariable, links, Nfactor, Nvariable)
    sendToVariable = sparse(toVariable, fromFactor, links, Nvariable, Nfactor)

    ### Set damping parameters
    alphaNew = fill(1.0, Nlink)
    alphaOld = fill(0.0, Nlink)
    bernoulliSample = randsubseq(collect(1:Nlink), dampProbability)
    @inbounds for i in bernoulliSample
        alphaNew[i] = 1.0 - dampAlpha
        alphaOld[i] = dampAlpha
    end

    ### Initialize marginal mean and variance vectors
    mean = fill(0.0, Nvariable)
    variance = fill(0.0, Nvariable)

    ### Iteration counters
    iterateFactor = collect(1:Nindirect)
    iterateVariable = collect(1:Nvariable)
    iterateMarginal = copy(iterateVariable)

    return ContinuousGraph(Nvariable, Nfactor, Nindirect, Nlink,
            meanVirtual, varianceVirtual, meanDirect, weightDirect, meanIndirect, varianceIndirect, coefficient,
            sendToVariable, sendToFactor, rowptr, colptr, colptrMarginal, alphaNew, alphaOld,
            iterateFactor, iterateVariable, iterateMarginal, dynamic),
           ContinuousInference(fromFactor, toVariable, meanFactorVariable, varianceFactorVariable,
            fromVariable, toFactor, meanVariableFactor, varianceVariableFactor, mean, variance)
end

########## Produce the graphical model ##########
function makeContinuousTreeGraph(system, virtualMean, virtualVariance, root)
    ### Number of factor and variable nodes
    Nfactor, Nvariable = size(system.jacobian)

    ### Find graph numbers, set the direct mean and variance, set factor numeration and pass through the rows
    Nlink = 0; Nindirect = 0
    meanDirect = fill(virtualMean / virtualVariance, Nvariable)
    weightDirect = fill(1 / virtualVariance, Nvariable)

    rowForward = [Int[] for i = 1:Nfactor]; rowBackward = [Int[] for i = 1:Nfactor]
    incomingToVariable = [Int[] for i = 1:Nvariable]
    @inbounds for i = 1:Nfactor
        NvariableInRow = system.jacobianTranspose.colptr[i + 1] - system.jacobianTranspose.colptr[i]
        if NvariableInRow == 1
            variable = system.jacobianTranspose.rowval[system.jacobianTranspose.colptr[i]]
            meanDirect[variable] = 0.0
            weightDirect[variable] = 0.0
        else
            Nlink += NvariableInRow
            Nindirect += 1
            for j = system.jacobianTranspose.colptr[i]:(system.jacobianTranspose.colptr[i + 1] - 1)
                row = system.jacobianTranspose.rowval[j]
                push!(rowForward[i], row)
                push!(rowBackward[i], row)
            end
        end
    end

    ## Pass through the columns
    iterateVariable = fill(0, Nvariable)
    counter = 0
    colForward = [Int[] for i = 1:Nvariable]; colBackward = [Int[] for i = 1:Nvariable]
    incomingToFactor = [Int[] for i = 1:Nfactor]
    @inbounds for col = 1:Nvariable
        for i = system.jacobian.colptr[col]:(system.jacobian.colptr[col + 1] - 1)
            row = system.jacobian.rowval[i]
            NvariableInRow = system.jacobianTranspose.colptr[row + 1] - system.jacobianTranspose.colptr[row]
            if NvariableInRow == 1
                meanDirect[col] += system.observation[row] * system.jacobian[row, col] / system.variance[row]
                weightDirect[col] += system.jacobian[row, col]^2 / system.variance[row]
            else
                push!(colForward[col], row)
                push!(colBackward[col], row)
            end
        end
        if length(colForward[col]) == 1 && col != root
            counter += 1
            iterateVariable[counter] = col
        end
    end
    resize!(iterateVariable, counter)

    ### Initialize data
    fromVariable = fill(0, Nlink)
    toFactor = fill(0, Nlink)
    meanVariableFactor = fill(0.0, Nlink)
    varianceVariableFactor = fill(0.0, Nlink)

    fromFactor = fill(0, Nlink)
    toVariable = fill(0, Nlink)
    meanFactorVariable = fill(0.0, Nlink)
    varianceFactorVariable = fill(0.0, Nlink)

    mean = fill(0, Nvariable)
    variance = fill(0, Nvariable)

    passFactorVariable = 0
    passVariableFactor = 0
    iterateFactor = Int64[]

    return ContinuousTreeGraph(Nvariable, Nfactor, Nlink, root, virtualMean, virtualVariance, meanDirect, weightDirect,
            rowForward, rowBackward, colForward, colBackward, incomingToFactor, incomingToVariable,
            iterateFactor, iterateVariable, passFactorVariable, passVariableFactor, true, true),
           ContinuousInference(fromFactor, toVariable, meanFactorVariable, varianceFactorVariable,
            fromVariable, toFactor, meanVariableFactor, varianceVariableFactor, mean, variance)
end


text:
module OneHotArrays

using Adapt
using ChainRulesCore
using CUDA
using LinearAlgebra
using MLUtils
using NNlib

export onehot, onehotbatch, onecold, OneHotArray,
  OneHotVector, OneHotMatrix, OneHotLike

include("array.jl")
include("onehot.jl")
include("linalg.jl")

end


text:
#
# Plot the results of a SIRD simulation
#

using Plots

function plotfit( md_input :: MDInput, md_traj :: MDTraj,
                  sird_input :: SIRDInput, sird_traj :: SIRDTraj,
                  filename :: String;
                  size=(400,400) )

  ENV["GKSwstype"]="nul"

  time = sird_traj.time
  sU = sird_traj.U
  sS = sird_traj.S
  sD = sird_traj.D
  sI = sird_traj.I
  nsteps = sird_traj.nsteps

  mU = md_traj.U
  mS = md_traj.S
  mD = md_traj.D
  mI = md_traj.I

  xmin = time[1]
  xmax = time[nsteps]

  plot(size=(400,400),framestyle=:box)

  scatter!(time,mU,label="",color=:blue     , marker = (:circle, 3, 1.0, :blue, stroke(0.1, 0.5, :black)))
  scatter!(time,mS,label="",color=:red      , marker = (:circle, 3, 1.0, :red, stroke(0.1, 0.5, :black)))
  scatter!(time,mI,label="",color=:darkgreen, marker = (:circle, 3, 1.0, :darkgreen, stroke(0.1, 0.5, :black)))
  scatter!(time,mD,label="",color=:black    , marker = (:circle, 3, 1.0, :black, stroke(0.1, 0.5, :black)))

  plot!(time,sU,label="",linewidth=2,color=:blue)
  plot!(time,sS,label="",linewidth=2,color=:red)
  plot!(time,sI,label="",linewidth=2,color=:darkgreen)
  plot!(time,sD,label="",linewidth=2,color=:black)

  plot!(xlabel="time")
  plot!(ylabel="Fraction of the population")
  plot!(xlim=[0,xmax],ylim=[0,1])

  plot!(rectangle(-0.01*xmax,0.4*xmax,0.76,1.02), opacity=0.9,label="",color=:white)

  fontsize=8
  x = 0.01*(xmax-xmin) + xmin
  d = 0.05 ; y = [ 0.99 + d ]
  annotate!(x,yd!(y,d),text("Final data:",:left,fontsize,:serif,:black))
  annotate!(x,yd!(y,d),text("Susceptible: $(@sprintf("%2.1f",sU[nsteps]*100))%",:left,fontsize,:serif,:blue))
  annotate!(x,yd!(y,d),text("Sick: $(@sprintf("%2.1f",sS[nsteps]*100))%",:left,fontsize,:serif,:red))
  annotate!(x,yd!(y,d),text("Immune: $(@sprintf("%2.1f",sI[nsteps]*100))%",:left,fontsize,:serif,:darkgreen))
  annotate!(x,yd!(y,d),text("Dead: $(@sprintf("%2.1f",sD[nsteps]*100))%",:left,fontsize,:serif,:black))

  plot!(rectangle(0.44*xmax, 1.01*xmax, 0.915, 1.02), opacity=0.9,label="",color=:white)

  x = 0.98*(xmax-xmin)+xmin
  d = 0.05 ; y = [ 0.99 + d ]
  s = @sprintf("kc = %7.5f ; kd = %7.5f",sird_input.kc,sird_input.kd)
  annotate!(x,yd!(y,d),text(s,:right,fontsize,:serif,:black))
  s = @sprintf("ki = %7.5f ; kiu = %7.5f",sird_input.ki,sird_input.kiu)
  annotate!(x,yd!(y,d),text(s,:right,fontsize,:serif,:black))

  savefig(filename)

end


text:
using ..CUDASupport
using ..TopOpt: @init_cuda
@init_cuda()
using ..GPUUtils
import ..TopOpt: whichdevice
using ..TopOpt: GPU

### Compliance

whichdevice(c::Compliance) = whichdevice(c.cell_comp)
whichdevice(c::MeanCompliance) = whichdevice(c.compliance)

hutch_rand!(v::CuArray) = v .= round.(CuArrays.CURAND.curand(size(v))) .* 2 .- 1

function Compliance(
    ::GPU,
    problem::StiffnessTopOptProblem{dim,T},
    solver::AbstractDisplacementSolver,
    ::Type{TI}=Int;
    tracing=false,
    logarithm=false,
    maxfevals=10^8,
) where {dim,T,TI}
    comp = T(0)
    cell_comp = zeros(CuVector{T}, getncells(problem.ch.dh.grid))
    grad = CuVector(
        fill(T(NaN), length(cell_comp) - sum(problem.black) - sum(problem.white))
    )
    topopt_trace = TopOptTrace{T,TI}()
    reuse = false
    fevals = TI(0)
    return Compliance(
        problem,
        solver,
        comp,
        cell_comp,
        grad,
        tracing,
        topopt_trace,
        reuse,
        fevals,
        logarithm,
        maxfevals,
    )
end
@define_cu(Compliance, :solver, :cell_comp, :grad, :cheqfilter)

function compute_compliance(
    cell_comp::CuVector{T}, grad, cell_dofs, Kes, u, black, white, varind, x, penalty, xmin
) where {T}
    args = (cell_comp, grad, cell_dofs, Kes, u, black, white, varind, x, penalty, xmin)
    callkernel(dev, comp_kernel1, args)
    CUDAdrv.synchronize(ctx)
    obj = compute_obj(cell_comp, x, varind, black, white, penalty, xmin)

    return obj
end

# CUDA kernels
function comp_kernel1(
    cell_comp::AbstractVector{T},
    grad,
    cell_dofs,
    Kes,
    u,
    black,
    white,
    varind,
    x,
    penalty,
    xmin,
) where {T}
    i = @thread_global_index()
    offset = @total_threads()
    @inbounds while i <= length(cell_comp)
        cell_comp[i] = zero(T)
        Ke = rawmatrix(Kes[i])
        for w in 1:size(Ke, 2)
            for v in 1:size(Ke, 1)
                if Ke isa Symmetric
                    cell_comp[i] += u[cell_dofs[v, i]] * Ke.data[v, w] * u[cell_dofs[w, i]]
                else
                    cell_comp[i] += u[cell_dofs[v, i]] * Ke[v, w] * u[cell_dofs[w, i]]
                end
            end
        end
        if !(black[i] || white[i])
            d = ForwardDiff.Dual{T}(x[varind[i]], one(T))
            if PENALTY_BEFORE_INTERPOLATION
                p = density(penalty(d), xmin)
            else
                p = penalty(density(d, xmin))
            end
            grad[varind[i]] = -p.partials[1] * cell_comp[i]
        end

        i += offset
    end
    return nothing
end

function compute_obj(
    cell_comp::AbstractVector{T},
    x,
    varind,
    black,
    white,
    penalty,
    xmin,
    ::Val{blocksize}=Val(80),
    ::Val{threads}=Val(256),
) where {T,blocksize,threads}
    result = similar(cell_comp, T, (blocksize,))
    args = (result, cell_comp, x, varind, black, white, penalty, xmin, Val(threads))
    @cuda blocks = blocksize threads = threads comp_kernel2(args...)
    CUDAnative.synchronize()
    obj = reduce(+, Array(result))
    return obj
end

function comp_kernel2(
    result,
    cell_comp::AbstractVector{T},
    x,
    varind,
    black,
    white,
    penalty,
    xmin,
    ::Val{LMEM},
) where {T,LMEM}
    @mapreduce_block(
        i,
        length(cell_comp),
        +,
        T,
        LMEM,
        result,
        begin
            w_comp(cell_comp[i], x[varind[i]], black[i], white[i], penalty, xmin)
        end
    )

    return nothing
end
@inline function w_comp(comp::T, x, black, white, penalty, xmin) where {T}
    if PENALTY_BEFORE_INTERPOLATION
        return ifelse(
            black, comp, ifelse(white, xmin * comp, (d = ForwardDiff.Dual{T}(x, one(T));
            p = density(penalty(d), xmin);
            p.value * comp))
        )
    else
        return ifelse(
            black, comp, ifelse(white, xmin * comp, (d = ForwardDiff.Dual{T}(x, one(T));
            p = penalty(density(d, xmin));
            p.value * comp))
        )
    end
end

### Volume

whichdevice(v::Volume) = whichdevice(v.cellvolumes)
@define_cu(Volume, :cellvolumes, :grad, :problem) # should be optimized to avoid replicating problem

function compute_volume(
    cellvolumes::CuVector{T},
    x,
    fixed_volume,
    varind,
    black,
    white,
    ::Val{blocksize}=Val(80),
    ::Val{threads}=Val(256),
) where {T,blocksize,threads}
    result = similar(cellvolumes, T, (blocksize,))
    args = (result, cellvolumes, x, varind, black, white, Val(threads))
    @cuda blocks = blocksize threads = threads volume_kernel(args...)
    CUDAnative.synchronize()
    vol = reduce(+, Array(result)) + fixed_volume
    return vol
end

function volume_kernel(
    result, cellvolumes::AbstractVector{T}, x, varind, black, white, ::Val{LMEM}
) where {T,LMEM}
    @mapreduce_block(
        i,
        length(cellvolumes),
        +,
        T,
        LMEM,
        result,
        begin
            if !(black[i]) && !(white[i])
                x[varind[i]] * cellvolumes[i]
            else
                zero(T)
            end
        end
    )

    return nothing
end


text:
# hash
hash(c::AbstractGray, hx::UInt) = hash(gray(c), hx)
hash(c::TransparentGray, hx::UInt) = hash(alpha(c), hash(gray(c), hx))
hash(c::AbstractRGB, hx::UInt) = hash(blue(c), hash(green(c), hash(red(c), hx)))
hash(c::TransparentRGB, hx::UInt) = hash(alpha(c), hash(blue(c), hash(green(c), hash(red(c), hx))))

Base.adjoint(c::Colorant) = c

# gamut{min,max}
gamutmax(::Type{T}) where {T<:HSV} = (360,1,1)
gamutmin(::Type{T}) where {T<:HSV} = (0,0,0)
gamutmax(::Type{T}) where {T<:HSL} = (360,1,1)
gamutmin(::Type{T}) where {T<:HSL} = (0,0,0)
gamutmax(::Type{T}) where {T<:Lab} = (100,128,128)
gamutmin(::Type{T}) where {T<:Lab} = (0,-127,-127)
gamutmax(::Type{T}) where {T<:LCHab} = (100,1,360)
gamutmin(::Type{T}) where {T<:LCHab} = (0,0,0)
gamutmax(::Type{T}) where {T<:YIQ} = (1,0.5226,0.5226)
gamutmin(::Type{T}) where {T<:YIQ} = (0,-0.5957,-0.5957)

gamutmax(::Type{T}) where {T<:AbstractGray} = (1,)
gamutmax(::Type{T}) where {T<:TransparentGray} = (1,1)
gamutmax(::Type{T}) where {T<:AbstractRGB} = (1,1,1)
gamutmax(::Type{T}) where {T<:TransparentRGB} = (1,1,1,1)
gamutmin(::Type{T}) where {T<:AbstractGray} = (0,)
gamutmin(::Type{T}) where {T<:TransparentGray} = (0,0)
gamutmin(::Type{T}) where {T<:AbstractRGB} = (0,0,0)
gamutmin(::Type{T}) where {T<:TransparentRGB} = (0,0,0,0)

# rand
for t in [Float16,Float32,Float64,N0f8,N0f16,N0f32]
    @eval _rand(::Type{T}) where {T<:Union{AbstractRGB{$t},AbstractGray{$t}}} =
          mapc(x->rand(eltype(T)), base_colorant_type(T)())
end

function _rand(::Type{T}) where T<:Colorant
    Gmax = gamutmax(T)
    Gmin = gamutmin(T)
    Mi = eltype(T) <: FixedPoint ? 1.0/typemax(eltype(T)) : 1.0
    A = rand(eltype(T), length(T))
    for j in eachindex(Gmax)
        A[j] = A[j] * (Mi * (Gmax[j]-Gmin[j])) + Gmin[j]
    end
    T(A...)
end

if Base.VERSION < v"1.6.0-DEV.1083"
    reinterpretc(::Type{T}, A::Array{C}) where {T<:Real,C<:AbstractGray} = reinterpret(T, A)
    reinterpretc(::Type{T}, A::Array{C}) where {T<:Real,C<:Colorant} = reshape(reinterpret(T, A), (sizeof(C)÷sizeof(T), size(A)...))
else
    reinterpretc(::Type{T}, A::Array{C}) where {T<:Real,C<:Colorant} = reinterpret(reshape, T, A)
end

function _rand(::Type{T}, sz::Dims) where T<:Colorant
    Mi = eltype(T) <: FixedPoint ? 1.0/typemax(eltype(T)) : 1.0
    A = Array{T}(undef, sz)
    Tr = eltype(T) <: FixedPoint ? FixedPointNumbers.rawtype(eltype(T)) : eltype(T)
    nchannels = sizeof(T)÷sizeof(eltype(T))
    Au = Random.UnsafeView(convert(Ptr{Tr}, pointer(A)), length(A)*nchannels)
    rand!(Au)

    Gmax = gamutmax(T)
    Gmin = gamutmin(T)
    s = (Gmax .- Gmin) .* Mi
    if !all(iszero, Gmin) || !all(==(1), s)
        Ar = reinterpretc(eltype(T), A)
        if T<:AbstractGray
            s1, offset1 = s[1], Gmin[1]
            for I in CartesianIndices(A)
                Ar[I] = Ar[I] * s1 + offset1
            end
        else
            for I in CartesianIndices(A)
                for j in eachindex(s)
                    Ar[j,I] = Ar[j,I] * s[j] + Gmin[j]
                end
            end
        end
    end
    return A
end

rand(::Type{T}, sz::Dims...) where {T<:Colorant} = _rand(ccolor(T, base_colorant_type(T){Float64}), sz...)

rand(::Type{Gray24}) = Gray24(rand(N0f8))
rand(::Type{Gray24}, sz::Dims) = Gray24.(rand(N0f8,sz))
rand(::Type{AGray32}) = AGray32(rand(N0f8),rand(N0f8))
rand(::Type{AGray32}, sz::Dims) = AGray32.(rand(N0f8,sz),rand(N0f8,sz))

# broadcast

# Mapping a function over color channels
"""
    mapc(f, rgb) -> rgbf
    mapc(f, rgb1, rgb2) -> rgbf

`mapc` applies the function `f` to each color channel of the input
color(s), returning an output color in the same colorspace.

# Examples:

    julia> mapc(x->clamp(x,0,1), RGB(-0.2,0.3,1.2))
    RGB{Float64}(0.0,0.3,1.0)

    julia> mapc(max, RGB(0.1,0.8,0.3), RGB(0.5,0.5,0.5))
    RGB{Float64}(0.5,0.8,0.5)

    julia> mapc(+, RGB(0.1,0.8,0.3), RGB(0.5,0.5,0.5))
    RGB{Float64}(0.6,1.3,0.8)
"""
@inline mapc(f, c::C) where {C<:ColorantN{1}} =
    base_colorant_type(C)(f(comp1(c)))
@inline mapc(f, c::C) where {C<:ColorantN{2}} =
    base_colorant_type(C)(f(comp1(c)), f(comp2(c)))
@inline mapc(f, c::C) where {C<:ColorantN{3}} =
    base_colorant_type(C)(f(comp1(c)), f(comp2(c)), f(comp3(c)))
@inline mapc(f, c::C) where {C<:ColorantN{4}} =
    base_colorant_type(C)(f(comp1(c)), f(comp2(c)), f(comp3(c)), f(comp4(c)))
@inline mapc(f, c::C) where {C<:ColorantN{5}} =
    base_colorant_type(C)(f(comp1(c)), f(comp2(c)), f(comp3(c)), f(comp4(c)), f(comp5(c)))

@inline mapc(f, x::Number) = f(x)

mapc(f::F, x, y) where F = _mapc(_same_colorspace(x,y), f, x, y)
_mapc(::Type{C}, f, x::ColorantN{1}, y::ColorantN{1}) where C =
    C(f(comp1(x), comp1(y)))
_mapc(::Type{C}, f, x::ColorantN{2}, y::ColorantN{2}) where C =
    C(f(comp1(x), comp1(y)), f(comp2(x), comp2(y)))
_mapc(::Type{C}, f, x::ColorantN{3}, y::ColorantN{3}) where C =
    C(f(comp1(x), comp1(y)), f(comp2(x), comp2(y)), f(comp3(x), comp3(y)))
_mapc(::Type{C}, f, x::ColorantN{4}, y::ColorantN{4}) where C =
    C(f(comp1(x), comp1(y)), f(comp2(x), comp2(y)), f(comp3(x), comp3(y)), f(comp4(x), comp4(y)))
_mapc(::Type{C}, f, x::ColorantN{5}, y::ColorantN{5}) where C =
    C(f(comp1(x), comp1(y)), f(comp2(x), comp2(y)), f(comp3(x), comp3(y)), f(comp4(x), comp4(y)), f(comp5(x), comp5(y)))

_same_colorspace(x::Colorant, y::Colorant) = _same_colorspace(base_colorant_type(x),
                                                              base_colorant_type(y))
_same_colorspace(::Type{C}, ::Type{C}) where {C<:Colorant} = C
@noinline _same_colorspace(::Type{C1}, ::Type{C2}) where {C1<:Colorant,C2<:Colorant} =
    throw(ArgumentError("$C1 and $C2 are from different colorspaces"))

"""
    reducec(op, v0, c)

Reduce across color channels of `c` with the binary operator
`op`. `v0` is the neutral element used to initiate the reduction. For
grayscale,

    reducec(op, v0, c::Gray) = op(v0, comp1(c))

whereas for RGB

    reducec(op, v0, c::RGB) = op(comp3(c), op(comp2(c), op(v0, comp1(c))))

If `c` has an alpha channel, it is always the last one to be folded into the reduction.
"""
@inline reducec(op, v0, c::ColorantN{1}) = op(v0, comp1(c))
@inline reducec(op, v0, c::ColorantN{2}) = op(comp2(c), op(v0, comp1(c)))
@inline reducec(op, v0, c::ColorantN{3}) = op(comp3(c), op(comp2(c), op(v0, comp1(c))))
@inline reducec(op, v0, c::ColorantN{4}) = op(comp4(c), op(comp3(c), op(comp2(c), op(v0, comp1(c)))))
@inline reducec(op, v0, c::ColorantN{5}) = op(comp5(c), op(comp4(c), op(comp3(c), op(comp2(c), op(v0, comp1(c))))))

@inline reducec(op, v0, x::Number) = op(v0, x)

"""
    mapreducec(f, op, v0, c)

Reduce across color channels of `c` with the binary operator `op`,
first applying `f` to each channel. `v0` is the neutral element used
to initiate the reduction. For grayscale,

    mapreducec(f, op, v0, c::Gray) = op(v0, f(comp1(c)))

whereas for RGB

    mapreducec(f, op, v0, c::RGB) = op(f(comp3(c)), op(f(comp2(c)), op(v0, f(comp1(c)))))

If `c` has an alpha channel, it is always the last one to be folded into the reduction.
"""
@inline mapreducec(f, op, v0, c::ColorantN{1}) = op(v0, f(comp1(c)))
@inline mapreducec(f, op, v0, c::ColorantN{2}) = op(f(comp2(c)), op(v0, f(comp1(c))))
@inline mapreducec(f, op, v0, c::ColorantN{3}) =
    op(f(comp3(c)), op(f(comp2(c)), op(v0, f(comp1(c)))))
@inline mapreducec(f, op, v0, c::ColorantN{4}) =
    op(f(comp4(c)), op(f(comp3(c)), op(f(comp2(c)), op(v0, f(comp1(c))))))
@inline mapreducec(f, op, v0, c::ColorantN{5}) =
    op(f(comp5(c)), op(f(comp4(c)), op(f(comp3(c)), op(f(comp2(c)), op(v0, f(comp1(c)))))))

@inline mapreducec(f, op, v0, x::Number) = op(v0, f(x))


text:
module HelpersTests

using Test
using Gridap.Helpers

end # module HelpersTests


text:
module Strapping

using Tables, StructTypes

struct Error <: Exception
    msg::String
end

"""
    Strapping.construct(T, tbl)
    Strapping.construct(Vector{T}, tbl)

Given a [Tables.jl](https://github.com/JuliaData/Tables.jl)-compatible input table source `tbl`,
construct an instance of `T` (single object, first method), or `Vector{T}` (list of objects, 2nd method).

The 1st method will throw an error if the input table is empty, and warn if there are more rows
than necessary to construct a single `T`.

The 2nd method will return an empty list for an empty input source, and construct as many `T` as are
found until the input table is exhausted.

`Strapping.construct` utilizes the [StructTypes.jl](https://github.com/JuliaData/StructTypes.jl) package
for determining the `StructTypes.StructType` trait of `T` and constructing an instance appropriately:
  * `StructTypes.Struct`/`StructTypes.Mutable`: field reflection will be used to retrieve values from the input table row, with field customizations respected, like excluded fields, field-specific keyword args, etc.
  * `StructTypes.DictType`: each column name/value of the table row will be used as a key/value pair to be passed to the `DictType` constructor
  * `StructTypes.ArrayType`: column values will be "collected" as an array to be passed to the `ArrayType` constructor
  * `StructTypes.StringType`/`StructTypes.NumberType`/`StructTypes.BoolType`/`StructTypes.NullType`: only the first value of the row will be passed to the scalar type constructor

Note that for `StructTypes.DictType` and `StructTypes.ArrayType`, "aggregate" value/eltypes are not allowed, since
the entire row is treated as key/value pairs or array elements. That means, for example, I can't have a table with rows like `tbl = [(a=1, b=2)]` and try to do `Strapping.construct(Dict{Symbol, Dict{Int, Int}}, tbl)`. It first attempts to map column names to the outer `Dict` keys, (`a` and `b`), but then tries to map the values `1` and `2`
to `Dict{Int, Int}` and fails.

For structs with `ArrayType` fields, the first row values will be used for other scalar fields, and subsequent rows
will be iterated for the `ArrayType` field values. For example, I may wish to construct a type like:
```julia
struct TestResult
    id::Int
    values::Vector{Float64}
end
StructTypes.StructType(::Type{TestResult}) = StructTypes.Struct()
StructTypes.idproperty(::Type{TestResult}) = :id
```

and my input table would look something like, `tbl = (id=[1, 1, 1], values=[3.14, 3.15, 3.16])`. I can then construct my type like:

```julia
julia> Strapping.construct(TestResult, tbl)
TestResult(1, [3.14, 3.15, 3.16])
```

Note that along with defining the `StructTypes.StructType` trait for `TestResult`, I also needed to define `StructTypes.idproperty` to signal which field of my struct is a "unique key" identifier. This enables Strapping to distinguish which rows belong to a particular instance of `TestResult`. This allows the slightly more complicated example of returning multiple `TestResult`s from a single table:

```julia
julia> tbl = (id=[1, 1, 1, 2, 2, 2], values=[3.14, 3.15, 3.16, 40.1, 0.01, 2.34])
(id = [1, 1, 1, 2, 2, 2], values = [3.14, 3.15, 3.16, 40.1, 0.01, 2.34])

julia> Strapping.construct(Vector{TestResult}, tbl)
2-element Array{TestResult,1}:
 TestResult(1, [3.14, 3.15, 3.16])
 TestResult(2, [40.1, 0.01, 2.34])
```

Here, we actually have _two_ `TestResult` objects in our `tbl`, and Strapping uses the `id` field to identify object owners for a row. Note that currently the table rows need to be sorted on the `idproperty` field, i.e. rows belonging to the same object must appear consecutively in the input table rows.

Now let's discuss "aggregate" type fields. Let's say I have a struct like:

```julia
struct Experiment
    id::Int
    name::String
    testresults::TestResult
end
StructTypes.StructType(::Type{Experiment}) = StructTypes.Struct()
StructTypes.idproperty(::Type{Experiment}) = :id
```

So my `Experiment` type also as an `id` field, in addition to a `name` field, and an "aggregate" field of `testresults`. How should the input table source account for `testresults`, which is itself a struct made up of its own `id` and `values` fields? The key here is "flattening" nested structs into a single set of table column names, and utilizing the `StructTypes.fieldprefix` function, which allows specifying a `Symbol` prefix to identify an aggregate field's columns in the table row. So, in the case of our `Experiment`, we can do:

```julia
StructTypes.fieldprefix(::Type{Experiment}, nm::Symbol) = nm == :testresults ? :testresults_ : :_
```

Note that this is the default definition, so we don't really need to define this, but for illustration purposes, we'll walk through it. We're saying that for the `:testresults` field name, we should expect its column names in the table row to start with `:testresults_`. So the table data for an `Experiment` instance, would look something like:

```julia
tbl = (id=[1, 1, 1], name=["exp1", "exp1", "exp1"], testresults_id=[1, 1, 1], testresults_values=[3.14, 3.15, 3.16])
```

This pattern generalizes to structs with multiple aggregate fields, or aggregate fields that themselves have aggregate fields (nested aggregates); in the nested case, the prefixes are concatenated, like `testresults_specifictestresult_id`.
"""
function construct end

function construct(::Type{T}, source; silencewarnings::Bool=false, kw...) where {T}
    rows = Tables.rows(source)
    state = iterate(rows)
    state === nothing && throw(Error("can't construct `$T` from empty source"))
    row, st = state
    x, state = construct(rows, st, row, T; kw...)
    state === nothing || (silencewarnings && println("warning: additional source rows left after reading `$T`"))
    return x
end

function construct(::Type{Vector{T}}, source; kw...) where {T}
    rows = Tables.rows(source)
    state = iterate(rows)
    A = Vector{T}(undef, 0)
    while state !== nothing
        row, st = state
        x, state = construct(rows, st, row, T; kw...)
        push!(A, x)
    end
    return A
end

# collection handler: will iterate results for collection fields
function construct(rows, st, row, ::Type{T}; kw...) where {T}
    x = construct(StructTypes.StructType(T), row, T; kw...)
    idprop = StructTypes.idproperty(T)
    if idprop !== :_
        id = Tables.getcolumn(row, idprop)
        state = iterate(rows, st)
        if state !== nothing
            while state !== nothing
                row, st = state
                Tables.getcolumn(row, idprop) == id || break
                construct!(StructTypes.StructType(T), row, x; kw...)
                state = iterate(rows, st)
            end
        end
    else
        state = iterate(rows, st)
    end
    return x, state
end

# aggregate handlers (don't take specific `col`/`nm` arguments)
# construct versions construct initial object
# construct! versions take existing object and append additional elements to collection fields
function construct(::StructTypes.Struct, row, ::Type{T}; kw...) where {T}
    coloffset = Ref{Int}(0)
    return StructTypes.construct(T) do i, nm, TT
        construct(StructTypes.StructType(TT), T, row, i, coloffset, Symbol(), nm, TT; kw...)
    end
end

function construct!(::StructTypes.Struct, row, x::T; kw...) where {T}
    coloffset = Ref{Int}(0)
    return StructTypes.foreachfield(x) do i, nm, TT, v
        construct!(StructTypes.StructType(TT), T, row, i, coloffset, Symbol(), nm, TT, v; kw...)
    end
end

function construct(::StructTypes.Mutable, row, ::Type{T}; kw...) where {T}
    x = T()
    coloffset = Ref{Int}(0)
    StructTypes.mapfields!(x) do i, nm, TT
        y = construct(StructTypes.StructType(TT), T, row, i, coloffset, Symbol(), nm, TT; kw...)
        return y
    end
    return x
end

function construct!(::StructTypes.Mutable, row, x::T; kw...) where {T}
    coloffset = Ref{Int}(0)
    return StructTypes.foreachfield(x) do i, nm, TT, v
        construct!(StructTypes.StructType(TT), T, row, i, coloffset, Symbol(), nm, TT, v; kw...)
    end
end

# default aggregate
construct(::StructTypes.Struct, row, ::Type{Any}; kw...) =
    construct(StructTypes.DictType(), row, Dict{String, Any}; kw...)

construct(::StructTypes.DictType, row, ::Type{T}; kw...) where {T} = construct(StructTypes.DictType(), row, T, Symbol, Any; kw...)
construct(::StructTypes.DictType, row, ::Type{T}; kw...) where {T <: NamedTuple} = construct(StructTypes.DictType(), row, T, Symbol, Any; kw...)
construct(::StructTypes.DictType, row, ::Type{Dict}; kw...) = construct(StructTypes.DictType(), row, Dict, String, Any; kw...)
construct(::StructTypes.DictType, row, ::Type{T}; kw...) where {T <: AbstractDict} = construct(StructTypes.DictType(), row, T, keytype(T), valtype(T); kw...)

function construct(::StructTypes.DictType, row, ::Type{T}, ::Type{K}, ::Type{V}; kw...) where {T, K, V}
    #TODO: formally disallow aggregate types as V?
    # the problem is we're already treating the DictType as the aggregate, so
    # it's impossible to distinguish between multiple aggregate type V
    # from a single resultset set of fields
    x = Dict{K, V}()
    for (i, nm) in enumerate(Tables.columnnames(row))
        val = construct(StructTypes.StructType(V), T, row, i, Ref{Int}(0), Symbol(), nm, V; kw...)
        if K == Symbol
            x[nm] = val
        else
            x[StructTypes.construct(K, String(nm))] = val
        end
    end
    return StructTypes.construct(T, x; kw...)
end

function construct!(::StructTypes.DictType, row, x::T; kw...) where {T}
    for (i, nm) in enumerate(Tables.columnnames(row))
        v = x[nm]
        V = typeof(v)
        construct!(StructTypes.StructType(V), T, row, i, Ref{Int}(0), Symbol(), nm, V, v; kw...)
    end
    return
end

construct(::StructTypes.ArrayType, row, ::Type{T}; kw...) where {T} = construct(StructTypes.ArrayType(), row, T, Base.IteratorEltype(T) == Base.HasEltype() ? eltype(T) : Any; kw...)
construct(::StructTypes.ArrayType, row, ::Type{T}, ::Type{eT}; kw...) where {T, eT} = constructarray(row, T, eT; kw...)
construct(::StructTypes.ArrayType, row, ::Type{Tuple}, ::Type{eT}; kw...) where {eT} = constructarray(row, Tuple, eT; kw...)

function constructarray(row, ::Type{T}, ::Type{eT}; kw...) where {T, eT}
    #TODO: disallow aggregate eT? same problem as DictType; we're treating the ArrayType as our aggregate
    nms = Tables.columnnames(row)
    N = length(nms)
    x = Vector{eT}(undef, N)
    for (i, nm) in enumerate(nms)
        x[i] = construct(StructTypes.StructType(eT), T, row, i, Ref{Int}(0), Symbol(), nm, eT; kw...)
    end
    return StructTypes.construct(T, x; kw...)
end

function construct!(::StructTypes.ArrayType, row, x::T; kw...) where {T}
    for (i, nm) in enumerate(Tables.columnnames(row))
        v = x[i]
        V = typeof(v)
        construct!(StructTypes.StructType(V), T, row, i, Ref{Int}(0), Symbol(), nm, V, v; kw...)
    end
    return
end

function construct(::StructTypes.ArrayType, row, ::Type{T}, ::Type{eT}; kw...) where {T <: Tuple, eT}
    return StructTypes.construct(T) do i, nm, TT
        construct(StructTypes.StructType(TT), T, row, i, Ref{Int}(0), Symbol(), nm, TT; kw...)
    end
end

# constructing a single scalar from a row
# for example, you want the result to be a single Int from: SELECT COUNT(*) FROM table
function construct(::StructTypes.StringType, row, ::Type{T}; kw...) where {T}
    return StructTypes.construct(T, Tables.getcolumn(row, 1))
end

function construct(::StructTypes.NumberType, row, ::Type{T}; kw...) where {T}
    return StructTypes.construct(T, StructTypes.numbertype(T)(Tables.getcolumn(row, 1)))
end

function construct(::StructTypes.BoolType, row, ::Type{T}; kw...) where {T}
    return StructTypes.construct(T, Tables.getcolumn(row, 1))
end

function construct(::StructTypes.NullType, row, ::Type{T}; kw...) where {T}
    return StructTypes.construct(T, Tables.getcolumn(row, 1))
end

## field construction: here we have a specific col::Int/nm::Symbol argument for a parent aggregate
# Struct field
function construct(::StructTypes.Struct, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T}
    prefix = Symbol(prefix, StructTypes.fieldprefix(PT, nm))
    off = 0
    x = StructTypes.construct(T) do i, nm, TT
        off += 1
        construct(StructTypes.StructType(TT), T, row, coloffset[] + col + i - 1, coloffset, prefix, nm, TT; kw...)
    end
    coloffset[] += off - 1
    return x
end

function construct!(::Union{StructTypes.Struct, StructTypes.Mutable}, PT, row, col, coloffset, prefix, nm, ::Type{T}, v; kw...) where {T}
    off = 0
    StructTypes.foreachfield(v) do i, nm, TT, v
        off += 1
    end
    coloffset[] += off - 1
    return
end

function construct(::StructTypes.Mutable, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T}
    prefix = Symbol(prefix, StructTypes.fieldprefix(PT, nm))
    off = 0
    x = T()
    StructTypes.mapfields!(x) do i, nm, TT
        off += 1
        construct(StructTypes.StructType(TT), T, row, coloffset[] + col + i - 1, coloffset, prefix, nm, TT; kw...)
    end
    coloffset[] += off - 1
    return x
end

construct(::StructTypes.DictType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T} = construct(StructTypes.DictType(), PT, row, col, coloffset, prefix, nm, T, Symbol, Any; kw...)
construct(::StructTypes.DictType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T <: NamedTuple} = construct(StructTypes.DictType(), PT, row, col, coloffset, prefix, nm, T, Symbol, Any; kw...)
construct(::StructTypes.DictType, PT, row, col, coloffset, prefix, nm, ::Type{Dict}; kw...) = construct(StructTypes.DictType(), PT, row, col, coloffset, prefix, nm, Dict, String, Any; kw...)
construct(::StructTypes.DictType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T <: AbstractDict} = construct(StructTypes.DictType(), PT, row, col, coloffset, prefix, nm, T, keytype(T), valtype(T); kw...)

function construct(::StructTypes.DictType, PT, row, col, coloffset, prefix, nm, ::Type{T}, ::Type{K}, ::Type{V}; kw...) where {T, K, V}
    prefix = String(Symbol(prefix, StructTypes.fieldprefix(PT, nm)))
    off = 0
    x = Dict{K, V}()
    for (i, nm) in enumerate(Tables.columnnames(row))
        if startswith(String(nm), prefix)
            val = construct(StructTypes.StructType(V), T, row, i, Ref{Int}(0), Symbol(), nm, V; kw...)
            off += 1
            if K == Symbol
                x[nm] = val
            else
                x[StructTypes.construct(K, String(nm))] = val
            end
        end
    end
    coloffset[] += off
    return StructTypes.construct(T, x; kw...)
end

function construct(::StructTypes.ArrayType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T}
    eT = Base.IteratorEltype(T) == Base.HasEltype() ? eltype(T) : Any
    return StructTypes.construct(T, [construct(StructTypes.StructType(eT), PT, row, col, coloffset, prefix, nm, eT; kw...)]; kw...)
end

function construct!(::StructTypes.ArrayType, PT, row, col, coloffset, prefix, nm, ::Type{T}, v; kw...) where {T}
    eT = Base.IteratorEltype(T) == Base.HasEltype() ? eltype(T) : Any
    push!(v, construct(StructTypes.StructType(eT), PT, row, col, coloffset, prefix, nm, eT; kw...))
    return
end

# for all other construct!, we ignore
construct!(ST, PT, row, col, coloffset, prefix, nm, T, v; kw...) = nothing

# scalar handlers (take a `col` argument)
getvalue(row, ::Type{T}, col::Int, nm::Symbol) where {T} = Tables.getcolumn(row, T, col, nm)
getvalue(row, ::Type{T}, col::Int, nm::Int) where {T} = Tables.getcolumn(row, nm)

construct(::StructTypes.Struct, PT, row, col, coloffset, prefix, nm, ::Type{Any}; kw...) = getvalue(row, Any, col, Symbol(prefix, nm))
construct(::StructTypes.Struct, PT, row, col, coloffset, prefix, nm, U::Union; kw...) = getvalue(row, U, col, Symbol(prefix, nm))
construct(::StructTypes.StringType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T} = StructTypes.construct(T, getvalue(row, T, col, Symbol(prefix, nm)))
construct(::StructTypes.NumberType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T} =
    StructTypes.construct(T, StructTypes.numbertype(T)(getvalue(row, T, col, Symbol(prefix, nm))))
construct(::StructTypes.BoolType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T} = StructTypes.construct(T, getvalue(row, T, col, Symbol(prefix, nm)))
construct(::StructTypes.NullType, PT, row, col, coloffset, prefix, nm, ::Type{T}; kw...) where {T} = StructTypes.construct(T, getvalue(row, T, col, Symbol(prefix, nm)))

# deconstruct

"""
    Strapping.deconstruct(x::T)
    Strapping.deconstruct(x::Vector{T})

The inverse of `Strapping.construct`, where an object instance `x::T` or `Vector` of objects `x::Vector{T}` is "deconstructed" into a Tables.jl-compatible row iterator. This works following the same patterns outlined in `Strapping.construct` with regards to `ArrayType` and aggregate fields. Specifically, `ArrayType` fields will cause multiple rows to be outputted, one row per collection element, with other scalar fields being repeated in each row. Similarly for aggregate fields, the field prefix will be used (`StructTypes.fieldprefix`) and nested aggregates will all be flattened into a single list of column names with aggregate prefixes.

In general, this allows outputting any "object" as a 2D table structure that could be stored in any Tables.jl-compatible sink format, e.g. csv file, sqlite table, mysql database table, feather file, etc.
"""
function deconstruct end

# single object or vector => Tables.rows iterator
struct DeconstructedRowsIterator{T, A}
    values::A
    lens::Vector{Int}
    len::Int
    names::Vector{Symbol}
    types::Vector{Type}
    fieldindices::Vector{Tuple{Int, Int}}
    fieldnames::Vector{Tuple{Symbol, Symbol}}
    lookup::Dict{Symbol, Int}
end

Tables.isrowtable(::Type{<:DeconstructedRowsIterator}) = true

# disabled because
# Tables.schema(x::DeconstructedRowsIterator) = Tables.Schema(x.names, x.types)

struct DeconstructedRow{T} <: Tables.AbstractRow
    x::T # a single object we're deconstructing
    index::Int # index of this specific row (may be > 1 for objects w/ collection fields)
    names::Vector{Symbol}
    fieldindices::Vector{Tuple{Int, Int}}
    fieldnames::Vector{Tuple{Symbol, Symbol}}
    lookup::Dict{Symbol, Int}
end

obj(x::DeconstructedRow) = getfield(x, :x)
ind(x::DeconstructedRow) = getfield(x, :index)
names(x::DeconstructedRow) = getfield(x, :names)
inds(x::DeconstructedRow) = getfield(x, :fieldindices)
nms(x::DeconstructedRow) = getfield(x, :fieldnames)
lookup(x::DeconstructedRow) = getfield(x, :lookup)

Tables.columnnames(row::DeconstructedRow) = names(row)
Tables.getcolumn(row::DeconstructedRow, ::Type{T}, i::Int, nm::Symbol) where {T} =
    getfieldvalue(obj(row), ind(row), inds(row)[i], nms(row)[i])
Tables.getcolumn(row::DeconstructedRow, i::Int) =
    getfieldvalue(obj(row), ind(row), inds(row)[i], nms(row)[i])
Tables.getcolumn(row::DeconstructedRow, nm::Symbol) = Tables.getcolumn(row, lookup(row)[nm])

mutable struct DeconstructClosure{PT}
    len::Int
    names::Vector{Symbol}
    types::Vector{Type}
    fieldindices::Vector{Tuple{Int, Int}}
    fieldnames::Vector{Tuple{Symbol, Symbol}}
    parentindex::Int
    parentname::Symbol
    prefix::Symbol
    j::Int
    nm::Symbol
    parenttype::Type{PT}
end

DeconstructClosure(PT) = DeconstructClosure(1, Symbol[], Type[], Tuple{Int, Int}[], Tuple{Symbol, Symbol}[], 0, Symbol(), Symbol(), 1, Symbol(), PT)

function (f::DeconstructClosure)(i, nm, TT, v; kw...)
    len = valuelength(StructTypes.StructType(TT), v)
    if len > f.len
        f.len = len
    end
    nametypeindex!(v, i, nm, f)
    return
end

mutable struct DeconstructLenClosure
    len::Int
end

function (f::DeconstructLenClosure)(i, nm, TT, v; kw...)
    len = valuelength(StructTypes.StructType(TT), v)
    if len > f.len
        f.len = len
    end
    return
end

deconstructobj!(x::T, c) where {T} = deconstructobj!(StructTypes.StructType(T), x, c)
deconstructobj!(ST, x, c) = StructTypes.foreachfield(c, x)

function deconstructobj!(::StructTypes.DictType, x, c)
    i = 1
    for (k, v) in StructTypes.keyvaluepairs(x)
        c(i, k, typeof(v), v)
        i += 1
    end
    return
end

deconstruct(x) = deconstruct([x])

function deconstruct(values::A) where {T, A <: AbstractVector{T}}
    c = DeconstructClosure(T)
    flens = DeconstructLenClosure(0)
    lens = Int[]
    len = 0
    i = 0
    for x in values
        if i == 0
            deconstructobj!(x, c)
            push!(lens, c.len)
            len += c.len
        else
            flens.len = 0
            deconstructobj!(x, flens)
            len += flens.len
            push!(lens, flens.len)
        end
        i += 1
    end
    lookup = Dict(nm => i for (i, nm) in enumerate(c.names))
    return DeconstructedRowsIterator{T, A}(values, lens, len, c.names, c.types, c.fieldindices, c.fieldnames, lookup)
end

Base.eltype(x::DeconstructedRowsIterator{T}) where {T} = DeconstructedRow{T}
Base.length(x::DeconstructedRowsIterator) = x.len

# single object w/ plain scalar fields: 1 obj => 1 row
# single object w/ scalar field and collection field: 1 obj => N rows, where N = length(collection)
# single object w/ aggregate object field: 1 obj => 1 row
# single object w/ aggregate object collection field: 1 obj => N rows, where N = length(collection)
@inline function Base.iterate(x::DeconstructedRowsIterator, (i, j, k)=(1, 1, 1))
    k > x.len && return nothing
    if j == x.lens[i]
        return DeconstructedRow(x.values[i], j, x.names, x.fieldindices, x.fieldnames, x.lookup), (i + 1, 1, k + 1)
    else
        return DeconstructedRow(x.values[i], j, x.names, x.fieldindices, x.fieldnames, x.lookup), (i, j + 1, k + 1)
    end
end

valuelength(ST, x) = 1
valuelength(::StructTypes.ArrayType, x) = length(x)

getfieldvalue(x::T, ind, (i, j), (nm1, nm2)) where {T} = getfieldvalue(StructTypes.StructType(T), x, ind, (i, j), (nm1, nm2))
getfieldvalue(ST, x, ind, (i, j), (nm1, nm2)) = x
getfieldvalue(::StructTypes.ArrayType, x, ind, (i, j), (nm1, nm2)) = getfieldvalue(x[ind], ind, (i, j), (nm1, nm2))
getfieldvalue(::Union{StructTypes.Struct, StructTypes.Mutable}, x, ind, (i, j), (nm1, nm2)) = getsubfieldvalue(Core.getfield(x, i), ind, (i, j), (nm1, nm2))
getfieldvalue(::StructTypes.DictType, x, ind, (i, j), (nm1, nm2)) = getsubfieldvalue(x[nm1], ind, (i, j), (nm1, nm2))

getsubfieldvalue(x::T, ind, (i, j), (nm1, nm2)) where {T} = getsubfieldvalue(StructTypes.StructType(T), x, ind, (i, j), (nm1, nm2))
getsubfieldvalue(ST, x, ind, (i, j), (nm1, nm2)) = x
getsubfieldvalue(::StructTypes.ArrayType, x, ind, (i, j), (nm1, nm2)) = getsubfieldvalue(x[ind], ind, (i, j), (nm1, nm2))
getsubfieldvalue(::Union{StructTypes.Struct, StructTypes.Mutable}, x, ind, (i, j), (nm1, nm2)) = getsubfieldvalue(Core.getfield(x, j), ind, (i, j), (nm1, nm2))
function getsubfieldvalue(::StructTypes.DictType, x, ind, (i, j), (nm1, nm2))
    getsubfieldvalue(x[nm2], ind, (i, j), (nm1, nm2))
end

rowtypeof(x::T) where {T} = rowtypeof(StructTypes.StructType(T), x)
rowtypeof(::StructTypes.ArrayType, x) = eltype(x)
rowtypeof(ST, x) = typeof(x)

nametypeindex!(x::T, i, nm, c) where {T} = nametypeindex!(StructTypes.StructType(T), x, i, nm, c)

function nametypeindex!(ST, x, i, nm, c)
    push!(c.names, Symbol(c.prefix, nm))
    push!(c.types, rowtypeof(x))
    push!(c.fieldindices, (c.parentindex > 0 ? c.parentindex : i, c.j))
    push!(c.fieldnames, (c.parentindex > 0 ? c.parentname : nm, nm))
    c.j += 1
    return
end

function nametypeindex!(::Union{StructTypes.Struct, StructTypes.Mutable, StructTypes.DictType}, x::T, i, nm, c) where {T}
    c2 = DeconstructClosure(T)
    c2.names = c.names
    c2.types = c.types
    c2.fieldindices = c.fieldindices
    c2.fieldnames = c.fieldnames
    c2.parentindex = i
    c2.parentname = nm
    c2.prefix = Symbol(c.prefix, StructTypes.fieldprefix(c.parenttype, nm))
    deconstructobj!(x, c2)
    return
end

function nametypeindex!(::StructTypes.ArrayType, x::T, i, nm, c) where {T}
    nametypeindex!(x[1], i, nm, c)
end

end # module


text:
module Dualization

using JuMP, MathOptInterface
const MOI  = MathOptInterface
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

const SVF = MOI.SingleVariable
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}
const SQF{T} = MOI.ScalarQuadraticFunction{T}
const VQF{T} = MOI.VectorQuadraticFunction{T}

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

import MathOptInterface: dual_set

include("structures.jl")
include("utils.jl")
include("dual_sets.jl")
include("supported.jl")
include("dual_names.jl")
include("objective_coefficients.jl")
include("add_dual_cone_constraint.jl")
include("dual_model_variables.jl")
include("dual_equality_constraints.jl")
include("dualize.jl")
include("MOI_wrapper.jl")

end # module


text:
# https://juliabyexample.helpmanual.io/
# https://www.geeksforgeeks.org/taking-input-from-users-in-julia/

println("What's your name?")

inputName = readline()

if inputName == ""
	inputName = "Jakob"
end

println("Type in a number, ", inputName)

inputVar1 = readline()

println("Type in another number, ", inputName)

inputVar2 = readline()

println("Now add $inputVar1 and $inputVar2 = ")

inputSum = readline()

if parse(Int64, inputSum) == parse(Int64, inputVar1) + parse(Int64, inputVar2)
    println("Good job $inputName !")
else
	println("Incorrect answer $inputName !")
	println("Correct answer is: ", parse(Int64, inputVar1) + parse(Int64, inputVar2))
end


text:
using Test
using MRI, SWI

data_path = "testData/small"
TEs = [4,8,12]

mag_nii = readmag("$data_path/Mag.nii")
hdr = mag_nii.header
phase_nii = readphase("$data_path/Phase.nii")

data = Data(mag_nii, phase_nii, hdr, TEs)

swi = calculateSWI(data)

printstyled("Tests Passed!"; color=:green)


text:
# ==========================================================================================
# COMPARISON OF CLASSIFICATION RATES BETWEEN 3 DIFFERENT LOCAL DISCRIMINANT BASES (LDB)
# ALGORITHMS
#
# This is the source code used to compare the differences in classification rates between 3
# LDB algorithms:
# ==========================================================================================

using Wavelets,
      WaveletsExt,
      Statistics,
      DataFrames,
      CSV,
      MLJ

include("utils.jl")

## ========= Setup =========
# TODO: Change wavelet filter
wt = wavelet(WT.coif2)
# TODO: Change decomposition levels
L = 5
# TODO: Change the value of K
k = 10
# TODO: Change the number of features extracted
n_features = 10
# TODO： Change the type of discriminant power
dp = BasisDiscriminantMeasure()
# TODO: Change the discriminant measures
dm_ldbk    = AsymmetricRelativeEntropy()
dm_ldbkash = SymmetricRelativeEntropy()
dm_ldbkemd = EarthMoverDistance()
# Create LDB objects
ldbk    = LocalDiscriminantBasis(wt = wt,
                                 max_dec_level = L,
                                 dm = dm_ldbk,
                                 en = TimeFrequency(),
                                 dp = dp,
                                 top_k = k,
                                 n_features = n_features)
ldbkash = LocalDiscriminantBasis(wt = wt,
                                 max_dec_level = L,
                                 dm = dm_ldbkash,
                                 en = ProbabilityDensity(),
                                 dp = dp,
                                 top_k = k,
                                 n_features = n_features)
ldbkemd = LocalDiscriminantBasis(wt = wt,
                                 max_dec_level = L,
                                 dm = dm_ldbkemd,
                                 en = Signatures(),
                                 dp = dp,
                                 top_k = k,
                                 n_features = n_features)
ldbs = Dict("ldbk" => ldbk, "ldbkash" => ldbkash, "ldbkemd" => ldbkemd)
# TODO: Make hyperparameter changes if necessary
LDA = @load LDA pkg=MultivariateStats
CT = @load DecisionTreeClassifier pkg=DecisionTrees
classifiers = Dict("LDA" => LDA(), "CT" => CT())
# TODO: Make changes to model evaluation measures if necessary
measures = [MisclassificationRate()]

## ========== Run experiments ==========
results = repeat_experiment(ldbs, classifiers, measurs)


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

JLLWrappers.@generate_main_file_header("ssht")
JLLWrappers.@generate_main_file("ssht", UUID("06c0f48c-727e-59bc-bc79-048e0b09db0c"))
end  # module ssht_jll


text:
#
# count the number of surreals created by generation n
#    https://oeis.org/A174972
#

function a(n::Integer)
    a = zero(n)
    for i=1:2^n
        a += i * binomial(2^n - 1, i-1)
    end
    return a
end

# 8 is as high as OEIS goes
# approximately doubles the number of digits each time
for n=zero(BigInt) : 9
    println("  n = $n, a(n) = ", a(n), ", approx(a(n)) = ", float(a(n)))
end


text:
using HDF5
using UpROOT
using DataFrames

function GetHitInformation(filename::String)::Union{Nothing, GroupedDataFrame}
    extension = filename[findlast(isequal('.'), filename):end]

    if extension == ".root"
        return GetROOTHitInformation(filename)
    elseif extension == ".hdf5"
        return GetHDF5HitInformation(filename)
    end
    return nothing
end

function GetHDF5HitInformation(filename::String)::GroupedDataFrame
    @info "Reading Geant4 Hits info from HDF5"
    fid = h5open(filename, "r")
    id::Array{Int64, 1} = read(fid, "ntuple/hits/iD/pages")
    x::Array{Float64, 1} = read(fid, "ntuple/hits/x/pages")
    y::Array{Float64, 1} = read(fid, "ntuple/hits/y/pages")
    z::Array{Float64, 1} = read(fid, "ntuple/hits/z/pages")
    eDep::Array{Float64, 1} = read(fid, "ntuple/hits/eDep/pages")

    df = DataFrame(ID = id, X = x, Y = y, Z = z, E = eDep)

    return groupby(df, :ID)
end

function GetROOTHitInformation(filename::String)::GroupedDataFrame
    # Functional, but might be crappy performance. Might look into lazy evaluation?
    @info "Reading Geant4 Hits info from ROOT $filename"
    file = TFile(filename)
    t = file["ntuple/hits"]

    df = DataFrame(ID = t.iD, X = t.x, Y = t.y, Z = t.z, E = t.eDep)

    return groupby(df, :ID)
end

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
################################################################################

@testset "Sort" begin

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

  # load small array into hidden Markov model
  include("data/smallArray.jl")

  # load small array on model data
  for ι ∈ 1:10
    push!(hmm.data, x[ι, :])
  end

  # sort model data
  HiddenMarkovModelReaders.sortHMM!(hmm)

  # testing
  @test hmm.data == [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [6898.439972, 1674.825584, 2142.319328, 923.535397, 562.254679, 542.897036, 544.597461, 447.280304, 338.127405, 249.294769, 331.028492, 335.807143, 421.090866, 236.912484, 241.648853, 423.595156], [8597.829269, 1833.756462, 1628.936111, 1163.014809, 697.425591, 718.742764, 431.683627, 524.160545, 522.474073, 404.779799, 286.859573, 291.428293, 371.560964, 297.910174, 426.403369, 289.336908], [8926.864639, 2024.670303, 1371.395947, 1277.803127, 978.470009, 456.044172, 580.524843, 546.832397, 342.843889, 706.63784, 408.624054, 324.00528, 337.768628, 331.042826, 444.800913, 449.033676], [9221.546188, 1511.961166, 1357.589783, 1006.436895, 532.545324, 605.403816, 468.057093, 388.927434, 439.460726, 281.329487, 469.301092, 309.418703, 340.09914, 381.352773, 227.883408, 402.3395], [9112.851349, 1626.732288, 1602.905491, 1432.479465, 663.49947, 804.251534, 739.028986, 512.910912, 397.43222, 389.500046, 407.156647, 517.501736, 347.87634, 241.200329, 416.171814, 341.232077], [9246.340912, 1466.947083, 1906.052313, 1209.67667, 606.815348, 520.541983, 571.563523, 557.675089, 473.436043, 402.273785, 345.846071, 505.062222, 376.109373, 439.998122, 406.717417, 398.954038], [10628.121246, 1498.386951, 1187.115562, 925.775265, 694.86956, 257.954501, 592.700003, 492.023813, 452.887017, 563.844071, 473.371099, 296.826092, 479.452673, 288.918503, 369.519865, 430.178178], [10362.170502, 2234.273811, 2244.926224, 647.264278, 655.394653, 571.437164, 662.987797, 486.684425, 395.458192, 392.33078, 349.11353, 403.472491, 409.692116, 500.905411, 547.84833, 789.415714], [11321.235947, 1362.566475, 944.705898, 696.605307, 411.700334, 517.204907, 383.636882, 489.580448, 301.713282, 371.661125, 476.757751, 266.741815, 154.978138, 412.402929, 279.845093, 382.712366], [11443.098831, 1500.468246, 1641.78447, 929.626314, 695.049788, 473.568182, 568.370424, 542.348984, 449.226664, 351.930719, 440.918041, 346.611481, 539.18642, 325.506171, 361.940255, 451.27547]]
  @test hmm.model == [[0.0; repeat([-1.0], size(x, 1))]]

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

end

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


text:
using DiscreteValueIteration
using POMDPModels
using POMDPModelTools
using POMDPPolicies
using POMDPs
using Test
import POMDPLinter: @show_requirements

# special grid world for some tests:
#  - reward(m, s, a, sp)
#  - constrained actions
mutable struct SpecialGridWorld <: MDP{GridWorldState, GridWorldAction}
    gw::LegacyGridWorld
end

POMDPs.discount(g::SpecialGridWorld) = discount(g.gw)
POMDPs.transition(g::SpecialGridWorld, s::GridWorldState, a::GridWorldAction) = transition(g.gw, s, a)
POMDPs.reward(g::SpecialGridWorld, s::GridWorldState, a::GridWorldAction, sp::GridWorldState) = reward(g.gw, s, a, sp)
POMDPs.stateindex(g::SpecialGridWorld, s::GridWorldState) = stateindex(g.gw, s)
POMDPs.actionindex(g::SpecialGridWorld, a::GridWorldAction) = actionindex(g.gw, a)
POMDPs.actions(g::SpecialGridWorld, s::GridWorldState) = actions(g.gw, s)
POMDPs.states(g::SpecialGridWorld) = states(g.gw)
POMDPs.actions(g::SpecialGridWorld) = actions(g.gw)

SpecialGridWorld() = SpecialGridWorld(LegacyGridWorld(sx=2, sy=3, rs = [GridWorldState(2,3)], rv = [10.0]))

# Let's extend actions to hard-code & limit to the
# particular feasible actions from each state....
function POMDPs.actions(mdp::SpecialGridWorld, s::GridWorldState)
    # up: 1, down: 2, left: 3, right: 4
    sidx = stateindex(mdp, s)
    if sidx == 1
        acts = [GridWorldAction(:left), GridWorldAction(:right)]
    elseif sidx == 2
        acts = [GridWorldAction(:up), GridWorldAction(:right)]
    elseif sidx == 3
        acts = [GridWorldAction(:left), GridWorldAction(:up)]
    elseif sidx == 4
        acts = [GridWorldAction(:left), GridWorldAction(:right)]
    elseif sidx == 5
        acts = [GridWorldAction(:left), GridWorldAction(:right)]
    elseif sidx == 6
        acts = [GridWorldAction(:up), GridWorldAction(:down), GridWorldAction(:left), GridWorldAction(:right)]
    elseif sidx == 7
        acts = [GridWorldAction(:up), GridWorldAction(:down), GridWorldAction(:left), GridWorldAction(:right)]
    end
    return acts
end

# Test basic value iteration functionality
@testset "all" begin
    @testset "policy" begin
        include("test_value_iteration_policy.jl") # test the policy object
    end
    @testset "basic" begin
        include("test_basic_value_iteration.jl")  # then the creation of a policy
    end
    @testset "basic disallowing actions" begin
        include("test_basic_value_iteration_disallowing_actions.jl") # then a complex form where states determine actions
    end

    @testset "sparse" begin
        include("test_sparse.jl")
    end

    println("Testing Requirements")
    @show_requirements solve(ValueIterationSolver(), LegacyGridWorld())
end