text:
"""
    Arrow.write(io::IO, tbl)
    Arrow.write(file::String, tbl)

Write any Tables.jl-compatible `tbl` out as arrow formatted data.
Providing an `io::IO` argument will cause the data to be written to it
in the "streaming" format, while providing a `file::String` argument
will result in the "file" format being written.

Multiple record batches will be written based on the number of
`Tables.partitions(tbl)` that are provided; by default, this is just
one for a given table, but some table sources support automatic
partitioning. Note you can turn multiple table objects into partitions
by doing `Tables.partitioner([tbl1, tbl2, ...])`, but do remember that
each table must have the exact same `Tables.Schema`.
"""
function write end

function write(file::String, tbl; debug::Bool=false)
    open(file, "w") do io
        write(io, tbl, true, debug)
    end
    return file
end

function write(io::IO, tbl; debug::Bool=false, file::Bool=false)
    return write(io, tbl, file, debug)
end

if isdefined(Tables, :partitions)
    parts = Tables.partitions
else
    parts(x) = (x,)
    parts(x::Tuple) = x
end

function write(io, source, writetofile, debug)
    if writetofile
        Base.write(io, "ARROW1\0\0")
    end
    msgs = Channel{Message}(Inf)
    # build messages
    sch = Ref{Tables.Schema}()
    firstcols = Ref{Any}()
    blocks = (Block[], Block[])
    dictid = Ref(0)
    dictencodings = Dict{Int, Tuple{Int, Type, Any}}()
    # start message writing from channel
@static if VERSION >= v"1.3-DEV"
    tsk = Threads.@spawn for msg in msgs
        Base.write(io, msg, blocks)
    end
else
    tsk = @async for msg in msgs
        Base.write(io, msg, blocks)
    end
end
    @sync for (i, tbl) in enumerate(parts(source))
        if i == 1
            cols = Tables.columns(tbl)
            sch[] = Tables.schema(cols)
            firstcols[] = cols
            for (i, col) in enumerate(Tables.Columns(cols))
                if col isa DictEncode
                    id = dictid[]
                    dictid[] += 1
                    values = unique(col)
                    dictencodings[i] = (id, encodingtype(length(values)), values)
                end
            end
            put!(msgs, makeschemamsg(sch[], cols, dictencodings))
            if !isempty(dictencodings)
                for (colidx, (id, T, values)) in dictencodings
                    dictsch = Tables.Schema((:col,), (eltype(values),))
                    put!(msgs, makedictionarybatchmsg(dictsch, (col=values,), id, false, debug))
                end
            end
            put!(msgs, makerecordbatchmsg(sch[], cols, dictencodings, debug))
        else
@static if VERSION >= v"1.3-DEV"
            Threads.@spawn begin
                try
                    cols = Tables.columns(tbl)
                    if !isempty(dictencodings)
                        for (colidx, (id, T, values)) in dictencodings
                            dictsch = Tables.Schema((:col,), (eltype(values),))
                            col = Tables.getcolumn(cols, colidx)
                            # get new values we haven't seen before for delta update
                            vals = setdiff(col, values)
                            put!(msgs, makedictionarybatchmsg(dictsch, (col=vals,), id, true, debug))
                            # add new values to existing set for future diffs
                            union!(values, vals)
                        end
                    end
                    put!(msgs, makerecordbatchmsg(sch[], cols, dictencodings, debug))
                catch e
                    showerror(stdout, e, catch_backtrace())
                    rethrow(e)
                end
            end
else
            @async begin
                try
                    cols = Tables.columns(tbl)
                    if !isempty(dictencodings)
                        for (colidx, (id, T, values)) in dictencodings
                            dictsch = Tables.Schema((:col,), (eltype(values),))
                            col = Tables.getcolumn(cols, colidx)
                            # get new values we haven't seen before for delta update
                            vals = setdiff(col, values)
                            put!(msgs, makedictionarybatchmsg(dictsch, (col=vals,), id, true, debug))
                            # add new values to existing set for future diffs
                            union!(values, vals)
                        end
                    end
                    put!(msgs, makerecordbatchmsg(sch[], cols, dictencodings, debug))
                catch e
                    showerror(stdout, e, catch_backtrace())
                    rethrow(e)
                end
            end
end
        end
    end
    close(msgs)
    wait(tsk)
    # write empty message
    if !writetofile
        Base.write(io, Message(UInt8[], nothing, nothing, 0, true, false), blocks)
    end
    if writetofile
        b = FlatBuffers.Builder(1024)
        schfoot = makeschema(b, sch[], firstcols[], dictencodings)
        if !isempty(blocks[1])
            N = length(blocks[1])
            Meta.footerStartRecordBatchesVector(b, N)
            for blk in Iterators.reverse(blocks[1])
                Meta.createBlock(b, blk.offset, blk.metaDataLength, blk.bodyLength)
            end
            recordbatches = FlatBuffers.endvector!(b, N)
        else
            recordbatches = FlatBuffers.UOffsetT(0)
        end
        if !isempty(blocks[2])
            N = length(blocks[2])
            Meta.footerStartDictionariesVector(b, N)
            for blk in Iterators.reverse(blocks[2])
                Meta.createBlock(b, blk.offset, blk.metaDataLength, blk.bodyLength)
            end
            dicts = FlatBuffers.endvector!(b, N)
        else
            dicts = FlatBuffers.UOffsetT(0)
        end
        Meta.footerStart(b)
        Meta.footerAddVersion(b, Meta.MetadataVersion.V4)
        Meta.footerAddSchema(b, schfoot)
        Meta.footerAddDictionaries(b, dicts)
        Meta.footerAddRecordBatches(b, recordbatches)
        foot = Meta.footerEnd(b)
        FlatBuffers.finish!(b, foot)
        footer = FlatBuffers.finishedbytes(b)
        Base.write(io, footer)
        Base.write(io, Int32(length(footer)))
        Base.write(io, "ARROW1")
    end
    return io
end

struct Message
    msgflatbuf
    columns
    dictencodings
    bodylen
    isrecordbatch::Bool
    blockmsg::Bool
end

struct Block
    offset::Int64
    metaDataLength::Int32
    bodyLength::Int64
end

struct DictEncoder{T, A, D} <: AbstractVector{T}
    values::A
    pool::D # Dict{eltype(A), I} where I is index type
end

DictEncoder(values::A, unique::B, ::Type{T}) where {A, B, T} =
    DictEncoder{T, A, Dict{eltype(A), T}}(values, Dict{eltype(A), T}(x => (i - 1) for (i, x) in enumerate(unique)))

Base.IndexStyle(::Type{<:DictEncoder}) = Base.IndexLinear()
Base.size(x::DictEncoder) = (length(x.values),)
Base.eltype(x::DictEncoder{T, A}) where {T, A} = T
Base.getindex(x::DictEncoder, i::Int) = x.pool[x.values[i]]

function Base.write(io::IO, msg::Message, blocks)
    metalen = padding(length(msg.msgflatbuf))
    if msg.blockmsg
        push!(blocks[msg.isrecordbatch ? 1 : 2], Block(position(io), metalen + 8, msg.bodylen))
    end
    # now write the final message spec out
    # continuation byte
    n = Base.write(io, 0xFFFFFFFF)
    # metadata length
    n += Base.write(io, Int32(metalen))
    # message flatbuffer
    n += Base.write(io, msg.msgflatbuf)
    n += writezeros(io, paddinglength(n))
    # message body
    if msg.columns !== nothing
        # write out buffers
        for i = 1:length(Tables.columnnames(msg.columns))
            col = Tables.getcolumn(msg.columns, i)
            if msg.dictencodings !== nothing && haskey(msg.dictencodings, i)
                _, T, vals = msg.dictencodings[i]
                col = DictEncoder(col, vals, T)
            end
            writebuffer(io, eltype(col) === Missing ? Missing : Base.nonmissingtype(eltype(col)), col)
        end
    end
    return n
end

function makemessage(b, headerType, header, columns=nothing, dictencodings=nothing, bodylen=0)
    # write the message flatbuffer object
    Meta.messageStart(b)
    Meta.messageAddVersion(b, Meta.MetadataVersion.V4)
    Meta.messageAddHeaderType(b, headerType)
    Meta.messageAddHeader(b, header)
    Meta.messageAddBodyLength(b, Int64(bodylen))
    # Meta.messageAddCustomMetadata(b, meta)
    # Meta.messageStartCustomMetadataVector(b, num_meta_elems)
    msg = Meta.messageEnd(b)
    FlatBuffers.finish!(b, msg)
    return Message(FlatBuffers.finishedbytes(b), columns, dictencodings, bodylen, headerType == Meta.RecordBatch, headerType == Meta.RecordBatch || headerType == Meta.DictionaryBatch)
end

function makeschema(b, sch::Tables.Schema{names, types}, columns, dictencodings) where {names, types}
    # build Field objects
    N = length(names)
    fieldoffsets = [fieldoffset(b, i, names[i], fieldtype(types, i), dictencodings) for i = 1:N]
    Meta.schemaStartFieldsVector(b, N)
    for off in Iterators.reverse(fieldoffsets)
        FlatBuffers.prependoffset!(b, off)
    end
    fields = FlatBuffers.endvector!(b, N)
    # write schema object
    Meta.schemaStart(b)
    Meta.schemaAddEndianness(b, Meta.Endianness.Little)
    Meta.schemaAddFields(b, fields)
    # Meta.schemaAddCustomMetadata(b, meta)
    return Meta.schemaEnd(b)
end

function makeschemamsg(sch::Tables.Schema{names, types}, columns, dictencodings) where {names, types}
    b = FlatBuffers.Builder(1024)
    schema = makeschema(b, sch, columns, dictencodings)
    return makemessage(b, Meta.Schema, schema)
end

function fieldoffset(b, colidx, name, T, dictencodings)
    nameoff = FlatBuffers.createstring!(b, String(name))
    nullable = T >: Missing
    # build dictionary
    if dictencodings !== nothing && haskey(dictencodings, colidx)
        id, encodingtype, _ = dictencodings[colidx]
        _, inttype, _ = arrowtype(b, encodingtype)
        Meta.dictionaryEncodingStart(b)
        Meta.dictionaryEncodingAddId(b, Int64(id))
        Meta.dictionaryEncodingAddIndexType(b, inttype)
        # TODO: support isOrdered?
        Meta.dictionaryEncodingAddIsOrdered(b, false)
        dict = Meta.dictionaryEncodingEnd(b)
    else
        dict = FlatBuffers.UOffsetT(0)
    end
    type, typeoff, children = arrowtype(b, maybemissing(T))
    if children !== nothing
        Meta.fieldStartChildrenVector(b, length(children))
        for off in Iterators.reverse(children)
            FlatBuffers.prependoffset!(b, off)
        end
        children = FlatBuffers.endvector!(b, length(children))
    else
        Meta.fieldStartChildrenVector(b, 0)
        children = FlatBuffers.endvector!(b, 0)
    end
    # build field object
    Meta.fieldStart(b)
    Meta.fieldAddName(b, nameoff)
    Meta.fieldAddNullable(b, nullable)
    Meta.fieldAddTypeType(b, type)
    Meta.fieldAddType(b, typeoff)
    Meta.fieldAddDictionary(b, dict)
    Meta.fieldAddChildren(b, children)
    # Meta.fieldAddCustomMetadata(b, meta)
    return Meta.fieldEnd(b)
end

struct FieldNode
    length::Int64
    null_count::Int64
end

struct Buffer
    offset::Int64
    length::Int64
end

function makerecordbatchmsg(sch::Tables.Schema{names, types}, columns, dictencodings, debug) where {names, types}
    b = FlatBuffers.Builder(1024)
    recordbatch, bodylen = makerecordbatch(b, sch, columns, dictencodings, debug)
    return makemessage(b, Meta.RecordBatch, recordbatch, columns, dictencodings, bodylen)
end

function makerecordbatch(b, sch::Tables.Schema{names, types}, columns, dictencodings, debug) where {names, types}
    nrows = Tables.rowcount(columns)
    debug && println("building record batch message for $nrows rows")

    fieldnodes = FieldNode[]
    fieldbuffers = Buffer[]
    makenodesbuffers!(1, types, columns, dictencodings, fieldnodes, fieldbuffers, 0)
    debug && @show fieldnodes, fieldbuffers

    # write field nodes objects
    FN = length(fieldnodes)
    Meta.recordBatchStartNodesVector(b, FN)
    for fn in Iterators.reverse(fieldnodes)
        Meta.createFieldNode(b, fn.length, fn.null_count)
    end
    nodes = FlatBuffers.endvector!(b, FN)

    # write buffer objects
    bodylen = 0
    BN = length(fieldbuffers)
    Meta.recordBatchStartBuffersVector(b, BN)
    for buf in Iterators.reverse(fieldbuffers)
        Meta.createBuffer(b, buf.offset, buf.length)
        bodylen += padding(buf.length)
    end
    buffers = FlatBuffers.endvector!(b, BN)
    debug && @show bodylen

    # write record batch object
    Meta.recordBatchStart(b)
    Meta.recordBatchAddLength(b, Int64(nrows))
    Meta.recordBatchAddNodes(b, nodes)
    Meta.recordBatchAddBuffers(b, buffers)
    return Meta.recordBatchEnd(b), bodylen
end

function makedictionarybatchmsg(sch::Tables.Schema{names, types}, columns, id, isdelta, debug) where {names, types}
    b = FlatBuffers.Builder(1024)
    recordbatch, bodylen = makerecordbatch(b, sch, columns, nothing, debug)
    Meta.dictionaryBatchStart(b)
    Meta.dictionaryBatchAddId(b, Int64(id))
    Meta.dictionaryBatchAddData(b, recordbatch)
    Meta.dictionaryBatchAddIsDelta(b, isdelta)
    dictionarybatch = Meta.dictionaryBatchEnd(b)
    return makemessage(b, Meta.DictionaryBatch, dictionarybatch, columns, nothing, bodylen)
end

function makenodesbuffers!(colidx, types, columns, dictencodings, fieldnodes, fieldbuffers, bufferoffset)
    colidx > fieldcount(types) && return
    T = fieldtype(types, colidx)
    col = Tables.getcolumn(columns, colidx)
    if dictencodings !== nothing && haskey(dictencodings, colidx)
        _, T, vals = dictencodings[colidx]
        col = DictEncoder(col, vals, T)
        len = _length(col)
        nc = nullcount(col)
        push!(fieldnodes, FieldNode(len, nc))
        # validity bitmap
        blen = nc == 0 ? 0 : bitpackedbytes(len)
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        # adjust buffer offset, make array buffer
        bufferoffset += blen
        blen = sizeof(T) * len
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        bufferoffset += padding(blen)
    else
        bufferoffset = makenodesbuffers!(maybemissing(T), col, fieldnodes, fieldbuffers, bufferoffset)
    end
    # make next column node/buffers
    return makenodesbuffers!(colidx + 1, types, columns, dictencodings, fieldnodes, fieldbuffers, bufferoffset)
end

function makenodesbuffers!(::Type{Missing}, col, fieldnodes, fieldbuffers, bufferoffset)
    len = _length(col)
    push!(fieldnodes, FieldNode(len, len))
    return bufferoffset
end

function writebuffer(io, ::Type{Missing}, col)
    return
end

function makenodesbuffers!(::Type{T}, col, fieldnodes, fieldbuffers, bufferoffset) where {T}
    len = _length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    # adjust buffer offset, make primitive array buffer
    bufferoffset += blen
    blen = len * sizeof(T)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    return bufferoffset + padding(blen)
end

makenodesbuffers!(::Type{Dates.Date}, col, fieldnodes, fieldbuffers, bufferoffset) =
    makenodesbuffers!(Date{Meta.DateUnit.DAY, Int32}, converter(Date{Meta.DateUnit.DAY, Int32}, col), fieldnodes, fieldbuffers, bufferoffset)
makenodesbuffers!(::Type{Dates.Time}, col, fieldnodes, fieldbuffers, bufferoffset) =
    makenodesbuffers!(Time{Meta.TimeUnit.NANOSECOND, Int64}, converter(Time{Meta.TimeUnit.NANOSECOND, Int64}, col), fieldnodes, fieldbuffers, bufferoffset)
makenodesbuffers!(::Type{Dates.DateTime}, col, fieldnodes, fieldbuffers, bufferoffset) =
    makenodesbuffers!(Date{Meta.DateUnit.MILLISECOND, Int64}, converter(Date{Meta.DateUnit.MILLISECOND, Int64}, col), fieldnodes, fieldbuffers, bufferoffset)
makenodesbuffers!(::Type{P}, col, fieldnodes, fieldbuffers, bufferoffset) where {P <: Dates.Period} =
    makenodesbuffers!(Duration{arrowperiodtype(P)}, converter(Duration{arrowperiodtype(P)}, col), fieldnodes, fieldbuffers, bufferoffset)

function writebitmap(io, col)
    nullcount(col) == 0 && return 0
    len = _length(col)
    i = 0
    n = 0
    st = iterate(col)
    for _ = 1:bitpackedbytes(len)
        b = 0x00
        for j = 1:8
            if (i + j) <= len
                x, state = st
                b = setbit(b, !ismissing(x), j)
                st = iterate(col, state)
            end
        end
        n += Base.write(io, b)
        i += 8
    end
    return n
end

function writebuffer(io, ::Type{T}, col) where {T}
    writebitmap(io, col)
    n = writearray(io, T, col)
    writezeros(io, paddinglength(n))
    return
end

writebuffer(io, ::Type{Dates.Date}, col) = writebuffer(io, Date{Meta.DateUnit.DAY, Int32}, converter(Date{Meta.DateUnit.DAY, Int32}, col))
writebuffer(io, ::Type{Dates.Time}, col) = writebuffer(io, Time{Meta.TimeUnit.NANOSECOND, Int64}, converter(Time{Meta.TimeUnit.NANOSECOND, Int64}, col))
writebuffer(io, ::Type{Dates.DateTime}, col) = writebuffer(io, Date{Meta.DateUnit.MILLISECOND, Int64}, converter(Date{Meta.DateUnit.MILLISECOND, Int64}, col))
writebuffer(io, ::Type{P}, col) where {P <: Dates.Period} = writebuffer(io, Duration{arrowperiodtype(P)}, converter(Duration{arrowperiodtype(P)}, col))

function makenodesbuffers!(::Type{T}, col, fieldnodes, fieldbuffers, bufferoffset) where {T <: Union{AbstractString, AbstractVector}}
    len = _length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    # adjust buffer offset, make array buffer
    bufferoffset += blen
    # TODO: support Large lists
    blen = sizeof(Int32) * (len + 1)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    bufferoffset += padding(blen)
    if T <: AbstractString || T <: AbstractVector{UInt8}
        blen = 0
        for x in col
            blen += sizeof(x)
        end
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        bufferoffset += padding(blen)
    else
        bufferoffset = makenodesbuffers!(maybemissing(eltype(T)), flatten(skipmissing(col)), fieldnodes, fieldbuffers, bufferoffset)
    end
    return bufferoffset
end

function writebuffer(io, ::Type{T}, col) where {T <: Union{AbstractString, AbstractVector}}
    writebitmap(io, col)
    # write offsets
    off::Int32 = 0
    len = T <: AbstractString ? sizeof : length
    n = 0
    for x in col
        n += Base.write(io, off)
        if x !== missing
            off += len(x)
        end
    end
    n += Base.write(io, off)
    writezeros(io, paddinglength(n))
    # write values array
    if T <: AbstractString || T <: AbstractVector{UInt8}
        n = 0
        for x in col
            if x !== missing
                n += Base.write(io, x)
            end
        end
        writezeros(io, paddinglength(n))
    else
        writebuffer(io, maybemissing(eltype(T)), flatten(skipmissing(col)))
    end
    return
end

function makenodesbuffers!(::Type{NTuple{N, T}}, col, fieldnodes, fieldbuffers, bufferoffset) where {N, T}
    len = _length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    bufferoffset += blen
    if T === UInt8
        blen = N * len
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        bufferoffset += padding(blen)
    else
        bufferoffset = makenodesbuffers!(maybemissing(T), flatten(coalesce(x, default(NTuple{N, T})) for x in col), fieldnodes, fieldbuffers, bufferoffset)
    end
    return bufferoffset
end

function writebuffer(io, ::Type{NTuple{N, T}}, col) where {N, T}
    writebitmap(io, col)
    # write values array
    if T === UInt8
        n = writearray(io, NTuple{N, T}, col)
        writezeros(io, paddinglength(n))
    else
        writebuffer(io, maybemissing(T), flatten(coalesce(x, default(NTuple{N, T})) for x in col))
    end
    return
end

function makenodesbuffers!(::Type{Pair{K, V}}, col, fieldnodes, fieldbuffers, bufferoffset) where {K, V}
    # Struct child node
    bufferoffset = makenodesbuffers!(Vector{KeyValue{K, V}}, ( KeyValue(k, v) for (k, v) in pairs(col) ), fieldnodes, fieldbuffers, bufferoffset)
    return bufferoffset
end

function writebuffer(io, ::Type{Pair{K, V}}, col) where {K, V}
    # write values array
    writebuffer(io, Vector{KeyValue{K, V}}, ( KeyValue(k, v) for (k, v) in pairs(col) ))
    return
end

function makenodesbuffers!(::Type{KeyValue{K, V}}, col, fieldnodes, fieldbuffers, bufferoffset) where {K, V}
    len = _length(col)
    push!(fieldnodes, FieldNode(len, 0))
    # validity bitmap
    push!(fieldbuffers, Buffer(bufferoffset, 0))
    # keys
    bufferoffset = makenodesbuffers!(maybemissing(K), (x.key for x in col), fieldnodes, fieldbuffers, bufferoffset)
    # values
    bufferoffset = makenodesbuffers!(maybemissing(V), (x.value for x in col), fieldnodes, fieldbuffers, bufferoffset)
    return bufferoffset
end

function writebuffer(io, ::Type{KeyValue{K, V}}, col) where {K, V}
    writebitmap(io, col)
    # write keys
    writebuffer(io, maybemissing(K),(x.key for x in col))
    # write values
    writebuffer(io, maybemissing(V),(x.value for x in col))
    return
end

function makenodesbuffers!(::Type{NamedTuple{names, types}}, col, fieldnodes, fieldbuffers, bufferoffset) where {names, types}
    len = _length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    bufferoffset += blen
    for i = 1:length(names)
        bufferoffset = makenodesbuffers!(maybemissing(fieldtype(types, i)), (getfield(x, names[i]) for x in col), fieldnodes, fieldbuffers, bufferoffset)
    end
    return bufferoffset
end

function writebuffer(io, ::Type{NamedTuple{names, types}}, col) where {names, types}
    writebitmap(io, col)
    # write values arrays
    for i = 1:length(names)
        writebuffer(io, maybemissing(fieldtype(types, i)), (getfield(x, names[i]) for x in col))
    end
    return
end

function makenodesbuffers!(::Type{UnionT{T, typeIds, U}}, col, fieldnodes, fieldbuffers, bufferoffset) where {T, typeIds, U}
    len = _length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    # typeIds buffer
    push!(fieldbuffers, Buffer(bufferoffset, len))
    bufferoffset += padding(len)
    if T == Meta.UnionMode.Dense
        # offsets buffer
        blen = sizeof(Int32) * len
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        bufferoffset += padding(blen)
        # value arrays
        for i = 1:fieldcount(U)
            S = fieldtype(U, i)
            bufferoffset = makenodesbuffers!(maybemissing(S), filtered(i == 1 ? Union{S, Missing} : maybemissing(S), col), fieldnodes, fieldbuffers, bufferoffset)
        end
    else
        # value arrays
        for i = 1:fieldcount(U)
            S = fieldtype(U, i)
            bufferoffset = makenodesbuffers!(maybemissing(S), replaced(S, col), fieldnodes, fieldbuffers, bufferoffset)
        end
    end
    return bufferoffset
end

isatypeid(x::T, ::Type{types}) where {T, types} = isatypeid(x, fieldtype(types, 1), types, 1)
isatypeid(x::T, ::Type{S}, ::Type{types}, i) where {T, S, types} = x isa S ? i : isatypeid(x, fieldtype(types, i + 1), types, i + 1)

function writebuffer(io, ::Type{UnionT{T, typeIds, U}}, col) where {T, typeIds, U}
    # typeIds buffer
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    n = 0
    for x in col
        typeid = x === missing ? Int8(0) : Int8(typeids[isatypeid(x, U)])
        n += Base.write(io, typeid)
    end
    writezeros(io, paddinglength(n))
    if T == Meta.UnionMode.Dense
        # offset buffer
        n = 0
        offs = zeros(Int32, fieldcount(U))
        for x in col
            idx = x === missing ? 1 : isatypeid(x, U)
            n += Base.write(io, offs[idx])
            offs[idx] += 1
        end
        writezeros(io, paddinglength(n))
        # write values arrays
        for i = 1:fieldcount(U)
            S = fieldtype(U, i)
            writebuffer(io, maybemissing(S), filtered(i == 1 ? Union{S, Missing} : maybemissing(S), col))
        end
    else
        # value arrays
        for i = 1:fieldcount(U)
            S = fieldtype(U, i)
            writebuffer(io, maybemissing(S), replaced(S, col))
        end
    end
    return
end


text:
using InvariantMeasures
using ValidatedNumerics
using LinearAlgebra: I, opnorm
using SparseArrays: sparse
using Test

@testset "Norms of powers" begin

import Random
Random.seed!(1) #ensures tests are deterministic

n = 9
m = 10
M = 0.2*randn(n, n)
R = 1e-14*rand(n, n)

LL = sparse(interval_from_midpoint_radius.(M, R))

Q = IntegralPreservingDiscretizedOperator(LL)

U = [ones(1,n-1); -Matrix(I, n-1,n-1)]

@test norms_of_powers(Linf, m, Q, zeros(Interval{Float64}, 1)) ≈ [opnorm(M^k*U,Inf) for k = 1:m]

@test norms_of_powers(L1, m, Q, zeros(Interval{Float64}, 1)) ≈ [opnorm(M^k*U,1) for k = 1:m]

e = ones(n)
f = adjoint(e) / n
intf = map(Interval, f)
Q = NonIntegralPreservingDiscretizedOperator(LL, e, f)

@test norms_of_powers(Linf, m, Q, intf) ≈ [opnorm((M+e*(f-f*M))^k*U,Inf) for k = 1:m]
@test norms_of_powers(L1, m, Q, intf) ≈ [opnorm((M+e*(f-f*M))^k*U,1) for k = 1:m]

@test refine_norms_of_powers([0.5, 1, 2, 0.001]) == [0.5, 0.25, 0.125, 0.001]
@test refine_norms_of_powers([0.5, 1, 2, 1e-3], 8) == [0.5, 0.25, 0.125, 1e-3, 0.5e-3, 0.25e-3, 0.125e-3, 1.0000000000000002e-6] # also tests correct rounding
@test refine_norms_of_powers([2,0.2,0.1],4) == [2, 0.2, 0.1, 0.04000000000000001]

end

ext: jl


text:
const DS = DataStructure101

export TreeNode, BinaryTree
export value, parent, left_child, right_child, level
export root, tree_repr

mutable struct TreeNode{T} <: DS.AbstractNode
    value::T
    parent::DS.AbstractNode
    left::DS.AbstractNode
    right::DS.AbstractNode

    TreeNode{T}(value::T) where T =new(
        value,
        DS.NullNode(),
        DS.NullNode(),
        DS.NullNode()
    )
end

value(tn::TreeNode) = tn.value

parent(tn::TreeNode) = tn.parent

left_child(tn::TreeNode) = tn.left

right_child(tn::TreeNode) = tn.right

parent(tn::DS.NullNode) = DS.NullNode()

left_child(tn::DS.NullNode) = DS.NullNode()

right_child(tn::DS.NullNode) = DS.NullNode()

function Base.show(io::IO, node::DS.TreeNode{T}) where T
    print(io, "TreeNode{$T}($(value(node)))")
end

function level(node::DS.AbstractNode)
    !(node isa DS.TreeNode) && (return -1)

    level = 0
    while !(parent(node) isa DS.NullNode)
        node = parent(node)
        level += 1
    end

    return level
end

mutable struct BinaryTree
    root::DS.AbstractNode
    length::Int
end

BinaryTree(value::Int64) = BinaryTree(TreeNode{Int64}(value), 1)

BinaryTree() = BinaryTree(DS.NullNode(), 0)

root(bt::BinaryTree) = bt.root

Base.length(bt::BinaryTree) = bt.length

Base.eltype(bt::BinaryTree) = Int64

tree_repr(node::DS.NullNode; tree_str="", level=0) = "$(tree_str)\n"

function tree_repr(node::TreeNode; tree_str="", level=0)
    tree_str = tree_repr(node.right, tree_str=tree_str, level=level+1)
    tree_str = "$(tree_str)\n$("\t"^level)TreeNode($(value(node)))"
    tree_str = tree_repr(node.left, tree_str=tree_str, level=level+1)

    return tree_str
end

function Base.show(io::IO, bt::BinaryTree)
    print(io, "BinaryTree($(tree_repr(root(bt))))")
end

function Base.getindex(bt::BinaryTree, i::Int)
    # calculate path
    path = []
    while i != 1
        push!(path, rem(i, 2))
        i = div(i, 2)
    end

    path = BitArray(reverse(path))

    # search node
    node = root(bt)
    for is_right in path
        if is_right
            node = right_child(node)
        else
            node = left_child(node)
        end
    end

    return node
end

function Base.setindex!(bt::BinaryTree, v::Int64, i::Int64)
    node = TreeNode{Int64}(v)
    if (root(bt) isa DS.NullNode && i == 1)
        bt.root = node
        bt.length += 1
        return
    end

    parent = Base.getindex(bt, div(i, 2))
    is_right = Bool(rem(i, 2))
    !(parent isa TreeNode) && (throw(BoundsError("parent of $(i)th node not found")))
    node.parent = parent
    is_right ? (parent.right = node) : (parent.left = node)

    bt.length += 1
end


text:
using Documenter, TimeseriesSurrogates
ENV["GKSwstype"] = "100"
push!(LOAD_PATH,"../src/")

PAGES = [
    "Overview" => "index.md",
    "What is a surrogate?" => "man/whatisasurrogate.md",
    "Surrogate algorithms" => [
        "Random shuffle (RS)" => "constrained/randomshuffle.md",
        "Fourier transform (FT)" => "constrained/fourier_surrogates.md",
        "Amplitude adjusted Fourier (AAFT)" => "constrained/amplitude_adjusted.md"
    ],
    "Visualising surrogate realizations and creating gifs" => "plots/visualize_surrogate.md",
    "Example systems" => "man/exampleprocesses.md"
]

makedocs(
    modules = [TimeseriesSurrogates],
    format = :html,
    sitename = "TimeseriesSurrogates.jl",
    authors = "Kristian Agasøster Haaga",
    pages = PAGES
    # Use clean URLs, unless built as a "local" build
    #html_prettyurls = !("local" in ARGS),
    #html_canonical = "https://kahaaga.github.io/TimeseriesSurrogates.jl/latest/"
)

deploydocs(
    repo   = "github.com/kahaaga/TimeseriesSurrogates.jl.git",
    julia  = "0.6",
    target = "build",
    deps = nothing,
    make = nothing,
    osname = "linux"
)


text:
struct DownsampledData{T<:TextData, F<:AbstractFloat} <: TextData
    downdata::T
    threshold::F
end

function downsample(data::TextData, threshold::AbstractFloat)
    DownsampledData(data, threshold)
end

function Base.iterate(data::DownsampledData)
    error("you didn't implement this yet, dummy")
end

struct MiniBatchData{T<:Iterators.PartitionIterator{<:TextData}}
    minidata::T
end

function minibatch(data::TextData, n::Integer)
    minidata = (Iterators.partition(data, n))
    MiniBatchData(minidata)
end

Base.length(data::MiniBatchData) = length(data.minidata)

function Base.iterate(data::MiniBatchData, state = 1)
    temp = iterate(data.minidata, state)
    if temp == nothing
        return nothing
    else
        (temp, newstate) = temp
        temp = Iterators.Stateful(temp)
        X, Y = iterate(temp)[1]
        for (i, j) in temp
            X = hcat(X, i)
            Y = vcat(Y, j)
        end
    end
    return ((X, Y), newstate)
end


text:
#####################################################
# Try Bouncy particle sampler with exact bounds     #
#####################################################

using ZigZagBoomerang, SparseArrays, LinearAlgebra, Random
#using CairoMakie
include("../../src/faberschauder.jl")
include("../../src/asvar.jl")
include("../../src/pdmp.jl")

const ZZB = ZigZagBoomerang
using CSV
using DataFrames
using Statistics
# Drift
b(x) = α * sin(x)
# First derivative
b′(x) = α * cos(x)
# Second derivative
b″(x) = -α * sin(x)

struct MyBound
    c::Float64
end
# overloading Poisson times
# WARNING: this function works only for diffusionn bridges and not processes. to generalize it
# change [T^(1.5)/2^((L - lvl(i, L)) * 1.5 + 2) * (α^2 + α) for i in eachindex(x)]
function ZZB.ab(x, θ, c::MyBound, F::BouncyParticle)
    a = c.c + dot([T^(1.5)/2^((L - lvl(i, L)) * 1.5 + 2) * (α^2 + α) for i in eachindex(x)], abs.(θ))
    b1 = dot(x, θ)
    b2 = dot(θ, θ)
    return a, b1, b2
end

function ZZB.adapt!(b::MyBound, x)
    b = MyBound(b.c * x)
end

#Do not reflect the first and the last coefficient
function ZZB.reflect!(∇ϕx, x, θ, F::BouncyParticle)
    θ[2:end-1] .-= (2*dot(∇ϕx[2:end-1], θ[2:end-1])/ZZB.normsq(∇ϕx[2:end-1]))*∇ϕx[2:end-1]
    θ
end
#Do not refresh the first and the last coefficient
function ZZB.refresh!(θ, F::BouncyParticle)
    for i in eachindex(θ)
        θ[i] = randn()
    end
    θ[1] = θ[end] = 0.0
    θ
end
ZZB.sλ(∇ϕx, θ, F) = ZZB.λ(∇ϕx, θ, F)
ZZB.sλ̄((a, b, c)::NTuple{3}, Δt) = a + ZZB.pos(b + c * Δt)

function ∇ϕ!(y, ξ, (L, T); K = 1)
    s = T * (rand())
    x = dotψ(ξ, s, L, T)
    bb = 2b(x) * b′(x) + b″(x)
    for i in eachindex(ξ)
        if i == (2 << L) + 1    # final point
            y[i] = 0.0
        elseif i == 1   # initial point
            y[i] = 0.0
        else
            l = lvl(i, L)
            k = (i - 1) ÷ (2 << l)
            δ = T / (1 << (L - l))
            if δ*k <= s <=   δ*(k + 1)
                y[i] = 0.5 * δ * Λ(s, L - l, T) * bb + ξ[i]
            else
                y[i] =  ξ[i]
            end
        end
    end
    y
end

function run_bps(df, T′)
    c_bps = MyBound(0.0)
    n = (2 << L) + 1
    ξ0, θ0 = randn(n), randn(n)
    θ0[end] = θ0[1] = 0.0 # fix final point
    u, v = 0.0, 0.0  # initial and fianl point
    ξ0[1] = u / sqrt(T)
    ξ0[end] = v / sqrt(T)
    λref_bps = 1.0
    Γ = sparse(1.0I, n, n)
    B = BouncyParticle(Γ, ξ0*0, λref_bps)
    bps_time = @elapsed((bps_trace, uT, acc) = pdmp(∇ϕ!, 0.0, ξ0, θ0, T′, c_bps, B, (L, T); adapt=false, factor=0.0))
    t_bps = getindex.(bps_trace.events, 1)
    x_bps = getindex.(bps_trace.events, 2)
    v_bps = getindex.(bps_trace.events, 3)
    burn1 = 100.0
    i1 = findfirst(x -> x > burn1, t_bps)
    xT_2 = sqrt(T)/2*getindex.(x_bps[i1:end], Int((length(ξ0)+1)/2))
    ave_xT_2 = cumsum(((xT_2[1:end-1] .+ xT_2[2:end])/2).*diff(t_bps[i1:end]))./(t_bps[i1+1:end] .- burn1)
    for batches in [50]
        ess = ess_pdmp_components(t_bps[i1:end] .- burn1, x_bps[i1:end], v_bps[i1:end], n_batches = batches)
        push!(df, Dict(:sampler => "BPS", :alpha => α, :T => T′, :nbatches => batches,
            :stat => "ess_mean", :y => sum(ess[2:end-1])/(length(ess)-2), :runtime => bps_time))
        push!(df, Dict(:sampler => "BPS", :alpha => α, :T => T′, :nbatches => batches,
            :stat => "ess_median", :y => median(ess[2:end-1]), :runtime => bps_time))
        push!(df, Dict(:sampler => "BPS", :alpha => α, :T => T′, :nbatches => batches,
            :stat => "ess_min", :y => minimum(ess[2:end-1]), :runtime => bps_time))
        push!(df, Dict(:sampler => "BPS", :alpha => α, :T => T′, :nbatches => batches,
            :stat => "ess_x_T2", :y => ess[Int((length(ξ0)+1)/2)], :runtime => bps_time))
    end
    return df, (ave_xT_2, t_bps)
end

function data_collection_bps(df)
    Random.seed!(0)
    ave = []
    for α′ in [0.1, 0.3, 0.5]
        global α = α′
        global L = 6
        global T = 50.0
        T′ = 25000
        df, x1 = run_bps(df, T′)
        push!(ave, x1)
    end
    return df, ave
end
df, bps = data_collection_bps(df)
CSV.write("./scripts/zz_diff_bridges/compare/benchamrk_bps_final.csv", df)
@save "convergence_bps.jld2" bps
error("")
# dt = 10.0
# xx = ZZB.trajectory(discretize(out1, dt))
#
# using Makie
# using CairoMakie
# S = T*(0:n)/(n+1)
# p1 = lines(S, [dotψ(xx.x[end], s, L, T) for s in S], linewidth=0.3)
# for ξ in xx.x[1:5:end]
#     lines!(p1, S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.1, alpha = 0.1)
# end
# display(p1)

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
using OptEvapCool
using StatsPlots
using Printf: @sprintf
using LaTeXStrings
using Plots.PlotMeasures

function otago_ghost_trap(duration = 0.24, input_dir = "")
    # Physical parameters
    Np = 1.6e6 # Initial atom count
    T₀ = 2e-6 # Initial temperature
    species = Rb87

    # Numerical parameters
    Nt = ceil(Int64, 1e5) # Test particles
    F = Np / Nt
    Nc = 3 # Target number of test particles per cell

    # Beam parameters
    P_h = 0.52 # Horizontal
    P_v = 0.09 # Vertical
    P_g = 0.29 # Ghost

    w_h = 60e-6 # Beam waist (m)
    w_v = 40e-6
    w_g = 40e-6

    dir_h = [0, 0, 1]
    dir_v = [0, -1, 0]
    dir_g = [0, -1, 0]

    λ = 1064e-9

    # Ghost beam focus
    foc_g = (t) -> [85e-6, 0, 0] + (t/0.24) * [-30e-6, 0, 0]

    beam_h = GaussianBeam([0,0,0], dir_h, P_h, w_h, λ)
    beam_v = GaussianBeam([0,0,0], dir_v, P_v, w_v, λ)
    beam_g = GaussianBeam(foc_g, dir_g, P_g, w_g, λ)

    # Trapping frequencies
    m = species.m
    κ = kappa(species)

    ωr²_v = 8κ*P_v / (m * π * w_v^4) # Radial trapping frequency
    ωr²_h = 8κ*P_h / (m * π * w_h^4)

    ωx = sqrt(ωr²_v + ωr²_h)
    ωy = sqrt(ωr²_h)
    ωz = sqrt(ωr²_v)

    # Cloud initialisation
    positions = harmonic_boltzmann_positions(Nt, m, T₀, ωx, ωy, ωz)
    velocities = boltzmann_velocities(Nt, m, T₀)

    # Translate to equilibrium position
    y₀ = - g / ωy^2
    centre = [0, y₀, 0]
    positions .+= centre

    # Three-body loss & background loss
    K = 4.3e-29 * 1e-12
    τ_bg = 180

    # Maximum timestep
    max_dt = 0.05 * 2π / max(ωx, ωy, ωz)

    # Gravity
    acc_grav! = acceleration(gravity)
    pot_grav! = potential(gravity)

    # GAUSSIAN BEAM SIMULATION
    acc_h! = acceleration(beam_h)
    acc_v! = acceleration(beam_v)
    acc_g! = acceleration(beam_g)

    pot_h! = potential(beam_h)
    pot_v! = potential(beam_v)
    pot_g! = potential(beam_g)

    function accel(p, s, t, o)
        fill!(o, 0.0)
        acc_h!(p, s, t, o)
        acc_v!(p, s, t, o)
        acc_g!(p, s, t, o)
        acc_grav!(p, s, t, o)
        return nothing
    end

    function poten(p, s, t, o)
        fill!(o, 0.0)
        pot_h!(p, s, t, o)
        pot_v!(p, s, t, o)
        pot_g!(p, s, t, o)
        pot_grav!(p, s, t, o)
        return nothing
    end

    sensor = GlobalSensor()
    measure = measurer(sensor, 0.00001, ωx, ωy, ωz, centre)

    T(t) = (length(sensor.ke) > 0) ? 2 * last(sensor.ke) / (3 * kB) : T₀
    evap = ellipsoid_evap(ωx, ωy, ωz, T, 1e-6, centre)

    conditions = SimulationConditions(species, F, positions, velocities,
        accel, poten, evap = evap, τ_bg = τ_bg, K = K)

    if input_dir == ""
        final_cloud = evolve(conditions, duration;
            Nc = Nc, max_dt = max_dt, measure = measure)
    else
        sensor = loadsensor("$(input_dir)/sensor-data.csv")
    end

    # Save sensor data (before plotting, in case plotting fails)
    ft = filetime()
    dir = "./results/$ft-otago-ghost-trap"
    mkpath(dir)
    savecsv(sensor, "$dir/sensor-data.csv")

    # PLOTTING
    default(fontfamily="Computer Modern",
        linewidth=2, framestyle=:box, label=nothing, grid=false)
    scalefontsizes()
    scalefontsizes(1.5)

    window_time = min(duration/2, 1e-2)
    time = sensor.time
    window_size = OptEvapCool.window_time_size(time, window_time)
    rolling_time = rollmean(time, window_size)

    # Colour palette
    col1 = RGB(0.0588, 0, 0.5882)
    #col2 = RGB(0.9000, 0, 0.3490)
    #col3 = RGB(1, 0.6510, 0)
    #col4 = RGB(0.0824, 0.8392, 0)

    # Temperature
    rolling_temp = 2 / (3 * kB) * rollmean(sensor.ke, window_size)

    temp_order = -6
    temp_yformatter(y) = @sprintf("%.2f",y/(10.0^temp_order))

    temperature_plt = plot(rolling_time, rolling_temp,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Temperature\ \ }({}\times10^{%$temp_order}\mathrm{K})",
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        minorticks = 5,
        yformatter = temp_yformatter,
        dpi = 300)

    # Number
    instant_Np = sensor.Nt .* sensor.F

    num_order = 6
    num_yformatter(y) = @sprintf("%.2f",y/(10.0^num_order))

    number_plt = plot(time, instant_Np,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Number\ \ }({}\times10^{%$num_order})",
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        minorticks = 5,
        yformatter = num_yformatter,
        dpi = 300)

    # Collrate
    rolling_timesteps = (
        time[window_size:end] - time[1:end - window_size + 1]
    )
    rolling_Nt = rollmean(sensor.Nt, window_size)
    rolling_Γ = (2 * rolling(sum, sensor.coll, window_size) ./
        (rolling_Nt .* rolling_timesteps)
    )

    collrate_plt = plot(rolling_time, rolling_Γ,
        xlabel = "Time (s)",
        ylabel = "Collision rate (Hz)",
        label = false,
        linecolor = col1,
        ls = :solid,
        dpi = 300)

    # Phase space density
    rolling_n0 = rollmean(sensor.n0 .* sensor.F, window_size)
    rolling_psd = rolling_n0 .* (
        (2*π*h̄^2 ./ (m * kB .* rolling_temp)) .^1.5
    )

    psd_plt = plot(rolling_time, rolling_psd,
        xlabel = "Time (s)",
        ylabel = "Phase space density",
        label = false,
        linecolor = col1,
        ylims = (0, Inf),
        minorticks = 5,
        ls = :solid,
        dpi = 300)

    # Number vs PSD
    rolling_Np = rollmean(instant_Np, window_size)

    numpsd_plt = plot(rolling_psd, rolling_Np,
        xlabel = "Phase space density",
        ylabel = L"\textrm{Number\ \ }({}\times10^{%$num_order})",
        yformatter = num_yformatter,
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, 1.8e6),
        #xlims = (0, 2.8),
        right_margin = 5mm,
        minorticks = 5,
        dpi = 300)
    scatter!([2.6], [1.4e5], color = "black", label = false,
        markersize = 8, markershape = :utriangle)
    #vline!([2.6], line = (:black, 5))

    # Save plots
    savefig(temperature_plt, "$dir/temp.png")
    savefig(collrate_plt, "$dir/collrate.png")
    savefig(number_plt, "$dir/number.png")
    savefig(psd_plt, "$dir/psd.png")
    savefig(numpsd_plt, "$dir/numpsd.png")

    return nothing
end


text:
using Makie

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    scene = Scene()
    scatter(scene, rand(50), rand(50), markersize = 0.01)
    a = axis(scene, range(0, stop = 1, length = 4), range(0, stop = 1, length = 4), textsize = 0.1, axisnames = ("", "", ""))
    tf = to_value(a, :tickfont2d)
    a[:tickfont2d] = map(x-> (0.07, x[2:end]...), tf)
    center!(scene)
    Makie.block(scene)
    return 0
end


text:
using Polyhedra
using LazySets
using MeshCat
using HierarchicalGeometry

vis = Visualizer()
render(vis)

ball = LazySets.Ball2(zeros(3),1.0)
ball2 = LazySets.Ball2([1.0,0.0,0.0],1.0)
lazy_set = UnionSet(ball,ball2)

bbox = overapproximate(ball)
hpoly = convert(LazySets.HPolytope,bbox)
vpoly = tovrep(hpoly)
poly = Polyhedra.polyhedron(hpoly)

# model = RegularPolyhedronOverapprox(3,8)
model = HierarchicalGeometry.equatorial_overapprox_model()
# model = HierarchicalGeometry.equatorial_overapprox_model([0.0],0.0:π/2:2π)
hpoly = overapproximate(lazy_set,model)
poly = Polyhedra.polyhedron(hpoly)

setobject!(vis, Polyhedra.Mesh(poly))
# open(vis)


text:
"""
    main_ptdf
A main file to generate a PTDF file

Parameters:
    input_path : path to input data directory describing a grid
                (not the pscopf_ files but branches.txt and buses.txt)
"""

using Dates
using Printf

root_path = dirname(@__DIR__)
push!(LOAD_PATH, root_path);
cd(root_path)
include(joinpath(root_path, "src", "PTDF.jl"));

#########################
# INPUT & PARAMS
#########################

input_path = ( length(ARGS) > 0 ? ARGS[1] :
                    joinpath(@__DIR__, "..", "data", "ptdf", "3buses_3branches") )
ref_bus_num = 1
distributed = true

function compute(network_p, ref_bus_num_p, distributed_p, i_cut_branch_p=nothing)
    if isnothing(i_cut_branch_p)
        network_l = network_p
        name = ""
    else
        cut_branch_l, network_l = PTDF.reduced_network(network_p, i_cut_branch_p)
        name = "_n1_"*cut_branch_l.name
    end

    ptdf_l = PTDF.compute_ptdf(network_l, ref_bus_num_p)
    if distributed_p
        ptdf_l = PTDF.distribute_slack(ptdf_l);
        # coeffs = Dict([ "poste_1_0" => .2,
        #                 "poste_2_0" => .8])
        # ptdf_l = PTDF.distribute_slack(ptdf_l, coeffs, network_l);
    end
    filename = @sprintf("pscopf_ptdf%s.txt", name)
    output_path = joinpath(input_path, filename)
    #use original network to print full ptdf containing the cut branch (with 0 coeffs)
    PTDF.write_PTDF(output_path, network_p, ptdf_l, distributed_p, ref_bus_num_p, i_cut_branch_p)
end

#########################
# EXECUTION
#########################
network = PTDF.read_network(input_path)

compute(network, ref_bus_num, distributed)
for (i_cut_branch,_) in network.branch_to_i
    compute(network, ref_bus_num, distributed, i_cut_branch)
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
function readgeno!(geno::Matrix{Float64}, curSNPSize::Int,
                   nPer::Int, SNPSize::Int, idxRead::Int,
                   bin2geno::Matrix{Float64},
                   rawdata::Matrix{Int8}, offset::Int,
                   flagAnnotate::Bool;
                   offsetSize::Vector{Int64} = Int64[])

  if !flagAnnotate
    offset = (idxRead - 1) * SNPSize;
  else
    offset = offsetSize[idxRead];
  end
  nrow = iceil(nPer / 4);

  for j = 1 : curSNPSize
    cumPer = 0;

    for i = 1 : nrow - 1
      curvalue = rawdata[i, offset + j];
      if curvalue >= 0
        curvalue += 1;
      else
        curvalue += 257;
      end
      for k = 1 : 4
        cumPer = cumPer + 1;
        geno[cumPer, j] = bin2geno[k, curvalue];
      end
    end

    curvalue = rawdata[nrow, offset + j];
    if curvalue >= 0
      curvalue += 1;
    else
      curvalue += 257;
    end
    for k = 1 : 4
      if cumPer < nPer
        cumPer = cumPer + 1;
        geno[cumPer, j] = bin2geno[k, curvalue];
      else
        break;
      end
    end

  end

end


text:
using Colors
using LightGraphs
using VisJS

c = distinguishable_colors(10)
label = ["Node $i" for i in 1:10]
g = WheelGraph(10)
VisJS.drawGraph(g,label=label,color=c)


text:
export testcase1a, testcase1a_point, testcase1b, testcase1c
export testcase1abis, testcase1ater
export testcase2
export testcase4a, testcase4b, testcase4c

# Test case 1.a
function testcase1a_point()
    x = PointE([Matrix{Float64}(I, 3, 3)], Float64[])
end

function testcase1a(; symmetric=false)
    if symmetric
        cmat = [Symmetric(sparse([1, 2, 3], [1, 2, 3], [1., 1., 1.]))]
        a1mats = [Symmetric(sparse([1], [1], [1.], 3, 3), :L)]
        a2mats = [Symmetric(sparse([2, 3], [2, 1], [1., 1.], 3, 3), :L)]
        a3mats = [Symmetric(sparse([2, 3], [1, 3], [1., 1.], 3, 3), :L)]
    else
        cmat = [sparse([1, 2, 3], [1, 2, 3], [1., 1., 1.])]
        a1mats = [sparse([1], [1], [1.], 3, 3)]
        a2mats = [sparse([1, 2, 3], [3, 2, 1], [1., 1., 1.], 3, 3)]
        a3mats = [sparse([1, 2, 3], [2, 1, 3], [1., 1., 1.], 3, 3)]
    end

    c = PointE(cmat, Vector{Float64}())
    a1 = PointE(a1mats, Vector{Float64}())
    a2 = PointE(a2mats, Vector{Float64}())
    a3 = PointE(a3mats, Vector{Float64}())

    A = [a1, a2, a3]

    b = Vector{Float64}([1, 1, 1])

    problem = SDCOContext(c, A, b)

    x = PointE([Matrix{Float64}(I, 3, 3)], Float64[])
    s = PointE([Matrix{Float64}(I, 3, 3)], Float64[])
    y = zeros(Float64, length(b))

    return problem, PointPrimalDual(x, y, s)
end

function testcase1abis()
    cmat = [sparse([1, 2, 3], [1, 2, 3], [1., 1., 1.])]
    c = PointE(cmat, Vector{Float64}())

    a1mats = [sparse([1], [1], [1.], 3, 3)]
    a1 = PointE(a1mats, Vector{Float64}())

    a2mats = [sparse([2, 3], [2, 1], [1., 1.], 3, 3)]
    a2 = PointE(a2mats, Vector{Float64}())

    a3mats = [sparse([2, 3], [1, 3], [1., 1.], 3, 3)]
    a3 = PointE(a3mats, Vector{Float64}())

    A = [a1, a2, a3]

    b = Vector{Float64}([1, 1, 1])

    problem = SDCOContext(c, A, b)

    x = PointE([Matrix{Float64}(I, 3, 3)], Float64[])

    y = ones(Float64, length(b)) * -1 * (1+sqrt(17)) / 4

    s = c - evaluate(A, y)

    return problem, PointPrimalDual(x, y, s)
end

function testcase1ater()
    cmat = [sparse([1, 2, 3], [1, 2, 3], [1., 1., 1.])]
    c = PointE(cmat, Vector{Float64}())

    a1mats = [sparse([1], [1], [1.], 3, 3)]
    a1 = PointE(a1mats, Vector{Float64}())

    a2mats = [sparse([2, 3], [2, 1], [1., 1.], 3, 3)]
    a2 = PointE(a2mats, Vector{Float64}())

    a3mats = [sparse([2, 3], [1, 3], [1., 1.], 3, 3)]
    a3 = PointE(a3mats, Vector{Float64}())

    A = [a1, a2, a3]

    b = Vector{Float64}([1, 1, 1])

    problem = SDCOContext(c, A, b)

    x = PointE([Matrix{Float64}(I, 3, 3)], Float64[])

    y = Float64[0, -1, -0.5]

    s = c - evaluate(A, y)

    return problem, PointPrimalDual(x, y, s)
end

function testcase1b()
    c = PointE(Dense{Float64}[], Float64[1., 1.])

    A = [PointE(Dense{Float64}[], Float64[1., 2.])]
    b = [1.]

    problem = SDCOContext(c, A, b)

    x = PointE(Dense{Float64}[], Float64[1/3, 1/3])
    y = [0.]
    s = PointE(Dense{Float64}[], Float64[1., 1.])
    return problem, PointPrimalDual(x, y, s)
end

function testcase1c(; symmetric=false)
    if symmetric
        cmat = [Symmetric(sparse(1.0I, 3, 3), :L), Symmetric(sparse(1.0I, 3, 3), :L)]

        mat1 = Symmetric(sparse([1], [1], [1.], 3, 3), :L)
        mat2 = Symmetric(sparse([2, 3], [2, 1], [1., 1.], 3, 3), :L)
        mat3 = Symmetric(sparse([2, 3], [1, 3], [1., 1.], 3, 3), :L)
        matnull = Symmetric(spzeros(3,3), :L)
    else
        cmat = [sparse(1.0I, 3, 3), sparse(1.0I, 3, 3)]

        mat1 = sparse([1], [1], [1.], 3, 3)
        mat2 = sparse([1, 2, 3], [3, 2, 1], [1., 1., 1.], 3, 3)
        mat3 = sparse([1, 2, 3], [2, 1, 3], [1., 1., 1.], 3, 3)
        matnull = spzeros(3,3)
    end

    c = PointE(cmat, Float64[1., 1.])

    a1 = PointE([mat1, matnull], Float64[0, 0])
    a2 = PointE([mat2, matnull], Float64[0, 0])
    a3 = PointE([mat3, matnull], Float64[0, 0])

    a4 = PointE([matnull, mat1], Float64[0, 0])
    a5 = PointE([matnull, mat2], Float64[0, 0])
    a6 = PointE([matnull, mat3], Float64[0, 0])

    a7 = PointE([matnull, matnull], Float64[1., 2.])
    A = [a1, a2, a3, a4, a5, a6, a7]

    b = Vector{Float64}([1, 1, 1, 1, 1, 1, 1])

    problem = SDCOContext(c, A, b)

    y0 = -(1+sqrt(17)) / 4
    y = y0 * ones(7)

    x = PointE([Matrix{Float64}(I, 3, 3), Matrix([1 0 0.25; 0 0.5 0; 0.25 0 1])], Float64[(5-sqrt(17))/2, (-3+sqrt(17))/4])

    s = c - evaluate(A, y)

    return problem, PointPrimalDual(x, y, s)
end

function testcase2(p, q, r; storage=:sparsefull)

    ## Objective
    B0 = rand(q, r)

    cmat::Sparse{Float64} = spzeros(q+r, q+r)
    if storage == :sparsesym
        for i=1:r, j=1:q
            cmat[q+i, j] = B0[j, i]
        end
        c = PointE([Symmetric(cmat, :L)], Float64[])
    elseif storage == :sparsefull
        for i=1:r, j=1:q
            cmat[q+i, j] = B0[j, i]
            cmat[j, i+q] = B0[j, i]
        end
        c = PointE([cmat], Float64[])
    else
        error("Unknown parameter $storage, choose from :sparsesym, :sparsefull")
    end

    if storage == :sparsesym
        A = PointE{Float64, SparseSym{Float64}}[]
    elseif storage == :sparsefull
        A = PointE{Float64, Sparse{Float64}}[]
    end
    b = Float64[]

    if storage == :sparsesym
        push!(A, PointE([Symmetric(sparse(-1.0I, q+r, q+r), :L)], Float64[]))
    elseif storage == :sparsefull
        push!(A, PointE([sparse(-1.0I, q+r, q+r)], Float64[]))
    end
    push!(b, -1.)

    for ctrind=1:p
        aimat = spzeros(q+r, q+r)
        Bi = rand(q, r)
        if storage == :sparsesym
            for i=1:r, j=1:q
                aimat[q+i, j] = Bi[j, i]
            end
            push!(A, PointE([Symmetric(aimat, :L)], Float64[]))
        elseif storage == :sparsefull
            for i=1:r, j=1:q
                aimat[q+i, j] = Bi[j, i]
                aimat[j, q+i] = Bi[j, i]
            end
            push!(A, PointE([aimat], Float64[]))
        end
        push!(b, 0.)
    end

    problem = SDCOContext(c, A, b)

    y = zeros(p+1)
    y[1] = norm(B0) + 1
    x = PointE([1/(q+r) * Matrix(1.0I, q+r, q+r)], Float64[])
    s = c - evaluate(A, y)

    z = PointPrimalDual(x, y, s)

    return problem, z
end

function testcase4a()
    cmat = [sparse([3], [3], [1.], 3, 3)]
    a1mats = [sparse([1, 2, 3], [2, 1, 3], [1., 1., 1.], 3, 3)]
    a2mats = [sparse([2], [2], [1.], 3, 3)]

    c = PointE(cmat, Vector{Float64}())
    a1 = PointE(a1mats, Vector{Float64}())
    a2 = PointE(a2mats, Vector{Float64}())

    A = [a1, a2]

    b = Vector{Float64}([-1, 0])

    return SDCOContext(c, A, b)
end

function testcase4b()
    cmat = [sparse(Float64[0 1; 1 0])]
    a1mats = [sparse(Float64[-1 0; 0 0])]
    a1mats = [sparse(Float64[0 0; 0 -1])]

    c = PointE(cmat, Vector{Float64}())
    a1 = PointE(a1mats, Vector{Float64}())
    a2 = PointE(a2mats, Vector{Float64}())

    A = [a1, a2]

    b = Vector{Float64}([-1, 0])

    return SDCOContext(c, A, b)
end

function testcase4c()
    cmat = [sparse(Float64[0 0; 0 0])]
    a1mats = [sparse(Float64[1 0; 0 0])]
    a1mats = [sparse(Float64[0 1; 1 0])]

    c = PointE(cmat, Vector{Float64}())
    a1 = PointE(a1mats, Vector{Float64}())
    a2 = PointE(a2mats, Vector{Float64}())

    A = [a1, a2]

    b = Vector{Float64}([0, 2])

    return SDCOContext(c, A, b)
end


text:
using Ju
using BenchmarkTools

function run_env(env, steps=10^3)
    reset!(env)
    while steps > 0
        if observe(env).isdone
            reset!(env)
        end
        env(sample(actionspace(env)))
        steps -= 1
    end
end

function benchmark_env(env, nsteps=10^3)
    println("\n", repeat('=', 50))
    println("\nreset!\n")
    display(@benchmark reset!($env))
    println("\n", repeat('=', 50))
    println("\nobserve\n")
    display(@benchmark observe($env))
    println("\n", repeat('=', 50))
    println("\nrun $nsteps steps\n")
    display(@benchmark run_env($env, $nsteps))
    println("\n", repeat('=', 50))
end

benchmark_env(CartPoleEnv())

# julia --project=benchmarks -e 'include("benchmarks/environments.jl")'

# ==================================================

# reset!

# BenchmarkTools.Trial:
#   memory estimate:  32 bytes
#   allocs estimate:  1
#   --------------
#   minimum time:     42.047 ns (0.00% GC)
#   median time:      43.184 ns (0.00% GC)
#   mean time:        53.883 ns (13.18% GC)
#   maximum time:     47.807 μs (99.89% GC)
#   --------------
#   samples:          10000
#   evals/sample:     991
# ==================================================

# observe

# BenchmarkTools.Trial:
#   memory estimate:  32 bytes
#   allocs estimate:  1
#   --------------
#   minimum time:     8.863 ns (0.00% GC)
#   median time:      11.440 ns (0.00% GC)
#   mean time:        19.778 ns (37.83% GC)
#   maximum time:     51.952 μs (99.97% GC)
#   --------------
#   samples:          10000
#   evals/sample:     999
# ==================================================

# run 1000 steps

# BenchmarkTools.Trial:
#   memory estimate:  31.25 KiB
#   allocs estimate:  1000
#   --------------
#   minimum time:     73.355 μs (0.00% GC)
#   median time:      74.090 μs (0.00% GC)
#   mean time:        86.768 μs (8.98% GC)
#   maximum time:     52.824 ms (99.83% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1
# ==================================================


text:
@testset "Torsion points" begin

  E43_a1 = @inferred EllipticCurve([0, 1, 1, 0, 0])

  curves_to_test_tor_short = [
  ([1, 2], 4),
  ([2, 3], 2),
  ([0, 1], 6),
  ([1, 0], 2),
  ([1, 1], 1)
  ]

  curves_to_test_tor_struc =
  [
  ([0, -1, 1, -7820, -263580], [1]),
  ([1, 0, 1, -2731, -55146], [2]),
  ([0, 1, 1, -9, -15], [3]),
  ([1, 1, 1, -80, 242], [4]),
  ([0, -1, 1, -10, -20], [5]),
  ([1, 0, 1, -36, -70], [6]),
  ([1, -1, 1, -3, 3], [7]),
  ([1, 1, 1, 35, -28], [8]),
  ([1, -1, 1, -14, 29], [9]),
  ([1, 0, 0, -45, 81], [10]),
  ([1, -1, 1, -122, 1721], [12]),
  ([1, 1, 1, -135, -660], [2, 2]),
  ([1, 1, 1, -10, -10], [2, 4]),
  ([1, 0, 1, -19, 26], [2, 6]),
  ([1, 0, 0, -1070, 7812], [2, 8]),
  ]

  Rx, x = PolynomialRing(QQ, "x")
  K1, a = number_field(x^2 - x - 1)
  K2, a = number_field(x^2 +2)

  curves_to_test_tor_struc_nf =
   [(map(K1, [1, 1, 1, -3, 1]), [15]),
    (map(K2, [1, 0, 0, 115, 561]), [10, 2])
   ]

  @testset "Point order computation" begin
    E = EllipticCurve([0, -1, 1, -7820, -263580])
    @test 1 == @inferred order(infinity(E))

    E = EllipticCurve([1, 0, 1, -2731, -55146])
    @test 2 == @inferred order(E([fmpq(-121, 4), fmpq(117, 8)]))

    E = EllipticCurve([0, 1, 1, -9, -15])
    @test 3 == @inferred order(E([5, 9]))

    E = EllipticCurve([1, 1, 1, -80, 242])
    @test 4 == @inferred order(E([5, -2]))

    E = EllipticCurve([0, -1, 1, -10, -20])
    @test 5 == @inferred order(E([5, 5]))

    E = EllipticCurve([1, 0, 1, -36, -70])
    @test 6 == @inferred order(E([-4, 5]))

    E = EllipticCurve([1, -1, 1, -3, 3])
    @test 7 == @inferred order(E([1, 0]))

    E = EllipticCurve([1, 1, 1, 35, -28])
    @test 8 == @inferred order(E([2, 6]))

    E = EllipticCurve([1, -1, 1, -14, 29])
    @test 9 == @inferred order(E([-3, 7]))

    E = EllipticCurve([1, 0, 0, -45, 81])
    @test 10 == @inferred order(E([0, 9]))

    E = EllipticCurve([1, -1, 1, -122, 1721])
    @test 12 == @inferred order(E([-9, 49]))
  end

  @testset "Torsion test" begin
    E = EllipticCurve([1, -1, 1, -19353, 958713])
    @test @inferred is_torsion_point(E([103,172]))
    @test @inferred !is_torsion_point(E([-121,1292]))
  end

  @testset "Torsion points (Lutz-Nagell)" begin
    for c in curves_to_test_tor_short
      E = EllipticCurve(c[1])
      T = @inferred torsion_points_lutz_nagell(E)
      @test c[2] == length(T)
    end
  end

  @testset "Torsion points (division polynomials" begin
    for c in curves_to_test_tor_short
      E = EllipticCurve(c[1])
      T = @inferred torsion_points_division_poly(E)
      @test c[2] == length(T)
    end

    for c in curves_to_test_tor_struc
      E = EllipticCurve(c[1])
      T = @inferred torsion_points_division_poly(E)
      @test prod(c[2]) == length(T)
    end
  end

  @testset "Torsion points" begin
    for c in curves_to_test_tor_short
      E = EllipticCurve(c[1])
      T = @inferred torsion_points(E)
      @test c[2] == length(T)
      # test the caching
      T = @inferred torsion_points(E)
      @test c[2] == length(T)
    end

    for c in curves_to_test_tor_struc
      E = EllipticCurve(c[1])
      T = @inferred torsion_points(E)
      @test prod(c[2]) == length(T)
      # test the caching
      T = @inferred torsion_points(E)
      @test prod(c[2]) == length(T)
    end
  end

  @testset "Torsion point structure" begin
    for c in curves_to_test_tor_struc
      E = EllipticCurve(c[1])
      T = @inferred torsion_structure(E)
      @test T[1] == c[2]
      @test order(T[2][1]) == T[1][1]
      if length(T[1]) == 2
        @test order(T[2][2]) == T[1][2]
      end
      # test the caching
      T = @inferred torsion_structure(E)
      @test T[1] == c[2]
      @test order(T[2][1]) == T[1][1]
      if length(T[1]) == 2
        @test order(T[2][2]) == T[1][2]
      end
    end
  end

  @testset "Torsion point structure number fields" begin
    for c in curves_to_test_tor_struc_nf
      E = EllipticCurve(c[1])
      T = @inferred torsion_structure(E)
      @test T[1] == c[2]
      @test order(T[2][1]) == T[1][1]
      if length(T[1]) == 2
        @test order(T[2][2]) == T[1][2]
      end
      # test the caching
      T = @inferred torsion_structure(E)
      @test T[1] == c[2]
      @test order(T[2][1]) == T[1][1]
      if length(T[1]) == 2
        @test order(T[2][2]) == T[1][2]
      end
    end
  end

  @testset "Division polynomial" begin

    Kx, x = PolynomialRing(base_field(E43_a1), "x")
    Kxy, y = PolynomialRing(Kx, "y")
    f = @inferred division_polynomial(E43_a1, 4, x, y)
    @test  f == (4*x^6 + 8*x^5 + 20*x^3 + 20*x^2 + 8*x - 2)*y + 2*x^6 + 4*x^5 + 10*x^3 + 10*x^2 + 4*x - 1

  end
end


text:
# A quick summmed-area table implementation
function sumtable!(table, arr)
    out = arr
    for n in 1:ndims(table)
        out = cumsum!(table, out, dims = n)
    end
    out
end

"Efficient region queries over a summed-area table. Note that
the `table` should be padded by one strip of zeros on the left
and top sides."
function regionquery end

function regionquery(table::Vector, ylo, yhi)
    table[yhi+1] - table[ylo]
end

function regionquery(table::Matrix, ylo, yhi, xlo, xhi)
    yhi += 1
    xhi += 1
    a = table[yhi, xhi]
    b = table[ylo, xhi]
    c = table[yhi, xlo]
    d = table[ylo, xlo]
    a - b - c + d
end

"Specify how to handle boundaries at the edges of the image"
abstract type EdgeCondition end
"Clamp causes smaller blurs to be performed at the edges"
struct Clamp <: EdgeCondition end
"Zero assumes implicit zeroes outside of the image region"
struct Zero  <: EdgeCondition end

#
# I don't like all this but it's the best I've got so far for
# handling multidimensional inputs efficiently.
# CartesianIndexes were slower than otherwise, and I didn't want
# to use generated functions.
#

@inline scale(::Zero,  r, ylo, yhi) = 1 / (2r + 1)^1
@inline scale(::Zero,  r, ylo, yhi, xlo, xhi) = 1 / (2r + 1)^2

@inline scale(::Clamp, r, ylo, yhi, xlo=0, xhi=0) =
    1 / ((yhi-ylo+1) * (xhi-xlo+1))

abstract type Blur end

struct BoxBlur{N} <: Blur
    table::Array{Float64, N}
    tableview::SubArray
    function BoxBlur{N}(sz) where {N}
        # Create a sum table padded with zeros on the starting side
        table = zeros(map(x->x+1, sz))
        tableview = view(table, map(x -> (1:x) .+ 1, sz)...)
        new(table, tableview)
    end
end
BoxBlur(sz) = BoxBlur{length(sz)}(sz)

# Preparing a blur happens once per image. We can subsequently
# blur the same image efficiently at multiple scales.
prepare!(b::BoxBlur, arr) = sumtable!(b.tableview, arr)

function blur!(b::BoxBlur{2}, out, r, edge)
    (sy, sx) = size(out)
    for x in 1:sx, y in 1:sy
        ylo, yhi = max(1, y-r), min(y+r, sy)
        xlo, xhi = max(1, x-r), min(x+r, sx)
        s = scale(edge, r, ylo, yhi, xlo, xhi)
        out[y, x] = s * regionquery(b.table, ylo, yhi, xlo, xhi)
    end
    out
end

function blur!(b::BoxBlur{1}, out, r, edge)
    (sy,) = size(out)
    for y in 1:sy
        ylo, yhi = max(1, y-r), min(y+r, sy)
        s = scale(edge, r, ylo, yhi)
        out[y] = s * regionquery(b.table, ylo, yhi)
    end
    out
end

struct IteratedBoxBlur{N} <: Blur
    box::BoxBlur
    target::Array{Float64, N}
    iters::Int
    IteratedBoxBlur{N}(sz, iters) where {N} = new(BoxBlur(sz), zeros(sz), iters)
end
IteratedBoxBlur(sz, iters) = IteratedBoxBlur{length(sz)}(sz, iters)

prepare!(b::IteratedBoxBlur, arr) = copy!(b.target, arr)

function blur!(b::IteratedBoxBlur, out, r, edge)
    copy!(out, b.target)
    for _ in 1:b.iters
        prepare!(b.box, out)
        blur!(b.box, out, r, edge)
    end
    out
end


text:
using MeshingBenchmarks
using Test

MeshingBenchmarks.benchmark()


text:
using ProfileView
using PolynomialBasis
using ImplicitDomainQuadrature
using BenchmarkTools
# using Revise
using CutCellDG
include("../../../test/useful_routines.jl")
# include("transformation-elasticity-solver.jl")

# TES = TransformationElasticitySolver

polyorder = 2
nelmts = 15
penaltyfactor = 1e3

################################################################################
width = 1.0              # mm
interfacecenter = [0.5, 0.5]
interfaceradius = 0.4
outerradius = 1.2
################################################################################

################################################################################
K1, K2 = 247.0, 192.0    # GPa
mu1, mu2 = 126.0, 87.0   # GPa
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
theta0 = -0.067
transfstress =
    CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)
################################################################################

meshwidth = [width, width]
numqp = required_quadrature_order(polyorder) + 2

levelsetbasis = HermiteTensorProductBasis(2)
quad = tensor_product_quadrature(2, 4)
elasticitybasis = LagrangeTensorProductBasis(2, polyorder)

basispts = interpolation_points(elasticitybasis)
dim, numpts = size(interpolation_points(levelsetbasis))
cgmesh = CutCellDG.CGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], numpts)
dgmesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basispts)

levelset = CutCellDG.LevelSet(
    x -> -circle_distance_function(x, interfacecenter, interfaceradius)[1],
    cgmesh,
    levelsetbasis,
    quad,
)

elementsize = CutCellDG.element_size(cgmesh)
minelmtsize = minimum(elementsize)
maxelmtsize = maximum(elementsize)
penalty = penaltyfactor / minelmtsize * 0.5 * (lambda1 + mu1 + lambda2 + mu2)
tol = minelmtsize^(polyorder + 1)
boundingradius = 1.5 * maxelmtsize

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)

@profview cellquads = CutCellDG.CellQuadratures(cutmesh,levelset,numqp)


text:
grid_size = 25
gamma_vals = linspace(0.05, 0.95, grid_size)
w_bar_vals = similar(gamma_vals)

mcm = McCallModel()

for (i, gamma) in enumerate(gamma_vals)
    mcm.gamma = gamma
    w_bar = compute_reservation_wage(mcm)
    w_bar_vals[i] = w_bar
end

plot(gamma_vals,
    w_bar_vals,
    lw=2,
    alpha=0.7,
    xlabel="job offer rate",
    ylabel="reservation wage",
    label=L"$\bar w$ as a function of $\gamma$")


text:
function rrule(::Type{ITensor}, x1::AllowAlias, x2::TensorStorage, x3)
  y = ITensor(x1, x2, x3)
  function ITensor_pullback(ȳ)
    x̄1 = NoTangent()
    x̄2 = ȳ.tensor.storage
    x̄3 = NoTangent()
    return (NoTangent(), x̄1, x̄2, x̄3)
  end
  return y, ITensor_pullback
end

function rrule(::Type{ITensor}, x1::AllowAlias, x2::Tensor)
  y = ITensor(x1, x2)
  function ITensor_pullback(ȳ)
    x̄1 = NoTangent()
    x̄2 = Tensor(x1, ȳ)
    return (NoTangent(), x̄1, x̄2)
  end
  return y, ITensor_pullback
end

function rrule(::Type{ITensor}, x1::Tensor)
  y = ITensor(x1)
  function ITensor_pullback(ȳ)
    x̄1 = Tensor(ȳ)
    return (NoTangent(), x̄1)
  end
  return y, ITensor_pullback
end

function rrule(::typeof(itensor), x1::Tensor)
  y = itensor(x1)
  function itensor_pullback(ȳ)
    x̄1 = tensor(ȳ)
    return (NoTangent(), x̄1)
  end
  return y, itensor_pullback
end

function rrule(f::Type{<:Tensor}, x1::AllowAlias, x2::ITensor)
  y = f(x1, x2)
  function Tensor_pullback(ȳ)
    x̄1 = NoTangent()
    x̄2 = ITensor(x1, ȳ)
    return (NoTangent(), x̄1, x̄2)
  end
  return y, Tensor_pullback
end

function rrule(f::Type{<:Tensor}, x1::ITensor)
  y = f(x1)
  function Tensor_pullback(ȳ)
    x̄1 = ITensor(ȳ)
    return (NoTangent(), x̄1)
  end
  return y, Tensor_pullback
end

function rrule(::typeof(tensor), x1::ITensor)
  y = tensor(x1)
  function tensor_pullback(ȳ)
    x̄1 = ITensor(typeof(storage(x1))(ȳ.storage.data), inds(x1))
    return (NoTangent(), x̄1)
  end
  return y, tensor_pullback
end

# Special case for contracting a pair of ITensors
function ChainRulesCore.rrule(::typeof(contract), x1::ITensor, x2::ITensor)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1, x̄2)
  end
  return y, contract_pullback
end

@non_differentiable ITensors.optimal_contraction_sequence(::Any)

function ChainRulesCore.rrule(::typeof(*), x1::Number, x2::ITensor)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1[], x̄2)
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(*), x1::ITensor, x2::Number)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1, x̄2[])
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(+), x1::ITensor, x2::ITensor)
  y = x1 + x2
  function add_pullback(ȳ)
    return (NoTangent(), ȳ, ȳ)
  end
  return y, add_pullback
end

function ChainRulesCore.rrule(::typeof(itensor), x::Array, a...)
  y = itensor(x, a...)
  function itensor_pullback(ȳ)
    uȳ = permute(unthunk(ȳ), a...)
    x̄ = reshape(array(uȳ), size(x))
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, itensor_pullback
end

function ChainRulesCore.rrule(::Type{ITensor}, x::Array{<:Number}, a...)
  y = ITensor(x, a...)
  function ITensor_pullback(ȳ)
    # TODO: define `Array(::ITensor)` directly
    uȳ = Array(unthunk(ȳ), a...)
    x̄ = reshape(uȳ, size(x))
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, ITensor_pullback
end

function ChainRulesCore.rrule(::Type{ITensor}, x::Number)
  y = ITensor(x)
  function ITensor_pullback(ȳ)
    x̄ = ȳ[]
    return (NoTangent(), x̄)
  end
  return y, ITensor_pullback
end

function ChainRulesCore.rrule(::typeof(dag), x)
  y = dag(x)
  function dag_pullback(ȳ)
    x̄ = dag(unthunk(ȳ))
    return (NoTangent(), x̄)
  end
  return y, dag_pullback
end

function ChainRulesCore.rrule(::typeof(permute), x::ITensor, a...)
  y = permute(x, a...)
  function permute_pullback(ȳ)
    x̄ = permute(unthunk(ȳ), inds(x))
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, permute_pullback
end

@non_differentiable combiner(::Indices)


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using RDatasets, Plots, Measures
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.24"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data = dataset("datasets", "iris")
	Text("Number of rows: $(nrow(data))")
end

# ╔═╡ 66321bda-250f-11eb-1fbd-df6e6c80e7cc
insertSpace(name) = begin
    i = findlast(isuppercase,name)
    name[1:i-1]*" "*name[i:end]
end

# ╔═╡ 663255dc-250f-11eb-2124-3b692cefccb8
begin
	featureNames = insertSpace.(string.(names(data)))[1:4]
	("Names of features:", featureNames)
end

# ╔═╡ 6632ece0-250f-11eb-337d-d9c2482dfc73
begin
	speciesNames = unique(data.Species)
	speciesFreqs = [sn => sum(data.Species .== sn) for sn in speciesNames]
	Text("Frequency per species: $(speciesFreqs)")
end

# ╔═╡ 66403c4c-250f-11eb-2b72-17cdf591ee90
begin
	default(msw = 0, ms = 3)

	scatters = [
		scatter(data[:,i], data[:,j], c=[:blue :red :green], group=data.Species,
			xlabel=featureNames[i], ylabel=featureNames[j], legend = i==1 && j==1)
		for i in 1:4, j in 1:4 ]

	plot(scatters..., size=(1200,800), margin = 4mm)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.24"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═66321bda-250f-11eb-1fbd-df6e6c80e7cc
# ╠═663255dc-250f-11eb-2124-3b692cefccb8
# ╠═6632ece0-250f-11eb-337d-d9c2482dfc73
# ╠═66403c4c-250f-11eb-2b72-17cdf591ee90
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
"""
    import_leadfield()

Return a leadfield.

Simple wrapper function that calls lower level functions and returns a valid leadfield.

#### Example

```julia
ldf = import_leadfield()
```

"""
function import_leadfield(;verbose::Bool=false)

    L                     = readLFT(verbose=verbose)
    positions, neighbours = readLOC(verbose=verbose)
    names, a, b           = readELP(verbose=verbose)

    L = convert(Array{Float64, 3}, L)

    if size(L, 1) != size(positions,2); error("Mismatched leadfield and locations"); end
    if size(L, 2) != size(positions,1); error("Mismatched leadfield and locations"); end
    if size(L, 3) != size(names,1); error("Mismatched leadfield and names"); end

    Leadfield(L, vec(positions[1,:]), vec(positions[2,:]), vec(positions[3,:]), vec(names))
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module MathProgBase

include(joinpath(dirname(@__FILE__),"SolverInterface","SolverInterface.jl"))
using .SolverInterface
using Compat

# deprecated name
if VERSION >= v"0.4-rc2"
    @eval @Base.deprecate_binding MathProgSolverInterface SolverInterface
else
    const MathProgSolverInterface = SolverInterface
end

include("defaultsolvers.jl")

include(joinpath(dirname(@__FILE__),"HighLevelInterface","HighLevelInterface.jl"))
using .HighLevelInterface
export linprog, mixintprog, quadprog, buildlp, solvelp

end


text:
include("euler/euler.jl")
using .Calculus: fibonacci_index, fibonacci_numbers

function compute(n::Int)::Int
    index = fibonacci_index(n)
    fibonacci = fibonacci_numbers(index, Int)
    last_sum, new_sum = 0, 0
    result = sum(fibonacci[1:5]) - 1
    for i ∈ 6:index + 1
        last_sum, new_sum = new_sum, new_sum + last_sum + fibonacci[i - 3]
        result += new_sum + fibonacci[i]
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
##### Beginning of file

Test.@test( isdir(RemoveLFS.package_directory()) )

Test.@test( isdir(RemoveLFS.package_directory("ci")) )

Test.@test( isdir(RemoveLFS.package_directory("ci", "travis")) )

Test.@test( isdir(RemoveLFS.package_directory(TestModuleA)) )

Test.@test( isdir(RemoveLFS.package_directory(TestModuleB)) )

Test.@test(
    isdir( RemoveLFS.package_directory(TestModuleB, "directory2",) )
    )

Test.@test(
    isdir(
        RemoveLFS.package_directory(
            TestModuleB, "directory2", "directory3",
            )
        )
    )

Test.@test_throws(
    ErrorException,RemoveLFS.package_directory(TestModuleC),
    )

##### End of file


text:
include("core/core.jl")
include("integration/integration.jl")
include("grammar/grammar.jl")
include("repl.jl")

using .Context;

println(INFO_STYLE, "Welcome! Use EQN format for input.")
println(INFO_STYLE, "Enter ? for help.")

if !isempty(ARGS)
    for arg in ARGS
        file(arg)
    end
end

while true
    repl()
end


text:
function toString(b::SymbolicType)
    b = Basic(b)
    a = ccall((:basic_str_julia, libsymengine), Cstring, (Ptr{Basic}, ), &b)
    string = unsafe_string(a)
    ccall((:basic_str_free, libsymengine), Void, (Cstring, ), a)
    return string
end

Base.show(io::IO, b::SymbolicType) = print(io, toString(b))

" show symengine logo "
type AsciiArt x end
function ascii_art()
    out = ccall((:ascii_art_str, libsymengine),  Ptr{UInt8},  ())
    AsciiArt(unsafe_string(out))
end

Base.show(io::IO, x::AsciiArt) = print(io, x.x)


text:
@init @require CUDAnative="be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative

    @Base.propagate_inbounds function __read_only_load(::CUDA, A::CUDAnative.CuDeviceArray, index::Integer)
        CUDAnative.ldg(A, index)
    end
end

@Base.propagate_inbounds __read_only_load(::CPU, A::Array, index::Integer) = A[index]

"""
    read_only_load(A, i)

Index the array `A` with the linear index `i`.  On the GPU this uses the
read-only texture cache (e.g., via `ldg` in CUDAnative) so the memory `A`
refers to must not be written to during the kernel execution.

See also: `Base.getindex`, `CUDAnative.ldg`
"""
@Base.propagate_inbounds read_only_load(A, index) = __read_only_load(backend(), A, index)

max_stars_repo_name: vchuravy/GPUifyLoops

max_issues_repo_name: vchuravy/GPUifyLoops

max_forks_repo_name: vchuravy/GPUifyLoops


text:
using LinearAlgebra
include("../src/EDKit.jl")
using .EDKit
using Test

@testset "Basic" begin
    L = 10
    b1 = TranslationalBasis(x -> sum(x) == 2, 0, L)
    b2 = TranslationalBasis(x -> sum(x) == 0, 1, L)
    b = DoubleBasis(b1, b2)
    mat = rand(2,2) |> Hermitian
    H = trans_inv_operator(mat, 1, b) |> Array
    @test size(H) == (5, 0)
end

@testset "Spin-1 Random" begin
    L = 6
    Sp = begin
        p3 = [1]
        p2 = [2, 4, 10]
        p1 = [3, 5, 7, 11, 13, 19]
        z0 = [6, 8, 12, 14, 16, 20, 22]
        m1 = [9, 15, 17, 21, 23, 25]
        m2 = [18, 24, 26]
        m3 = [27]
        mat = zeros(ComplexF64, 27, 27)
        mat[m2, m3] .= rand(ComplexF64, 3, 1)
        mat[m1, m2] .= rand(ComplexF64, 6, 3)
        mat[z0, m1] .= rand(ComplexF64, 7, 6)
        mat[p1, z0] .= rand(ComplexF64, 6, 7)
        mat[p2, p1] .= rand(ComplexF64, 3, 6)
        mat[p3, p2] .= rand(ComplexF64, 1, 3)
        mat
    end
    XY = spin((1, "+-"), (1, "-+"), (1, "z1"), (1, "1z"), D=3)
    A = 0.0
    for n = 1:2L
        b1 = ProjectedBasis(x -> sum(x)==(n-1), L, base=3)
        b2 = ProjectedBasis(x -> sum(x)==n, L, base=3)
        basis = DoubleBasis(b1, b2)
        e1, v1 = trans_inv_operator(XY, 2, b1) |> Array |> Hermitian |> eigen
        e2, v2 = trans_inv_operator(XY, 2, b2) |> Array |> Hermitian |> eigen
        op = trans_inv_operator(Sp, 3, basis) |> Array
        A += abs2.(v1' * op * v2) |> sum
    end
    e, v = trans_inv_operator(XY, 2, L) |> Array |> Hermitian |> eigen
    op = trans_inv_operator(Sp, 3, L) |> Array
    B = abs2.(v' * op * v) |> sum
    @test A ≈ B atol = 1e-7
end


text:
using Serialization
# using DataStructures
# const List = LinkedList

@nospecialize
include("compact.jl")

function _typed_list(::Type{T}, args::T...) where {T}
    foldr(args, init = nil(T)) do e, last
        cons(e, last)
    end
end

function compress_impl!(xs::List{T}) where {T}
    args = Any[]
    maxcomplexity = 0
    for x in collect(xs)
        meta, arg = compress_impl!(x)
        push!(args, arg)
        maxcomplexity += meta.complexity
    end
    encoded = Call(Constructor{_typed_list}(), (compress(T), args...))
    may_cache_call(SimpleMeta(maxcomplexity), encoded)
end

struct Buf{N}
    units::NTuple{N,UInt8}
end

Base.show(io::IO, buf::Buf{N}) where {N} = Base.show(io, "Buf{$N}()")

struct TypeLevel{T,BufData} end

function from_type(::Type{TypeLevel{T,BufData}}) where {T,BufData}
    compressed = deserialize(IOBuffer(UInt8[BufData.units...]))
    decompress(compressed)::T
end

function to_type(x::T) where {T}
    io = IOBuffer()
    serialize(io, compress(x))
    seek(io, 0)
    TypeLevel{T,Buf(Tuple(take!(io)))}
end

function Base.show(io::IO, ::T) where {T<:TypeLevel}
    print(io, "TypeEncoding(")
    Base.show(io, from_type(T))
    print(io, ")")
end

# xs = to_type(
#     quote
#         f
#         $(list(1, 2, 3))
#     end
# )()

# println(xs)
# println(CallPool)
# println(RefValPool)
# println(ExprPool)

@specialize


text:
include("...\\..\\..\\..\\..\\add_load_paths.jl")

using StaticNets, DynNets
using PyPlot
pygui(true)
using JLD2

import Utilities:meanSq, maxLargeVal

@load("./data/congress_covoting_US/juliaEstimates.jld", estParSS_T, obsMat_T,changeStats_T, stats_T)

#estParSS_T = reshape(estParSS_T, 1, 2, 73)

# Estimate SD-PML-ERGM of the paper

indTvPar = BitArray([false, true])
indTargPar =indTvPar #BitArray([true,true])

obsT = changeStats_T #[Real.(changeStats_t) for changeStats_t in changeStats_T]

model = DynNets.SdErgmDirBinGlobalPseudo(obsT, DynNets.fooGasPar, indTvPar,"")

estPar,convFlag,UM,startPoint,hess_opt = DynNets.estimate(model;UM =meanSq(estParSS_T[:,:],1), indTvPar = indTvPar, indTargPar = indTargPar,hess_opt_flag = true)

targetAllTv = false
 gasParSampled  = Array{Array{Float64,2}}(undef, Nsample)
 flag_cov_issue = falses(Nsample)
 gasParEst= fill(fill(Float64[],2), Nsample)
 startPointEst = fill(fill(0.0,2), Nsample)
 staticEst = fill(fill(0.0,2), Nsample)
 filtPar_T_Nsample = zeros(Nterms,T,Nsample)
 pVals_Nsample = zeros(Nterms,Nsample)
 scoreAutoc_Nsample = zeros(Nterms,Nsample)
 convFlag = falses(Nsample)

T_train = size(obsT)[1]

restrCase = false
if restrCase
    # restricted parameters
    f_T(x) = DynNets.logLike_T(model::DynNets.SdErgmDirBinGlobalPseudo,
                                obsT, x, model.indTvPar)
    allpar0 = DynNets.array_2_vec_all_par(model, estPar, model.indTvPar)
    vec_of_f_t =[x -> DynNets.logLike_t(model::DynNets.SdErgmDirBinGlobalPseudo,
                                        obsT[1:t], x, model.indTvPar)
                 for t in 1:T_train]
else
    #Not restricted parameters
    f_T(x) = DynNets.logLike_T(model::DynNets.SdErgmDirBinGlobalPseudo,
                                obsT,DynNets.restrictGasPar(model, x, model.indTvPar),
                                model.indTvPar)
    allpar0 = DynNets.unRestrictGasPar(model,
                                        DynNets.array_2_vec_all_par(model, estPar, model.indTvPar),
                                        indTvPar )
    vec_of_f_t =[x -> DynNets.logLike_t(model::DynNets.SdErgmDirBinGlobalPseudo,
                                        obsT[1:t],
                                        DynNets.restrictGasPar(model, x, model.indTvPar),
                                         model.indTvPar)
                                        for t in 1:T_train]
end

f_T(allpar0)
covA0 =  ForwardDiff.hessian(f_T, allpar0)./T_train
# run once to precompile the functions
[f(allpar0) for f in vec_of_f_t]
vec_g_t = [ForwardDiff.gradient(f, allpar0) for f in vec_of_f_t[2:end]]
global covB0 = zeros(size(covA0))
for g_t in vec_g_t
    global covB0 += g_t * g_t'
end
global covB0 = covB0./ T_train
diagCorrect = 1e-3
covMatHat = pinv(covA0) * covB0 * pinv(covA0)
covMatHat = Symmetric(covMatHat)
covMatHat_corr = (covMatHat .+ Diagonal(diagCorrect.*ones(size(covMatHat)[1])))
parDistr = MvNormal(allpar0, diag(covMatHat_corr))
samp = rand(parDistr, 6000)
indsToKeep =    (0 .< samp[end-1, :] .< 1) .& (samp[end, :] .>0) .&
                (0 .< samp[2, :] .< 1) .& (samp[3, :] .>0)
sum(indsToKeep)
hist(samp[end,:])
N_samp_par = 1000
sampled_gas_par = samp[:, indsToKeep][:, 1:N_samp_par]

gasParSampled[n] = sampled_gas_par
    gasParEst[n] = estPar # store gas parameters
    startPointEst[n] = startPoint
    gasParVec = zeros(sum(indTvPar)*3); for i=0:(sum(indTvPar)-1) gasParVec[1+3i : 3(i+1)] = estPar[indTvPar][i+1]; end
    constParVec = zeros(sum(.!indTvPar)); for i=1:sum(.!indTvPar) constParVec[i] = estPar[.!indTvPar][i][1]; end
    gasFiltPar , pseudolike = DynNets.score_driven_filter_or_dgp(model,gasParVec,indTvPar;
                              vConstPar = constParVec)#,ftot_0= startPoint )

    filtPar_T_Nsample[:,:,n] = gasFiltPar

plot(gasFiltPar)

##  Load Samles of many SD static parameters, filter the data once for each set
 # of static pars  and plot
using Utilities,AReg,StaticNets,DynNets , JLD,MLBase,StatsBase,CSV, RCall
 using PyCall; pygui(); using PyPlot
 using JLD2, GLM
 Nsample = 1
 dgpType = "sin"
 T = 50
 N = 50
 Nconf_bands = 500
 Nterms = 2
 Nsteps1 ,Nsteps2 = 0,1
 load_fold = "./data/estimatesTest/sdergmTest/conf_bands/"
 @load(load_fold*"conf_bands_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE_target_$(targetAllTv)_N_conf_$(Nconf_bands)_robust.jld" ,
     stats_T, changeStats_T,estParSS_T,sampledMat_T ,parMatDgp_T,
     Nsample,T,N,filtPar_T_Nsample,gasParEst,convFlag,pVals_Nsample,scoreAutoc_Nsample,staticEst,gasParSampled)

gasFiltPar_conf_all = fill(zeros(N_samp_par,Nterms,T),Nsample)
 tmp = Array{Array{Real,2},2}(undef,T,Nsample); for t=1:T,s=1:Nsample tmp[t,s] =  Real.(changeStats_T[t][s]);end;#changeStats_T = tmp
 # Iterate across all SD static pars vectors
 for n=1
     @show(n)
    indTvPar = BitArray([true,true])
    model = DynNets.SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,indTvPar,"")

    N_samp_par = size(gasParSampled[n])[2]
    # Works only when all parameters are time varying
    gasFiltPar_conf = zeros(N_samp_par,Nterms,T)

    for i= 1:N_samp_par
        @show(i)
        samPar = gasParSampled[n][:,i]
        gasParVec = samPar
        tmp, ~ = DynNets.score_driven_filter_or_dgp(model,gasParVec,indTvPar )
            gasFiltPar_conf[i,:,:] = tmp
    end
    gasFiltPar_conf_all[n] = gasFiltPar_conf
 end
using Statistics
# Compute confidence bands
quant_vals = [0.95, 0.05]
 conf_band = zeros(length(quant_vals),T,2)
 for t=1:T
     for k=1:Nterms
      filt_t = gasFiltPar_conf_all[n][:,k,t]
      conf_band[:,t,k] = Statistics.quantile(filt_t[.!isnan.(filt_t)],quant_vals)
    end
 end

# Plot DGP filtered Path and confidence bands
figure()
     legTex = ["DGP";"SDERGM"; "95%"]
     n=1
     indTvPar = BitArray([true,true])
     model = DynNets.SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,indTvPar,"")
     estPar = gasParEst[n] # store gas parameters
     gasFiltPar  = filtPar_T_Nsample[:,:,n]
     parInd = 1
     subplot(1,2,1);
     plot(1:T,parMatDgp_T[parInd,:],"k",linewidth=4)
     plot(1:T,gasFiltPar[parInd,:],"r")
     plt.fill_between(1:T, conf_band[1,:,parInd], y2 =conf_band[2,:,parInd],color =(0.9, 0.2 , 0.2, 0.1)  )#, color='b', alpha=.1)
     #N_bands = length(gasFiltPar_conf_all[n])
     # for i=1:N_bands
     #     plot(1:T,gasFiltPar_conf_all[n][i][parInd,:],"r",alpha = 0.02)
     # end
    parInd = 2
    subplot(1,2,2);
    plot(1:T,parMatDgp_T[parInd,:],"k",linewidth=4)
    plot(1:T,gasFiltPar[parInd,:],"r")
     plt.fill_between(1:T, conf_band[1,:,parInd], y2 =conf_band[2,:,parInd],color =(0.9, 0.2 , 0.2, 0.1)  )#, color='b', alpha=.1)
    namePar1 = "Number of Links"
    namePar2 = "GWESP"
    subplot(1,2,1);    title(namePar1); legend(legTex)
    subplot(1,2,2);  title(namePar2 ); legend(legTex)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
export iauPpp
"""
**  P-vector addition.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  vector/matrix support function.
**
**  Given:
**     a        double[3]      first p-vector
**     b        double[3]      second p-vector
**
**  Returned:
**     apb      double[3]      a + b
**
**  Note:
**     It is permissible to re-use the same array for any of the
**     arguments.
**
**  This revision:  2013 June 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPpp(double a[3], double b[3], double apb[3])
function iauPpp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})
   apb = zeros(Float64, 3)
   ccall((:iauPpp, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
         convert(Array{Float64, 1}, a),
         convert(Array{Float64, 1}, b),
         apb)

   return SVector{3}(apb)
end


text:
struct ConstantRateJump{F1,F2} <: AbstractJump
  rate::F1
  affect!::F2
end

struct VariableRateJump{R,F,I,T,T2} <: AbstractJump
  rate::R
  affect!::F
  idxs::I
  rootfind::Bool
  interp_points::Int
  save_positions::Tuple{Bool,Bool}
  abstol::T
  reltol::T2
end

"""
Salis H., Kaznessis Y.,  Accurate hybrid stochastic simulation of a system of
coupled chemical or biochemical reactions, Journal of Chemical Physics, 122 (5),
DOI:10.1063/1.1835951
"""
VariableRateJump(rate,affect!;
                   idxs = nothing,
                   rootfind=true,
                   save_positions=(true,true),
                   interp_points=10,
                   abstol=1e-12,reltol=0) = VariableRateJump(
                              rate,affect!,idxs,
                              rootfind,interp_points,
                              save_positions,abstol,reltol)

struct RegularJump{iip,R,C,MD}
    rate::R
    c::C
    numjumps::Int
    mark_dist::MD
    function RegularJump{iip}(rate,c,numjumps::Int; mark_dist = nothing) where iip
      new{iip,typeof(rate),typeof(c),typeof(mark_dist)}(rate,c,numjumps,mark_dist)
    end
end

DiffEqBase.isinplace(::RegularJump{iip,R,C,MD}) where {iip,R,C,MD} = iip

RegularJump(rate,c,numjumps::Int; kwargs...) = RegularJump{DiffEqBase.isinplace(rate,4)}(rate,c,numjumps;kwargs...)

# deprecate old call
function RegularJump(rate,c,dc::AbstractMatrix; constant_c=false, mark_dist = nothing)
  @warn("The RegularJump interface has changed to be matrix-free. See the documentation for more details.")
  function _c(du,u,p,t,counts,mark)
    c(dc,u,p,t,mark)
    mul!(du,dc,counts)
  end
  RegularJump{true}(rate,_c,size(dc,2);mark_dist=mark_dist)
end

struct MassActionJump{T,S,U} <: AbstractJump
  scaled_rates::T
  reactant_stoch::S
  net_stoch::U

  function MassActionJump{T,S,U}(rates::T, rs_in::S, ns::U, scale_rates::Bool, useiszero::Bool, nocopy::Bool) where {T <: AbstractVector, S, U}
    sr  = nocopy ? rates : copy(rates)
    rs = nocopy ? rs_in : copy(rs_in)
    for i in eachindex(rs)
      if useiszero && (length(rs[i]) == 1) && iszero(rs[i][1][1])
        rs[i] = typeof(rs[i])()
      end
    end

    if scale_rates && !isempty(sr)
      scalerates!(sr, rs)
    end
    new(sr, rs, ns)
  end
  function MassActionJump{T,S,U}(rate::T, rs_in::S, ns::U, scale_rates::Bool, useiszero::Bool, nocopy::Bool) where {T <: Number, S, U}
    rs = rs_in
    if useiszero && (length(rs) == 1) && iszero(rs[1][1])
      rs = typeof(rs)()
    end
    sr = scale_rates ? scalerate(rate, rs) : rate
    new(sr, rs, ns)
  end

end
MassActionJump(usr::T, rs::S, ns::U; scale_rates = true, useiszero = true, nocopy=false) where {T,S,U} = MassActionJump{T,S,U}(usr, rs, ns, scale_rates, useiszero, nocopy)

@inline get_num_majumps(maj::MassActionJump) = length(maj.scaled_rates)
@inline get_num_majumps(maj::Nothing) = 0

struct JumpSet{T1,T2,T3,T4} <: AbstractJump
  variable_jumps::T1
  constant_jumps::T2
  regular_jump::T3
  massaction_jump::T4
end
JumpSet(vj, cj, rj, maj::MassActionJump{S,T,U}) where {S <: Number, T, U} = JumpSet(vj, cj, rj, check_majump_type(maj))

JumpSet(jump::ConstantRateJump) = JumpSet((),(jump,),nothing,nothing)
JumpSet(jump::VariableRateJump) = JumpSet((jump,),(),nothing,nothing)
JumpSet(jump::RegularJump)      = JumpSet((),(),jump,nothing)
JumpSet(jump::MassActionJump)   = JumpSet((),(),nothing,jump)
JumpSet() = JumpSet((),(),nothing,nothing)
JumpSet(jb::Nothing) = JumpSet()

# For Varargs, use recursion to make it type-stable
JumpSet(jumps::AbstractJump...) = JumpSet(split_jumps((), (), nothing, nothing, jumps...)...)

# handle vector of mass action jumps
function JumpSet(vjs, cjs, rj, majv::Vector{T}) where {T <: MassActionJump}
  if isempty(majv)
    error("JumpSets do not accept empty mass action jump collections; use \"nothing\" instead.")
  end

  maj = setup_majump_to_merge(majv[1].scaled_rates, majv[1].reactant_stoch, majv[1].net_stoch)
  for i = 2:length(majv)
    massaction_jump_combine(maj, majv[i])
  end

  JumpSet(vjs, cjs, rj, maj)
end

@inline get_num_majumps(jset::JumpSet) = get_num_majumps(jset.massaction_jump)

@inline split_jumps(vj, cj, rj, maj) = vj, cj, rj, maj
@inline split_jumps(vj, cj, rj, maj, v::VariableRateJump, args...) = split_jumps((vj..., v), cj, rj, maj, args...)
@inline split_jumps(vj, cj, rj, maj, c::ConstantRateJump, args...) = split_jumps(vj, (cj..., c), rj, maj, args...)
@inline split_jumps(vj, cj, rj, maj, c::RegularJump, args...) = split_jumps(vj, cj, regular_jump_combine(rj,c), maj, args...)
@inline split_jumps(vj, cj, rj, maj, c::MassActionJump, args...) = split_jumps(vj, cj, rj, massaction_jump_combine(maj,c), args...)
@inline split_jumps(vj, cj, rj, maj, j::JumpSet, args...) = split_jumps((vj...,j.variable_jumps...),
                                                                        (cj..., j.constant_jumps...),
                                                                        regular_jump_combine(rj,j.regular_jump),
                                                                        massaction_jump_combine(maj,j.massaction_jump), args...)

regular_jump_combine(rj1::RegularJump,rj2::Nothing) = rj1
regular_jump_combine(rj1::Nothing,rj2::RegularJump) = rj2
regular_jump_combine(rj1::Nothing,rj2::Nothing) = rj1
regular_jump_combine(rj1::RegularJump,rj2::RegularJump) = error("Only one regular jump is allowed in a JumpSet")

# functionality to merge two mass action jumps together
check_majump_type(maj::MassActionJump) = maj
check_majump_type(maj::MassActionJump{S,T,U}) where {S <: Number, T, U} = setup_majump_to_merge(maj.scaled_rates, maj.reactant_stoch, maj.net_stoch)

# if given containers of rates and stoichiometry directly create a jump
function setup_majump_to_merge(sr::T, rs::AbstractVector{S}, ns::AbstractVector{U}) where {T <: AbstractVector, S <: AbstractArray, U <: AbstractArray}
  MassActionJump(sr, rs, ns; scale_rates=false)
end

# if just given the data for one jump (and not in a container) wrap in a vector
function setup_majump_to_merge(sr::T, rs::S, ns::U) where {T <: Number, S <: AbstractArray, U <: AbstractArray}
  MassActionJump([sr], [rs], [ns]; scale_rates=false)
end

# when given a collection of reactions to add to maj
function majump_merge!(maj::MassActionJump{U,V,W}, sr::U, rs::V, ns::W) where {U <: AbstractVector, V <: AbstractVector, W <: AbstractVector}
  append!(maj.scaled_rates, sr)
  append!(maj.reactant_stoch, rs)
  append!(maj.net_stoch, ns)
  maj
end

# when given a single jump's worth of data to add to maj
function majump_merge!(maj::MassActionJump{U,V,W}, sr::T, rs::S1, ns::S2) where {T <: Number, S1 <: AbstractArray, S2 <: AbstractArray, U <: AbstractVector{T}, V <: AbstractVector{S1}, W <: AbstractVector{S2}}
  push!(maj.scaled_rates, sr)
  push!(maj.reactant_stoch, rs)
  push!(maj.net_stoch, ns)
  maj
end

# when maj only stores a single jump's worth of data (and not in a collection)
# create a new jump with the merged data stored in vectors
function majump_merge!(maj::MassActionJump{T,S,U}, sr::T, rs::S, ns::U) where {T <: Number, S <: AbstractArray, U <: AbstractArray}
  MassActionJump([maj.scaled_rates, sr], [maj.reactant_stoch, rs], [maj.net_stoch, ns]; scale_rates=false)
end

massaction_jump_combine(maj1::MassActionJump, maj2::Nothing) = maj1
massaction_jump_combine(maj1::Nothing, maj2::MassActionJump) = maj2
massaction_jump_combine(maj1::Nothing, maj2::Nothing) = maj1
massaction_jump_combine(maj1::MassActionJump, maj2::MassActionJump) = majump_merge!(maj1, maj2.scaled_rates, maj2.reactant_stoch, maj2.net_stoch)

##### helper methods for unpacking rates and affects! from constant jumps #####
function get_jump_info_tuples(constant_jumps)
  if (constant_jumps !== nothing) && !isempty(constant_jumps)
    rates    = ((c.rate for c in constant_jumps)...,)
    affects! = ((c.affect! for c in constant_jumps)...,)
  else
    rates    = ()
    affects! = ()
  end

  rates, affects!
end

function get_jump_info_fwrappers(u, p, t, constant_jumps)
  RateWrapper   = FunctionWrappers.FunctionWrapper{typeof(t),Tuple{typeof(u), typeof(p), typeof(t)}}
  AffectWrapper = FunctionWrappers.FunctionWrapper{Nothing,Tuple{Any}}

  if (constant_jumps !== nothing) && !isempty(constant_jumps)
    rates    = [RateWrapper(c.rate) for c in constant_jumps]
    affects! = [AffectWrapper(x->(c.affect!(x);nothing)) for c in constant_jumps]
  else
    rates    = Vector{RateWrapper}()
    affects! = Vector{AffectWrapper}()
  end

  rates, affects!
end


text:
module Poly

using Reexport

include("FinPolynomials.jl")
include("SumProdPolynomials.jl")

@reexport using .FinPolynomials
@reexport using .SumProdPolynomials

end


text:
using StableRNGs
using Test
using Statistics
using BayesianLinearRegression

m = @model X begin
    β ~ Normal() |> iid(size(X,2))
    y ~ For(eachrow(X)) do x
        Normal(x' * β, 1)
    end
    return collect(y)
end

rng = StableRNG(42)
X = randn(rng, 20, 2)
truth = simulate(rng, m(X=X))
y = value(truth)
post = dynamicHMC(rng, m(X=X) |  (y=y,))
pred = predictive(m,:β)

lr = BayesianLinReg(X,y)
lr.updateNoise = false
lr.updatePrior = false
fit!(lr)

@testset "dynamicHMC" begin
    @test mean(v.β for v in post) ≈ lr.weights  rtol = 0.05
    @test_nowarn [rand(rng, pred(;X=X, p...)) for p in post];
    #particles(post)
end

m2 = @model X begin
    N = size(X,1)
    k = size(X,2)
    β ~ Normal() |> iid(k)
    yhat = X * β
    y ~ For(N) do j
            Normal(yhat[j], 1)
        end
end;

jointdist = m2(X=X)

@testset "logdensity, symlogdensity" begin
    # @test_nowarn symlogdensity(m2).evalf(3)
    @test logdensity_def(jointdist, (β = truth.trace.β, y=y)) ≈ -8.015157812948065
    # @test logdensity_def(jointdist, truth, codegen) ≈ -28.551921801470904
end


text:
include("../../fancy_toys.jl")

version = v"1.0.1"

name = "CUTENSOR_CUDA$(cuda_version.major)$(cuda_version.minor)"

products = [
    LibraryProduct(["libcutensor"], :libcutensor),
]

dependencies = [Dependency(PackageSpec(name="CUDA_jll", version=cuda_version))]

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

if @isdefined(sources_linux) && should_build_platform("x86_64-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux, script,
                   [Linux(:x86_64)], products, dependencies)
end

if @isdefined(sources_windows) && should_build_platform("x86_64-w64-mingw32")
    build_tarballs(ARGS, name, version, sources_windows, script,
                   [Windows(:x86_64)], products, dependencies)
end

if @isdefined(sources_macos) && should_build_platform("x86_64-apple-darwin14")
    build_tarballs(ARGS, name, version, sources_macos, script,
                   [MacOS(:x86_64)], products, dependencies)
end