text:
#=
@inline function anchor!(stream::BufferedStreams.BufferedInputStream, p, immobilize = true)
    stream.anchor = p
    stream.immobilized = immobilize
    return stream
end

@inline function upanchor!(stream::BufferedStreams.BufferedInputStream)
    @assert stream.anchor != 0 "upanchor! called with no anchor set"
    anchor = stream.anchor
    stream.anchor = 0
    stream.immobilized = false
    return anchor
end

function ensure_margin!(stream::BufferedStreams.BufferedInputStream)
    if stream.position * 20 > length(stream.buffer) * 19
        BufferedStreams.shiftdata!(stream)
    end
    return nothing
end

@inline function resize_and_copy!(dst::Vector{UInt8}, src::Vector{UInt8}, r::UnitRange{Int})
    return resize_and_copy!(dst, 1, src, r)
end

@inline function resize_and_copy!(dst::Vector{UInt8}, dstart::Int, src::Vector{UInt8}, r::UnitRange{Int})
    rlen = length(r)
    if length(dst) != dstart + rlen - 1
        resize!(dst, dstart + rlen - 1)
    end
    copyto!(dst, dstart, src, first(r), rlen)
    return dst
end

function generate_index_function(record_type, machine, init_code, actions; kwargs...)
    kwargs = Dict(kwargs)
    context = Automa.CodeGenContext(
        generator = get(kwargs, :generator, :goto),
        checkbounds = get(kwargs, :checkbounds, false),
        loopunroll = get(kwargs, :loopunroll, 0)
    )
    quote
        function index!(record::$(record_type))
            data = record.data
            p = 1
            p_end = p_eof = sizeof(data)
            initialize!(record)
            $(init_code)
            cs = $(machine.start_state)
            $(Automa.generate_exec_code(context, machine, actions))
            if cs != 0
                throw(ArgumentError(string("failed to index ", $(record_type), " ~>", repr(String(data[p:min(p+7,p_end)])))))
            end
            @assert isfilled(record)
            return record
        end
    end
end

function generate_readheader_function(reader_type, metainfo_type, machine, init_code, actions, finish_code=:())
    quote
        function readheader!(reader::$(reader_type))
            _readheader!(reader, reader.state)
        end

        function _readheader!(reader::$(reader_type), state::State)
            stream = state.stream
            ensure_margin!(stream)
            cs = state.cs
            linenum = state.linenum
            data = stream.buffer
            p = stream.position
            p_end = stream.available
            p_eof = -1
            finish_header = false
            record = $(metainfo_type)()

            $(init_code)

            while true
                $(Automa.generate_exec_code(Automa.CodeGenContext(generator=:table), machine, actions))

                state.cs = cs
                state.finished = cs == 0
                state.linenum = linenum
                stream.position = p

                if cs < 0
                    error("$($(reader_type)) file format error on line ", linenum)
                elseif finish_header
                    $(finish_code)
                    break
                elseif p > p_eof ≥ 0
                    error("incomplete $($(reader_type)) input on line ", linenum)
                else
                    hits_eof = BufferedStreams.fillbuffer!(stream) == 0
                    p = stream.position
                    p_end = stream.available
                    if hits_eof
                        p_eof = p_end
                    end
                end
            end
        end
    end
end

function generate_read_function(reader_type, machine, init_code, actions; kwargs...)
    kwargs = Dict(kwargs)
    context = Automa.CodeGenContext(
        generator=get(kwargs, :generator, :goto),
        checkbounds=get(kwargs, :checkbounds, false),
        loopunroll=get(kwargs, :loopunroll, 0)
    )
    quote
        function Base.read!(reader::$(reader_type), record::eltype($(reader_type)))::eltype($(reader_type))
            return _read!(reader, reader.state, record)
        end

        function _read!(reader::$(reader_type), state::State, record::eltype($(reader_type)))
            stream = state.stream
            ensure_margin!(stream)
            cs = state.cs
            linenum = state.linenum
            data = stream.buffer
            p = stream.position
            p_end = stream.available
            p_eof = -1
            found_record = false
            initialize!(record)

            $(init_code)

            if state.finished
                throw(EOFError())
            end

            while true
                $(Automa.generate_exec_code(context, machine, actions))

                state.cs = cs
                state.finished |= cs == 0
                state.linenum = linenum
                stream.position = p

                if cs < 0
                    error($(reader_type), " file format error on line ", linenum, " ~>", repr(String(data[p:min(p+7,p_end)])))
                elseif found_record
                    break
                elseif cs == 0
                    throw(EOFError())
                elseif p > p_eof ≥ 0
                    error("incomplete $($(reader_type)) input on line ", linenum)
                elseif BufferedStreams.available_bytes(stream) < 64
                    hits_eof = BufferedStreams.fillbuffer!(stream) == 0
                    p = stream.position
                    p_end = stream.available
                    if hits_eof
                        p_eof = p_end
                    end
                end
            end

            @assert isfilled(record)
            return record
        end
    end
end
=#

# Automa.jl generated readrecord! and readmetainfo! functions
# ========================================

# file   = header . body
# header = metainfo*
# body   = record*
const sam_metainfo_machine, sam_record_machine, sam_header_machine, sam_body_machine = (function ()

    isinteractive() && info("compiling SAM")

    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    alt = Automa.RegExp.alt
    opt = Automa.RegExp.opt
    any = Automa.RegExp.any

    metainfo = let
        tag = re"[A-Z][A-Z]" \ cat("CO")
        tag.actions[:enter] = [:pos1]
        tag.actions[:exit]  = [:metainfo_tag]

        dict = let
            key = re"[A-Za-z][A-Za-z0-9]"
            key.actions[:enter] = [:pos2]
            key.actions[:exit]  = [:metainfo_dict_key]
            val = re"[ -~]+"
            val.actions[:enter] = [:pos2]
            val.actions[:exit]  = [:metainfo_dict_val]
            keyval = cat(key, ':', val)

            cat(keyval, rep(cat('\t', keyval)))
        end
        dict.actions[:enter] = [:pos1]
        dict.actions[:exit]  = [:metainfo_val]

        co = cat("CO")
        co.actions[:enter] = [:pos1]
        co.actions[:exit]  = [:metainfo_tag]

        comment = re"[^\r\n]*"
        comment.actions[:enter] = [:pos1]
        comment.actions[:exit]  = [:metainfo_val]

        cat('@', alt(cat(tag, '\t', dict), cat(co, '\t', comment)))
    end
    metainfo.actions[:enter] = [:mark]
    metainfo.actions[:exit]  = [:metainfo]

    record = let
        qname = re"[!-?A-~]+"
        qname.actions[:enter] = [:pos]
        qname.actions[:exit]  = [:record_qname]

        flag = re"[0-9]+"
        flag.actions[:enter] = [:pos]
        flag.actions[:exit]  = [:record_flag]

        rname = re"\*|[!-()+-<>-~][!-~]*"
        rname.actions[:enter] = [:pos]
        rname.actions[:exit]  = [:record_rname]

        pos = re"[0-9]+"
        pos.actions[:enter] = [:pos]
        pos.actions[:exit]  = [:record_pos]

        mapq = re"[0-9]+"
        mapq.actions[:enter] = [:pos]
        mapq.actions[:exit]  = [:record_mapq]

        cigar = re"\*|([0-9]+[MIDNSHPX=])+"
        cigar.actions[:enter] = [:pos]
        cigar.actions[:exit]  = [:record_cigar]

        rnext = re"\*|=|[!-()+-<>-~][!-~]*"
        rnext.actions[:enter] = [:pos]
        rnext.actions[:exit]  = [:record_rnext]

        pnext = re"[0-9]+"
        pnext.actions[:enter] = [:pos]
        pnext.actions[:exit]  = [:record_pnext]

        tlen = re"[-+]?[0-9]+"
        tlen.actions[:enter] = [:pos]
        tlen.actions[:exit]  = [:record_tlen]

        seq = re"\*|[A-Za-z=.]+"
        seq.actions[:enter] = [:pos]
        seq.actions[:exit]  = [:record_seq]

        qual = re"[!-~]+"
        qual.actions[:enter] = [:pos]
        qual.actions[:exit]  = [:record_qual]

        field = let
            tag = re"[A-Za-z][A-Za-z0-9]"
            val = alt(
                re"A:[!-~]",
                re"i:[-+]?[0-9]+",
                re"f:[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?",
                re"Z:[ !-~]*",
                re"H:([0-9A-F][0-9A-F])*",
                re"B:[cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+")

            cat(tag, ':', val)
        end
        field.actions[:enter] = [:pos]
        field.actions[:exit]  = [:record_field]

        cat(
            qname, '\t',
            flag,  '\t',
            rname, '\t',
            pos,   '\t',
            mapq,  '\t',
            cigar, '\t',
            rnext, '\t',
            pnext, '\t',
            tlen,  '\t',
            seq,   '\t',
            qual,
            rep(cat('\t', field)))
    end
    record.actions[:enter] = [:mark]
    record.actions[:exit]  = [:record]

    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(re"\r?", lf)
    end

    header′ = rep(cat(metainfo, newline))
    header′.actions[:exit] = [:header]
    header = cat(header′, opt(any() \ cat('@')))  # look ahead

    body = rep(cat(record, newline))

    return map(Automa.compile, (metainfo, record, header, body))
end)()

# write("sam_metainfo_machine.dot", Automa.machine2dot(sam_metainfo_machine))
# run(`dot -Tsvg -o sam_metainfo_machine.svg sam_metainfo_machine.dot`)
#
# write("sam_record_machine.dot", Automa.machine2dot(sam_record_machine))
# run(`dot -Tsvg -o sam_record_machine.svg sam_record_machine.dot`)
#
# write("sam_header_machine.dot", Automa.machine2dot(sam_header_machine))
# run(`dot -Tsvg -o sam_header_machine.svg sam_header_machine.dot`)
#
# write("sam_body_machine.dot", Automa.machine2dot(sam_body_machine))
# run(`dot -Tsvg -o sam_body_machine.svg sam_body_machine.dot`)

function appendfrom!(dst, dpos, src, spos, n)
    if length(dst) < dpos + n - 1
        resize!(dst, dpos + n - 1)
    end
    copyto!(dst, dpos, src, spos, n)
    return dst
end

const sam_metainfo_actions = Dict(
    :mark => :(@mark),
    :pos1  => :(pos1 = @relpos(p)),
    :pos2  => :(pos2 = @relpos(p)),
    :countline => :(linenum += 1),
    :metainfo_tag => :(metainfo.tag = pos1:@relpos(p-1)),
    :metainfo_val => :(metainfo.val = pos1:@relpos(p-1)),
    :metainfo_dict_key => :(push!(metainfo.dictkey, pos2:@relpos(p-1))),
    :metainfo_dict_val => :(push!(metainfo.dictval, pos2:@relpos(p-1))),
    :metainfo => quote
        appendfrom!(metainfo.data, 1, data, @markpos(), p-@markpos())
        # metainfo.filled = 1:@relpos(p-1)
        # metainfo.filled = 1:p-1
        # metainfo.filled = markpos:p
        metainfo.filled = @markpos():p-1
        found_metainfo = true
        @escape
    end
)

sam_metainfo_context = Automa.CodeGenContext(
    generator = :goto,
    checkbounds = false,
    loopunroll = 0
)

sam_metainfo_initcode = quote
    pos1 = 0
    pos2 = 0
    filled = 0
    found_metainfo = false
    initialize!(metainfo)
    cs, linenum = state
end

sam_metainfo_loopcode = quote
    if cs < 0
        throw(ArgumentError("malformed metainfo at line $(linenum)"))
    end
    found_metainfo && @goto __return__
end

sam_metainfo_returncode = quote
    return cs, linenum, found_metainfo
end

Automa.Stream.generate_reader(
    :readmetainfo!,
    sam_metainfo_machine,
    arguments = (:(metainfo::MetaInfo), :(state::Tuple{Int,Int})),
    actions = sam_metainfo_actions,
    context = sam_metainfo_context,
    initcode = sam_metainfo_initcode,
    loopcode = sam_metainfo_loopcode,
    returncode = sam_metainfo_returncode
) |> eval

const sam_header_actions = merge(
    sam_metainfo_actions,
    Dict(
        :metainfo => quote
            let markpos = @markpos()
                # appendfrom!(metainfo.data, filled + 1, data, markpos, n)
                # appendfrom!(metainfo.data, 1, data, markpos, @relpos(p-1))
                appendfrom!(metainfo.data, 1, data, markpos, p-markpos)

                # metainfo.filled = markpos:@relpos(p-1)
                # metainfo.filled = markpos:p
                metainfo.filled = markpos:p-1

                @info ":metainfo" p pos markpos @relpos(markpos)

                # filled += n
            end
        end,
        :header => quote
            finish_header = true
            @escape
        end
    )
)

sam_header_context = Automa.CodeGenContext(
    generator = :goto,
    checkbounds = false,
    loopunroll = 0
)

sam_header_initcode = quote
    @info "sam_header_initcode" state

    pos = 0
    # filled = 0
    found_metainfo = false
    finish_header = false
    metainfo = MetaInfo()
    cs, linenum = state
end

sam_header_loopcode = quote
    @info "sam_header_loopcode" state p pos metainfo.filled String(copy(metainfo.data))

    if cs < 0
        throw(ArgumentError("malformed metainfo at line $(linenum)"))
    end

    @assert isfilled(metainfo) #TODO: move to loopcode?

    push!(reader.header.metainfo, metainfo) #TODO: move to loopcode?

    metainfo = MetaInfo()

    finish_header && @goto __return__
end

sam_header_returncode = quote
    return cs, linenum, finish_header
end

Automa.Stream.generate_reader(
    :readheader!,
    sam_header_machine,
    arguments = (:(reader::Reader), :(state::Tuple{Int,Int})),
    actions = sam_header_actions,
    context = sam_header_context,
    initcode = sam_header_initcode,
    loopcode = sam_header_loopcode,
    returncode = sam_header_returncode
) |> eval

#=
generate_index_function(
    MetaInfo,
    sam_metainfo_machine,
    :(pos1 = pos2 = offset = 0),
    sam_metainfo_actions
) |> eval

generate_readheader_function(
    Reader,
    MetaInfo,
    sam_header_machine,
    :(pos1 = pos2 = offset = 0),
    merge(sam_metainfo_actions, Dict(
        :metainfo => quote
            resize_and_copy!(record.data, data, upanchor!(stream):p-1)
            record.filled = (offset+1:p-1) .- offset
            @assert isfilled(record)
            push!(reader.header.metainfo, record)
            ensure_margin!(stream)
            record = MetaInfo()
        end,
        :header => :(finish_header = true; @escape),
        :countline => :(linenum += 1),
        :mark => :(anchor!(stream, p); offset = p - 1))),
    quote
        if !eof(stream)
            stream.position -= 1  # cancel look-ahead
        end
    end
) |> eval
=#

const sam_record_actions = Dict(
    :mark => :(@mark),
    :pos => :(pos = @relpos(p)),
    :countline => :(linenum += 1),
    :record_qname => :(record.qname = pos:@relpos(p-1)),
    :record_flag  => :(record.flag  = pos:@relpos(p-1)),
    :record_rname => :(record.rname = pos:@relpos(p-1)),
    :record_pos   => :(record.pos   = pos:@relpos(p-1)),
    :record_mapq  => :(record.mapq  = pos:@relpos(p-1)),
    :record_cigar => :(record.cigar = pos:@relpos(p-1)),
    :record_rnext => :(record.rnext = pos:@relpos(p-1)),
    :record_pnext => :(record.pnext = pos:@relpos(p-1)),
    :record_tlen  => :(record.tlen  = pos:@relpos(p-1)),
    :record_seq   => :(record.seq   = pos:@relpos(p-1)),
    :record_qual  => :(record.qual  = pos:@relpos(p-1)),
    :record_field => :(push!(record.fields, pos:@relpos(p-1))),
    :record       => quote
        let markpos = @markpos()
            # appendfrom!(record.data, 1, data, markpos, p-@markpos)
            appendfrom!(record.data, 1, data, markpos, p-markpos)

            # record.filled = markpos:@relpos(p-1)
            # record.filled = @relpos(markpos):@relpos(p-1)
            # record.filled = markpos:p
            record.filled = markpos:p-1

            found_record = true

            @info ":record" p pos markpos @relpos(markpos)
            @escape
        end
    end
)

sam_record_context = Automa.CodeGenContext(
    generator = :goto,
    checkbounds = false,
    loopunroll = 0
)

sam_record_initcode = quote
    @info "sam_record_initcode" state

    pos = 0
    # filled = 0
    found_record = false
    initialize!(record)
    cs, linenum = state
end

sam_record_loopcode = quote
    @info "sam_record_loopcode" state p pos record.filled String(copy(record.data))

    if cs < 0
        throw(ArgumentError("malformed SAM file at line $(linenum)"))
    end
    found_record && @goto __return__
end

sam_record_returncode = quote
    return cs, linenum, found_record
end

Automa.Stream.generate_reader(
    :readrecord!,
    sam_body_machine,
    arguments = (:(record::Record), :(state::Tuple{Int,Int})),
    actions = sam_record_actions,
    context = sam_record_context,
    initcode = sam_record_initcode,
    loopcode = sam_record_loopcode,
    returncode = sam_record_returncode
) |> eval

#=
generate_index_function(
    Record,
    sam_record_machine,
    :(pos = offset = 0),
    sam_record_actions
) |> eval

generate_read_function(
    Reader,
    sam_body_machine,
    :(pos = offset = 0),
    merge(sam_record_actions, Dict(
        :record    => quote
            resize_and_copy!(record.data, data, upanchor!(stream):p-1)
            record.filled = (offset+1:p-1) .- offset
            found_record = true
            @escape
        end,
        :countline => :(linenum += 1),
        :mark    => :(anchor!(stream, p); offset = p - 1))
    )
) |> eval
=#


text:
abstract type AbstractTrace end
abstract type AbstractLayout end

mutable struct GenericTrace{T <: AbstractDict{Symbol,Any}} <: AbstractTrace
    fields::T
end

function GenericTrace(kind::Union{AbstractString,Symbol},
                      fields=Dict{Symbol,Any}(); kwargs...)
    # use setindex! methods below to handle `_` substitution
    fields[:type] = kind
    gt = GenericTrace(fields)
    foreach(x -> setindex!(gt, x[2], x[1]), kwargs)
    if Symbol(kind) in [:contour, :contourcarpet, :heatmap, :heatmapgl]
        if !haskey(gt, :transpose)
            gt.transpose = true
        end
    end
    gt
end

function _layout_defaults()
    Dict{Symbol,Any}(
        :margin => Dict(:l => 50, :r => 50, :t => 60, :b => 50),
        :template => templates[templates.default],
    )
end

mutable struct Layout{T <: AbstractDict{Symbol,Any}} <: AbstractLayout
    fields::T
    subplots::Subplots

    function Layout{T}(fields::T; kwargs...) where T
        l = new{T}(merge(_layout_defaults(), fields), Subplots())
        foreach(x -> setindex!(l, x[2], x[1]), kwargs)
        l
    end
end

Layout(fields::T=Dict{Symbol,Any}(); kwargs...) where {T <: AbstractDict{Symbol,Any}} =
    Layout{T}(fields; kwargs...)

kind(gt::GenericTrace) = get(gt, :type, "scatter")
kind(l::Layout) = "layout"

# -------------------------------------------- #
# Specific types of trace or layout attributes #
# -------------------------------------------- #

function attr(fields::AbstractDict=Dict{Symbol,Any}(); kwargs...)
    # use setindex! methods below to handle `_` substitution
    s = PlotlyAttribute(fields)
    for (k, v) in kwargs
        s[k] = v
    end
    s
end
attr(x::PlotlyAttribute; kw...) = attr(;x..., kw...)

mutable struct PlotlyFrame{T <: AbstractDict{Symbol,Any}} <: AbstractPlotlyAttribute
    fields::T
    function PlotlyFrame{T}(fields::T) where T
        !(Symbol("name") in keys(fields)) && @warn("Frame should have a :name field for expected behavior")
        new{T}(fields)
    end
end

function frame(fields=Dict{Symbol,Any}(); kwargs...)
    for (k, v) in kwargs
        fields[k] = v
    end
    PlotlyFrame{Dict{Symbol,Any}}(fields)
end

abstract type AbstractLayoutAttribute <: AbstractPlotlyAttribute end
abstract type AbstractShape <: AbstractLayoutAttribute end

kind(::AbstractPlotlyAttribute) = "PlotlyAttribute"

# TODO: maybe loosen some day
const _Scalar = Union{DateTime,Date,Number,AbstractString,Symbol}

# ------ #
# Shapes #
# ------ #

mutable struct Shape <: AbstractLayoutAttribute
    fields::AbstractDict{Symbol}
end

function Shape(kind::AbstractString, fields=Dict{Symbol,Any}(); kwargs...)
    # use setindex! methods below to handle `_` substitution
    fields[:type] = kind
    s = Shape(fields)
    foreach(x -> setindex!(s, x[2], x[1]), kwargs)
    s
end

# helper method needed below
_rep(x, n) = take(cycle(x), n)

# line, circle, and rect share same x0, x1, y0, y1 args. Define methods for
# them here
for t in [:line, :circle, :rect]
    str_t = string(t)
    @eval $t(d::AbstractDict=Dict{Symbol,Any}(), ;kwargs...) =
        Shape($str_t, d; kwargs...)
    eval(Expr(:export, t))

    @eval function $(t)(x0::_Scalar, x1::_Scalar, y0::_Scalar, y1::_Scalar,
                        fields::AbstractDict=Dict{Symbol,Any}(); kwargs...)
        $(t)(fields; x0=x0, x1=x1, y0=y0, y1=y1, kwargs...)
    end

    @eval function $(t)(x0::Union{AbstractVector,_Scalar},
                        x1::Union{AbstractVector,_Scalar},
                        y0::Union{AbstractVector,_Scalar},
                        y1::Union{AbstractVector,_Scalar},
                        fields::AbstractDict=Dict{Symbol,Any}(); kwargs...)
        n = reduce(max, map(length, (x0, x1, y0, y1)))
        f(_x0, _x1, _y0, _y1) = $(t)(_x0, _x1, _y0, _y1, copy(fields); kwargs...)
        map(f, _rep(x0, n), _rep(x1, n), _rep(y0, n), _rep(y1, n))
    end
end

@doc "Draw a line through the points (x0, y0) and (x1, y2)" line

@doc """
Draw a circle from ((`x0`+`x1`)/2, (`y0`+`y1`)/2)) with radius
 (|(`x0`+`x1`)/2 - `x0`|, |(`y0`+`y1`)/2 -`y0`)|) """ circle

@doc """
Draw a rectangle linking (`x0`,`y0`), (`x1`,`y0`),
(`x1`,`y1`), (`x0`,`y1`), (`x0`,`y0`)""" rect

"Draw an arbitrary svg path"
path(p::AbstractString; kwargs...) = Shape("path"; path=p, kwargs...)

export path

# derived shapes

vline(x, ymin, ymax, fields::AbstractDict=Dict{Symbol,Any}(); kwargs...) =
    line(x, x, ymin, ymax, fields; kwargs...)

"""
`vline(x, fields::AbstractDict=Dict{Symbol,Any}(); kwargs...)`

Draw vertical lines at each point in `x` that span the height of the plot
"""
vline(x, fields::AbstractDict=Dict{Symbol,Any}(); kwargs...) =
    vline(x, 0, 1, fields; xref="x", yref="paper", kwargs...)

hline(y, xmin, xmax, fields::AbstractDict=Dict{Symbol,Any}(); kwargs...) =
    line(xmin, xmax, y, y, fields; kwargs...)

"""
`hline(y, fields::AbstractDict=Dict{Symbol,Any}(); kwargs...)`

Draw horizontal lines at each point in `y` that span the width of the plot
"""
hline(y, fields::AbstractDict=Dict{Symbol,Any}(); kwargs...) =
    hline(y, 0, 1, fields; xref="paper", yref="y", kwargs...)

# ---------------------------------------- #
# Implementation of getindex and setindex! #
# ---------------------------------------- #

const HasFields = Union{GenericTrace,Layout,Shape,PlotlyAttribute,PlotlyFrame}
const _LikeAssociative = Union{PlotlyAttribute,AbstractDict}

_symbol_dict(hf::HasFields) = _symbol_dict(hf.fields)

#= NOTE: Generate this list with the following code
using JSON, PlotlyJS, PlotlyBase
d = JSON.parsefile(Pkg.dir("PlotlyJS", "deps", "plotschema.json"))
d = PlotlyBase._symbol_dict(d)

nms = Set{Symbol}()
function add_to_names!(d::AbstractDict)
    map(add_to_names!, collect(keys(d)))
    map(add_to_names!, collect(values(d)))
    nothing
end
add_to_names!(s::Symbol) = push!(nms, s)
add_to_names!(x) = nothing

add_to_names!(d[:schema][:layout][:layoutAttributes])
for (_, v) in d[:schema][:traces]
    add_to_names!(v)
end

_UNDERSCORE_ATTRS = collect(
    filter(
        x-> contains(string(x), "_") && !startswith(string(x), "_"),
        nms
    )
) =#
const _UNDERSCORE_ATTRS = [:error_x, :copy_ystyle, :error_z, :plot_bgcolor,
                           :paper_bgcolor, :copy_zstyle, :error_y, :hr_name]

_isempty(x) = isempty(x)
_isempty(::Union{Nothing,Missing}) = true
_isempty(x::Bool) = x
_isempty(x::Union{Symbol,String}) = false

function Base.merge(hf::HasFields, d::Dict)
    out = deepcopy(hf)
    for (k, v) in d
        out[k] = d
    end
    out
end

function Base.merge!(hf1::HasFields, hf2::HasFields)
    for (k, v) in hf2.fields
        hf1[k] = v
    end
    hf1
end

Base.merge(d::Dict{Symbol}, hf2::HasFields) = merge(d, hf2.fields)
Base.merge!(d::Dict{Symbol}, hf2::HasFields) = merge!(d, hf2.fields)
Base.pairs(hf::HasFields) = pairs(hf.fields)
Base.keys(hf::HasFields) = keys(hf.fields)
Base.values(hf::HasFields) = values(hf.fields)

function setifempty!(hf::HasFields, key::Symbol, value)
    if _isempty(hf[key])
        hf[key] = value
    end
end

Base.haskey(hf::HasFields, k::Symbol) = haskey(hf.fields, k)

Base.merge(hf1::T, hf2::T) where {T <: HasFields} =
    merge!(deepcopy(hf1), hf2)

Base.isempty(hf::HasFields) = isempty(hf.fields)

function Base.get(hf::HasFields, k::Symbol, default)
    out = getindex(hf, k)
    (out == Dict()) ? default : out
end

Base.iterate(hf::HasFields) = iterate(hf.fields)
Base.iterate(hf::HasFields, x) = iterate(hf.fields, x)

==(hf1::T, hf2::T) where {T <: HasFields} = hf1.fields == hf2.fields

# NOTE: there is another method in dataframes_api.jl
_obtain_setindex_val(container::Any, val::Any, key::_Maybe{Symbol}=missing) = val
_obtain_setindex_val(container::Dict, val::Any, key::_Maybe{Symbol}=missing) = haskey(container, val) ? container[val] : val
_obtain_setindex_val(container::Any, func::Function, key::_Maybe{Symbol}=missing) = func(container)
_obtain_setindex_val(container::Missing, func::Function, key::_Maybe{Symbol}=missing) = func

# no container
Base.setindex!(gt::HasFields, val, key::String...) = setindex!(gt, val, missing, key...)

# methods that allow you to do `obj["first.second.third"] = val`
function Base.setindex!(gt::HasFields, val, container, key::String)
    if in(Symbol(key), _UNDERSCORE_ATTRS)
        return gt.fields[Symbol(key)] = _obtain_setindex_val(container, val)
    else
        return setindex!(gt, val, container, map(Symbol, split(key, ['.', '_']))...)
    end
end

Base.setindex!(gt::HasFields, val, container, keys::String...) =
    setindex!(gt, val, container, map(Symbol, keys)...)

# Now for deep setindex. The deepest the json schema ever goes is 4 levels deep
# so we will simply write out the setindex calls for 4 levels by hand. If the
# schema gets deeper in the future we can @generate them with @nexpr

# no container
Base.setindex!(gt::HasFields, val, key::Symbol...) = setindex!(gt, val, missing, key...)

function Base.setindex!(gt::HasFields, val, container, key::Symbol)
    # check if single key has underscores, if so split as str and call above
    # unless it is one of the special attribute names with an underscore
    if occursin("_", string(key))
        # check for title
        if key in [:xaxis_title, :yaxis_title, :zaxis_title, :title] && (typeof(val) in [String, Symbol])
            key = Symbol(key, :_text)
        end
        if !in(key, _UNDERSCORE_ATTRS)
            return setindex!(gt, val, container, string(key))
        end
    end
    gt.fields[key] = _obtain_setindex_val(container, val, key)
end

function Base.setindex!(gt::HasFields, val, container, k1::Symbol, k2::Symbol)
    d1 = get(gt.fields, k1, Dict())
    si_val = _obtain_setindex_val(container, val)
    d1[k2] = si_val
    gt.fields[k1] = d1
    si_val
end

function Base.setindex!(gt::HasFields, val, container, k1::Symbol, k2::Symbol, k3::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    si_val = _obtain_setindex_val(container, val)
    d2[k3] = si_val
    d1[k2] = d2
    gt.fields[k1] = d1
    si_val
end

function Base.setindex!(gt::HasFields, val, container, k1::Symbol, k2::Symbol,
                        k3::Symbol, k4::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    d3 = get(d2, k3, Dict())
    si_val = _obtain_setindex_val(container, val)
    d3[k4] = si_val
    d2[k3] = d3
    d1[k2] = d2
    gt.fields[k1] = d1
    si_val
end

function Base.setindex!(gt::HasFields, val, container, k1::Symbol, k2::Symbol,
    k3::Symbol, k4::Symbol, k5::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    d3 = get(d2, k3, Dict())
    d4 = get(d3, k4, Dict())
    si_val = _obtain_setindex_val(container, val)
    d4[k5] = si_val
    d3[k4] = d4
    d2[k3] = d3
    d1[k2] = d2
    gt.fields[k1] = d1
    si_val
end

#= NOTE: I need to special case instances when `val` is Associatve like so that
         I can partially update something that already exists.

Example:

hf = Layout(font_size=10)
val = Layout(font_family="Helvetica") =#

Base.setindex!(gt::HasFields, val::_LikeAssociative, key::Symbol...) = setindex!(gt, val, missing, key...)

function Base.setindex!(gt::HasFields, val::_LikeAssociative, container, key::Symbol)
    if occursin("_", string(key))

        if !in(key, _UNDERSCORE_ATTRS)
            return setindex!(gt, val, string(key))
        end
    end

    if key === :geojson
        gt.fields[key] = val
        return
    end

    for (k, v) in val
        setindex!(gt, v, container, key, k)
    end
end

function Base.setindex!(gt::HasFields, val::_LikeAssociative, container, k1::Symbol,
                        k2::Symbol)
    for (k, v) in val
        setindex!(gt, v, container, k1, k2, k)
    end
end

function Base.setindex!(gt::HasFields, val::_LikeAssociative, container, k1::Symbol,
                        k2::Symbol, k3::Symbol)
    for (k, v) in val
        setindex!(gt, v, container, k1, k2, k3, k)
    end
end

function Base.setproperty!(hf::HF, p::Symbol, val) where HF <: HasFields
    if hasfield(HF, p)
        return setfield!(hf, p, val)
    end
    setindex!(hf, val, p)
end

# now on to the simpler getindex methods. They will try to get the desired
# key, but if it doesn't exist an empty dict is returned
function Base.getindex(gt::HasFields, key::String)
    if in(Symbol(key), _UNDERSCORE_ATTRS)
        gt.fields[Symbol(key)]
    else
        getindex(gt, map(Symbol, split(key, ['.', '_']))...)
    end
end

Base.getindex(gt::HasFields, keys::String...) =
    getindex(gt, map(Symbol, keys)...)

function Base.getindex(gt::HasFields, key::Symbol)
    if occursin("_", string(key))
        if !in(key, _UNDERSCORE_ATTRS)
            return getindex(gt, string(key))
        end
    end
    get(gt.fields, key, Dict())
end

function Base.getindex(gt::HasFields, k1::Symbol, k2::Symbol)
    d1 = get(gt.fields, k1, Dict())
    get(d1, k2, Dict())
end

function Base.getindex(gt::HasFields, k1::Symbol, k2::Symbol, k3::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    get(d2, k3, Dict())
end

function Base.getindex(gt::HasFields, k1::Symbol, k2::Symbol,
                       k3::Symbol, k4::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    d3 = get(d2, k3, Dict())
    get(d3, k4, Dict())
end

function Base.getindex(gt::HasFields, k1::Symbol, k2::Symbol,
    k3::Symbol, k4::Symbol, k5::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    d3 = get(d2, k3, Dict())
    d4 = get(d3, k4, Dict())
    get(d4, k5, Dict())
end

function Base.getproperty(gt::HF, p::Symbol) where HF <: HasFields
    if hasfield(HF, p)
        return getfield(gt, p)
    end
    getindex(gt, p)
end

# Now to the pop! methods
function Base.pop!(gt::HasFields, key::String)
    if in(Symbol(key), _UNDERSCORE_ATTRS)
        pop!(gt.fields, Symbol(key))
    else
        pop!(gt, map(Symbol, split(key, ['.', '_']))...)
    end
end

Base.pop!(gt::HasFields, keys::String...) =
    pop!(gt, map(Symbol, keys)...)

function Base.pop!(gt::HasFields, key::Symbol)
    if occursin("_", string(key))
        if !in(key, _UNDERSCORE_ATTRS)
            return pop!(gt, string(key))
        end
    end
    pop!(gt.fields, key, Dict())
end

function Base.pop!(gt::HasFields, k1::Symbol, k2::Symbol)
    d1 = get(gt.fields, k1, Dict())
    pop!(d1, k2, Dict())
end

function Base.pop!(gt::HasFields, k1::Symbol, k2::Symbol, k3::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    pop!(d2, k3, Dict())
end

function Base.pop!(gt::HasFields, k1::Symbol, k2::Symbol,
                       k3::Symbol, k4::Symbol)
    d1 = get(gt.fields, k1, Dict())
    d2 = get(d1, k2, Dict())
    d3 = get(d2, k3, Dict())
    pop!(d3, k4, Dict())
end

# Function used to have meaningful display of traces and layouts
function _describe(x::HasFields)
    fields = sort(map(String, collect(keys(x.fields))))
    n_fields = length(fields)
    if n_fields == 0
        return "$(kind(x)) with no fields"
    elseif n_fields == 1
        return "$(kind(x)) with field $(fields[1])"
    elseif n_fields == 2
        return "$(kind(x)) with fields $(fields[1]) and $(fields[2])"
    else
        return "$(kind(x)) with fields $(join(fields, ", ", ", and "))"
    end
end

Base.show(io::IO, ::MIME"text/plain", g::HasFields) =
    println(io, _describe(g))

ext: jl


text:
struct StandardRidge{T} <: AbstractLinearModel
    regularization_coeff::T
end

"""
    StandardRidge(regularization_coeff)
    StandardRidge(;regularization_coeff=0.0)

Ridge regression training for all the models in the library. The
```regularization_coeff``` is the regularization, it can be passed as an arg or kwarg.
"""
function StandardRidge(;regularization_coeff=0.0)
    return StandardRidge(regularization_coeff)
end

#default training - OLS
function _train(states, target_data, sr::StandardRidge=StandardRidge(0.0))
    typeof(target_data[1])(sr.regularization_coeff)
    output_layer = ((states*states'+sr.regularization_coeff*I) \ (states*target_data'))'
    #output_layer = (target_data*states')*inv(states*states'+sr.regularization_coeff*I)
    return OutputLayer(sr, output_layer, size(target_data, 1), target_data[:,end])
end

#mlj interface
struct LinearModel{T,S,K} <: AbstractLinearModel
    regression::T
    solver::S
    regression_kwargs::K
end

"""
    LinearModel(;regression=LinearRegression,
        solver=Analytical(),
        regression_kwargs=(;))

Linear regression training based on
[MLJLinearModels](https://juliaai.github.io/MLJLinearModels.jl/stable/) for all the
models in the library. All the parameters have to be passed into ```regression_kwargs```,
apart from the solver choice. MLJLinearModels.jl needs to be called in order
to use these models.
"""
function LinearModel(
    ;regression=LinearRegression,
    solver=Analytical(),
    regression_kwargs=(;)
)
    return LinearModel(regression, solver, regression_kwargs)
end

function LinearModel(
    regression;
    solver=Analytical(),
    regression_kwargs=(;)
)
    return LinearModel(regression, solver, regression_kwargs)
end

function _train(states, target_data, linear::LinearModel)
    out_size = size(target_data, 1)
    output_layer = zeros(size(target_data, 1), size(states, 1))
    for i in 1:size(target_data, 1)
        regressor = linear.regression(; fit_intercept = false, linear.regression_kwargs...)
        output_layer[i,:] = MLJLinearModels.fit(regressor, states',
        target_data[i,:], solver = linear.solver)
    end

    return OutputLayer(linear, output_layer, out_size, target_data[:,end])
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


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

using Markdown
using InteractiveUtils

# ╔═╡ 79789782-8d97-11eb-112c-fd8c2c294527
using Pkg, DrWatson

# ╔═╡ dd81850e-8d97-11eb-0cff-8bf4f38f0660
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 3cb90c0e-8db1-11eb-324a-0f66762c31f6
md" ## Clip-11-09-14s.jl"

# ╔═╡ dd81af5c-8d97-11eb-10e4-73daebcd8875
begin
    df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame)
    df[!, "treatment"] = 1 .+ df.prosoc_left + 2 * df.condition
	[count(==(element), df.treatment) for element in unique(df.treatment) ]
end

# ╔═╡ dd82254a-8d97-11eb-3f3d-59b7642304f5
stan11_4 = "
data{
	int N;
	int K;
    int pulled_left[N];
    int treatment[N];
    int actor[N];
}
parameters{
    vector[7] a;
    vector[4] b;
}
model{
    vector[N] p;
    b ~ normal( 0 , 0.5 );
    a ~ normal( 0 , 1.5 );
    for ( i in 1:N ) {
        p[i] = a[actor[i]] + b[treatment[i]];
        p[i] = inv_logit(p[i]);
    }
    pulled_left ~ binomial( 1 , p );
}
generated quantities{
    vector[N] log_lik;
    vector[N] p;
    for ( i in 1:N ) {
        p[i] = a[actor[i]] + b[treatment[i]];
        p[i] = inv_logit(p[i]);
    }
    for ( i in 1:N ) log_lik[i] = binomial_lpmf( pulled_left[i] | 1 , p[i] );
}
";

# ╔═╡ dd8c8ce4-8d97-11eb-1249-03382bfe3feb
begin
	m11_4s = SampleModel("m11.4s", stan11_4)
	data = (N = size(df, 1), K = length(unique(df.treatment)),
		pulled_left = df.pulled_left, actor=df.actor,
		treatment = df.treatment)
	rc11_4s = stan_sample(m11_4s; data)
end;

# ╔═╡ dd96b924-8d97-11eb-0d47-ffa72b39a875
if success(rc11_4s)
	post11_4s_df = read_samples(m11_4s, :dataframe)
	PRECIS(post11_4s_df[:, 1:11])
end

# ╔═╡ 83d6ef06-8e48-11eb-3cf2-f93f08c47e5d
md" ##### rethinking result:"

# ╔═╡ 75bc9c60-8e48-11eb-2eaa-5fac58b91570
md"
```
      mean   sd  5.5% 94.5% n_eff Rhat4
a[1] -0.43 0.33 -0.95  0.12   672  1.00
a[2]  3.93 0.76  2.83  5.22  1301  1.00
a[3] -0.73 0.34 -1.29 -0.18   845  1.00
a[4] -0.73 0.34 -1.31 -0.20   702  1.00
a[5] -0.43 0.33 -0.95  0.09   656  1.00
a[6]  0.50 0.34 -0.05  1.05   665  1.00
a[7]  1.98 0.43  1.29  2.65   683  1.01
b[1] -0.06 0.29 -0.50  0.41   602  1.00
b[2]  0.46 0.30 -0.01  0.95   669  1.00
b[3] -0.41 0.30 -0.89  0.08   586  1.01
b[4]  0.35 0.29 -0.10  0.82   651  1.00
```
"

# ╔═╡ 41ccd42e-8db0-11eb-1ed5-9f897b54a7cc
begin
	nt11_4s = read_samples(m11_4s)
	p_left = [logistic.(nt11_4s.a[i]) for i in 1:7]
	mean.([p_left[i] for i in 1:7])
end

# ╔═╡ 5830c0b0-9001-11eb-1a21-b9b059a0e0a0
begin
	name = "m11.4s"
	a_pars = Symbol.(["a.$i" for i in 1:7])
	s_a, p_a = plot_logistic_coef(post11_4s_df, a_pars, name)
	p_a
end

# ╔═╡ eab67364-8fff-11eb-3dee-f55af5f7b21c
s_a[1]

# ╔═╡ 098f6a8c-9004-11eb-347d-2144dd33fd8f
begin
	b_pars = Symbol.(["b.$i" for i in 1:4])
	s_b, p_b = plot_model_coef(m11_4s, b_pars)
	p_b
end

# ╔═╡ 997c1663-9ec3-4bb9-ad2a-23b5539b8fb6
begin
	dfb = DataFrame()
	diff_b_pars = [:db13, :db24]
	dfb[!, :db13] = post11_4s_df[:, "b.1"] - post11_4s_df[:, "b.3"]
	dfb[!, :db24] = post11_4s_df[:, "b.2"] - post11_4s_df[:, "b.4"]
	s_b_diff1, p_b_diff1 = plot_model_coef(dfb, diff_b_pars;
		mname = "m11.4s")
	p_b_diff1
end

# ╔═╡ 2456f812-8db1-11eb-2e37-8f274b4509cc
md" ## End of clip-11-09-14s.jl"

# ╔═╡ Cell order:
# ╟─3cb90c0e-8db1-11eb-324a-0f66762c31f6
# ╠═79789782-8d97-11eb-112c-fd8c2c294527
# ╠═dd81850e-8d97-11eb-0cff-8bf4f38f0660
# ╠═dd81af5c-8d97-11eb-10e4-73daebcd8875
# ╠═dd82254a-8d97-11eb-3f3d-59b7642304f5
# ╠═dd8c8ce4-8d97-11eb-1249-03382bfe3feb
# ╠═dd96b924-8d97-11eb-0d47-ffa72b39a875
# ╟─83d6ef06-8e48-11eb-3cf2-f93f08c47e5d
# ╟─75bc9c60-8e48-11eb-2eaa-5fac58b91570
# ╠═41ccd42e-8db0-11eb-1ed5-9f897b54a7cc
# ╠═5830c0b0-9001-11eb-1a21-b9b059a0e0a0
# ╠═eab67364-8fff-11eb-3dee-f55af5f7b21c
# ╠═098f6a8c-9004-11eb-347d-2144dd33fd8f
# ╠═997c1663-9ec3-4bb9-ad2a-23b5539b8fb6
# ╟─2456f812-8db1-11eb-2e37-8f274b4509cc


text:
using LinearAlgebra
using Combinatorics

A = [3 1 2 2;
     1 3 8 2;
     2 3 5 1]

B = [1 0 2;
     1 3 6;
     2 2 0;
     0 1 1]

function caychy_binet(A, B)
    comb = collect(combinations(1:size(A)[2], size(A)[1]))
    res = 0
    for S in comb
        res += det(A[:,S])*det(B[S,:])
    end

    res
end

println("det(AB)=", det(A*B))
println("cauchy-binet: ", caychy_binet(A, B))


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

print("Interact Sample")

using Interact
using Plots
using Rotations

# +
function f(a)
    xs = -pi:0.01:pi
    ys = @. sin(xs - a)
    plot(xs, ys)
end

a_slider=slider(-3:0.1:3,value=0.2, label="a_slider")
p=map(f, observe(a_slider))
map(display,[a_slider, p])
# +
p=ones(3)/sqrt(3)
scatter3d(
    [p[1]],[p[2]],[p[3]],
    xlabel="x",ylabel="y",zlabel="z",
    color=:blue
)

function plot_rotated_pts(x_angle,y_angle,z_angle)
    xyz_rad = deg2rad.([x_angle,y_angle,z_angle])
    r=RotXYZ(xyz_rad...)
    q=r*p
    scatter3d!([q[1]],[q[2]],[q[3]],color=:green,legend=false)
    xlims!((-1.0,1.0)); ylims!((-1.0,1.0)); zlims!((-1.0,1.0))
end

x_axis=slider(-180:180, label="x_axis")
y_axis=slider(-180:180, label="y_axis")
z_axis=slider(-180:180, label="z_axis")

angle_plot=map(
    (x_deg,y_deg,z_deg)->plot_rotated_pts(x_deg,y_deg,z_deg),
    map(observe, [x_axis,y_axis,z_axis])...,
)
    #observe(x_axis),observe(y_axis),observe(z_axis)

map(display,[x_axis, y_axis, z_axis, angle_plot]);


text:
# Pure Julia implementation of decompression of zlib and gzip
# compressed data, as specified by RFC 1950-1952.
# Author: Gunnar Farneback
# See https://groups.google.com/forum/#!searchin/julia-users/decompress/julia-users/HbyFaJjfSs0/IhI1oqYbAtEJs
#
# Cut and paste from gist https://gist.github.com/GunnarFarneback/8254567 apart from small change to take into
# account concatenated gzips
#

# Huffman codes are internally represented by Vector{Vector{Int}},
# where code[k] are a vector of the values with code words of length
# k. Codes are assigned in order from shorter to longer codes and in
# the order listed. E.g.
# [[], [2, 7], [1, 3, 5], [6, 4]]
# would be the code
# 00    - 2
# 01    - 7
# 100   - 1
# 101   - 3
# 110   - 5
# 1110  - 6
# 1111  - 4

const fixed_literal_or_length_table = (Vector{Int})[Int[],
                                                    Int[],
                                                    Int[],
                                                    Int[],
                                                    Int[],
                                                    Int[],
				                    Int[256:279],
				                    Int[0:143, 280:287],
				                    Int[144:255]]

const fixed_distance_table = (Vector{Int})[Int[],
                                           Int[],
                                           Int[],
                                           Int[],
                                           Int[0:31]]

type InflateData
    bytes::Vector{Uint8}
    current_byte::Int
    bytepos::Int
    bitpos::Int
    literal_or_length_code::Vector{Vector{Int}}
    distance_code::Vector{Vector{Int}}
end

function InflateData(source::Vector{Uint8})
    InflateData(source, 0, 1, 0, fixed_literal_or_length_table, fixed_distance_table)
end

function getbit(in::InflateData)
    if in.bitpos == 0
        in.current_byte = int(in.bytes[in.bytepos])
    end
    b = in.current_byte & 1
    in.bitpos += 1
    if in.bitpos == 8
        in.bitpos = 0
        in.bytepos += 1
    else
        in.current_byte >>= 1
    end
    b
end

function getbits(in::InflateData, n::Int)
    b = 0
    for i = 0:(n-1)
        b |= getbit(in) << i
    end
    b
end

function skip_bits_to_byte_boundary(in::InflateData)
    if in.bitpos > 0
        in.bytepos += 1
        in.bitpos = 0
    end
    nothing
end

# It is the responsibility of the caller to make sure that bitpos is
# at zero, e.g. by calling skip_bits_to_byte_boundary.
function get_aligned_byte(in::InflateData)
    byte = in.bytes[in.bytepos]
    in.bytepos += 1
    byte
end

function get_value_from_code(in::InflateData, code::Vector{Vector{Int}})
    v = 0
    for i = 1:length(code)
        v = (v << 1) | getbit(in)
        if v < length(code[i])
            return code[i][1 + v]
        end
        v -= length(code[i])
    end
    throw("incomplete code table")
end

get_literal_or_length(in) = get_value_from_code(in, in.literal_or_length_code)

const base_length = [11, 13, 15, 17, 19, 23, 27, 31, 35, 43, 51, 59, 67, 83, 99, 115, 131, 163, 195, 227]
const extra_length_bits = [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5]

function getlength(in::InflateData, v::Int)
    if v <= 264
        return v - 254
    elseif v <= 284
        return base_length[v - 264] + getbits(in, extra_length_bits[v - 264])
    else
        return 258
    end
end

function getdist(in::InflateData)
    b = get_value_from_code(in, in.distance_code)
    if b <= 3
        return b + 1
    else
        extra_bits = fld(b - 2, 2)
        return 1 + ((2 + b % 2) << extra_bits) + getbits(in, extra_bits)
    end
end

function transform_code_lengths_to_code(code_lengths::Vector{Int})
    code = Vector{Int}[]
    for i = 1:length(code_lengths)
        n = code_lengths[i]
        if n > 0
            while n > length(code)
                push!(code, Int[])
            end
            push!(code[n], i - 1)
        end
    end
    return code
end

const order = [16, 17, 18, 0, 8, 7, 9, 6, 10, 5, 11, 4, 12, 3, 13, 2, 14, 1, 15]

function read_code_tables(in::InflateData)
    hlit = getbits(in, 5) + 257
    hdist = getbits(in, 5) + 1
    hclen = getbits(in, 4) + 4
    code_length_code_lengths = zeros(Int, 19)
    for i = 1:hclen
        code_length_code_lengths[1 + order[i]] = getbits(in, 3)
    end
    code_length_code = transform_code_lengths_to_code(code_length_code_lengths)
    code_lengths = zeros(Int, hlit + hdist)
    i = 1
    while i <= hlit + hdist
        c = get_value_from_code(in, code_length_code)
        n = 1
        l = 0
        if c < 16
            l = c
        elseif c == 16
            n = 3 + getbits(in, 2)
            l = code_lengths[i-1]
        elseif c == 17
            n = 3 + getbits(in, 3)
        elseif c == 18
            n = 11 + getbits(in, 7)
        else
            throw("invalid code length code")
        end
        code_lengths[i:(i+n-1)] = l
        i += n
    end
    in.literal_or_length_code = transform_code_lengths_to_code(code_lengths[1:hlit])
    in.distance_code = transform_code_lengths_to_code(code_lengths[(hlit+1):end])
end
# Decompress data given in deflate format, as specified by RFC 1951.
function inflate(in::InflateData)
    dest = Uint8[]
    final = false
    while !final
        final = getbits(in, 1) == 1
        compression_mode = getbits(in, 2)
        if compression_mode == 0
            skip_bits_to_byte_boundary(in)
            len = getbits(in, 16)
            nlen = getbits(in, 16)
            if len $ nlen != 0xffff
                throw("corrupted data")
            end
            for i = 1:len
                push!(dest, get_aligned_byte(in))
            end
            continue
        elseif compression_mode == 1
            in.literal_or_length_code = fixed_literal_or_length_table
            in.distance_code = fixed_distance_table
        elseif compression_mode == 2
	    read_code_tables(in)
        else
            throw("invalid block compression mode 3")
        end

        while true
            v = get_literal_or_length(in)
            if v < 256
                push!(dest, uint8(v))
            elseif v == 256
                break
            else
                length = getlength(in, v)
                distance = getdist(in)
                for i = 1:length
                    push!(dest, dest[end - distance + 1])
                end
            end
        end
    end

    return dest
end

function compute_adler_checksum(x::Vector{Uint8})
    s1 = 1
    s2 = 0
    for b in x
        s1 += b
        if s1 >= 65521
            s1 -= 65521
        end
        s2 += s1
        if s2 >= 65521
            s2 -= 65521
        end
    end
    return (s2 << 16) | s1
end

const crc_table = zeros(Uint32, 256)

function make_crc_table()
    for n = 1:256
        c = uint32(n - 1)
        for k = 1:8
            if bool(c & 0x00000001)
                c = 0xedb88320 $ (c >> 1)
            else
                c >>= 1
            end
        end
        crc_table[n] = c
    end
end

function crc(x::Vector{Uint8})
    if crc_table[1] == 0
        make_crc_table()
    end
    c = 0xffffffff
    for b in x
        c = crc_table[1 + ((c $ b) & 0xff)] $ (c >> 8)
    end
    c $ 0xffffffff
end

function read_zero_terminated_data(in)
    s = Uint8[]
    while true
        c = get_aligned_byte(in)
        push!(s, c)
        if c == 0
            break
        end
    end
end

# Decompress Zlib compressed data, as specified by RFC 1950.
function decompress(source::Vector{Uint8})
    in = InflateData(source)

    CMF = get_aligned_byte(in)
    FLG = get_aligned_byte(in)
    CM = CMF & 0x0f
    CINFO = CM >> 4
    FLEVEL = FLG >> 6
    FDICT = (FLG >> 5) & 0x01
    if CM != 8
        throw("unsupported compression method")
    end
    if CINFO > 7
        throw("invalid LZ77 window size")
    end
    if FDICT != 0
        throw("preset dictionary not supported")
    end
    if mod((uint(CMF) << 8) | FLG, 31) != 0
        throw("header checksum error")
    end

    dest = inflate(in)

    skip_bits_to_byte_boundary(in)
    ADLER = 0
    for i = [24, 16, 8, 0]
        ADLER |= int(get_aligned_byte(in)) << i
    end
    if compute_adler_checksum(dest) != ADLER
        throw("corrupted data")
    end

    dest
end

# Decompress GZIP compressed data, as specified by RFC 1952.

function check_gunzip_block_header( in::InflateData )
    ID1 = get_aligned_byte(in)
    ID2 = get_aligned_byte(in)
    if ID1 != 0x1f || ID2 != 0x8b
        throw("not gzipped data")
    end
    CM = get_aligned_byte(in)
    if CM != 8
        throw("unsupported compression method")
    end
    FLG = get_aligned_byte(in)
    MTIME = getbits(in, 32)
    XFL = get_aligned_byte(in)
    OS = get_aligned_byte(in)

    if bool(FLG & 0x04)
        # FLG.FEXTRA
        xlen = getbits(in, 16)
        in.bytepos += xlen
    end

    if bool(FLG & 0x08)
        # FLG.FNAME
        name = read_zero_terminated_data(in)
    end

    if bool(FLG & 0x10)
        # FLG.FCOMMENT
        comment = read_zero_terminated_data(in)
    end

    if bool(FLG & 0x02)
        # FLG.FHCRC
        crc16 = getbits(in, 16)
        if check_crc && crc16 != (crc(in.bytes[1:(in.bytepos-3)]) & 0xffff)
            throw("corrupted data, header crc check failed")
        end
    end
end

function check_gunzip_block_tail(in::InflateData,dest::Vector{Uint8}; check_crc::Bool = true)
    skip_bits_to_byte_boundary(in)
    crc32 = getbits(in, 32)
    if check_crc && crc32 != crc(dest)
        throw("corrupted data, crc check failed")
    end
    isize = getbits(in, 32)
    if isize != length(dest)
        throw("corrupted data, length check failed")
    end
end

function gunzip( source::Vector{Uint8}; check_crc::Bool = true)
    in = InflateData(source)
    destall = Uint8[]
    while in.bytepos < length(source)
        check_gunzip_block_header(in)
        dest=inflate(in)
        append!(destall,dest)
        check_gunzip_block_tail(in,dest,check_crc=check_crc)
    end
    return destall
end

function gunzip_old(source::Vector{Uint8}; check_crc::Bool = true)
    in = InflateData(source)

    ID1 = get_aligned_byte(in)
    ID2 = get_aligned_byte(in)
    if ID1 != 0x1f || ID2 != 0x8b
        throw("not gzipped data")
    end
    CM = get_aligned_byte(in)
    if CM != 8
        throw("unsupported compression method")
    end
    FLG = get_aligned_byte(in)
    MTIME = getbits(in, 32)
    XFL = get_aligned_byte(in)
    OS = get_aligned_byte(in)

    if bool(FLG & 0x04)
        # FLG.FEXTRA
        xlen = getbits(in, 16)
        in.bytepos += xlen
    end

    if bool(FLG & 0x08)
        # FLG.FNAME
        name = read_zero_terminated_data(in)
    end

    if bool(FLG & 0x10)
        # FLG.FCOMMENT
        comment = read_zero_terminated_data(in)
    end

    if bool(FLG & 0x02)
        # FLG.FHCRC
        crc16 = getbits(in, 16)
        if check_crc && crc16 != (crc(in.bytes[1:(in.bytepos-3)]) & 0xffff)
            throw("corrupted data, header crc check failed")
        end
    end

    dest = inflate(in)

    skip_bits_to_byte_boundary(in)
    crc32 = getbits(in, 32)
    if check_crc && crc32 != crc(dest)
        throw("corrupted data, crc check failed")
    end
    isize = getbits(in, 32)
    if isize != length(dest)
        throw("corrupted data, length check failed")
    end

    dest
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# * https://discourse.julialang.org/t/for-monte-carlo-simulation-with-same-code-same-algorithm-how-fast-is-julia-compared-with-fortran/67021/41
# * https://discourse.julialang.org/t/for-monte-carlo-simulation-with-same-code-same-algorithm-how-fast-is-julia-compared-with-fortran/67021/43
# * https://github.com/genkuroki/public/blob/main/0018/%40defkwargs.ipynb
# * https://github.com/genkuroki/public/blob/main/0018/%40defunpack.ipynb
# * https://github.com/genkuroki/public/blob/main/0018/How%20to%20define%20unpacking%20macros.ipynb
# * https://github.com/genkuroki/public/blob/main/0018/How%20to%20define%20unpacking%20macros%20Part%202.ipynb

# %%
"""
The `DefUnPack` module only exports the `@defunpack` macro,
which defines a macro unpacking properties of an object.

__Simple Intended Usage:__
```
julia> using .DefUnPack

julia> struct Foo{A, B, C} a::A; b::B; c::C end

julia> @defunpack all_Foo Foo
@unpackall_Foo (macro with 1 method)

help?> @unpackall_Foo
@unpackall_Foo(obj) unpacks the fields (:a, :b, :c) of obj.

julia> @macroexpand @unpackall_Foo foo
:((a, b, c) = (foo.a, foo.b, foo.c))

julia> @unpackall_Foo Foo(1, 2.0, "three")
(1, 2.0, "three")

julia> a, b, c
(1, 2.0, "three")
```
"""
module DefUnPack

export @defunpack

"""
    @defunpack(name::Symbol, expr)

defines the macro named `Symbol(:unpack, name)`
which unpacks the properties specified by `expr` of an object.

Let `val` be the value of `expr`. Then the list of the unpacking properties is set to

* `val` if `val` is a tuple of symbols,
* `fieldnames(val)` if `val` is a type,
* `propertynames(val)` otherwise.

__Example:__

```
julia> @defunpack _cat_and_dog (:cat, :dog)
@unpack_cat_and_dog (macro with 1 method)

help?> @unpack_cat_and_dog
@unpack_cat_and_dog(obj) unpacks the fields (:cat, :dog) of obj.

julia> @macroexpand @unpack_cat_and_dog x
:((cat, dog) = (x.cat, x.dog))

julia> @unpack_cat_and_dog (dog = "bowwow", mouse = "squeak", cat = "meow")
("meow", "bowwow")

julia> cat, dog
("meow", "bowwow")

julia> struct Foo{A, B, C} a::A; b::B; c::C end

julia> @defunpack all_Foo Foo
@unpackall_Foo (macro with 1 method)

help?> @unpackall_Foo
@unpackall_Foo(obj) unpacks the fields (:a, :b, :c) of obj.

julia> @macroexpand @unpackall_Foo foo
:((a, b, c) = (foo.a, foo.b, foo.c))

julia> @unpackall_Foo Foo(1, 2.0, "three")
(1, 2.0, "three")

julia> a, b, c
(1, 2.0, "three")

julia> baz = (p = "one", q = 2.0, r = 3)

julia> @defunpack all_baz baz
@unpackall_baz (macro with 1 method)

help?> @unpackall_baz
@unpackall_baz(obj) unpacks the fields (:p, :q, :r) of obj.

julia> @macroexpand @unpackall_baz baz
:((p, q, r) = (baz.p, baz.q, baz.r))

julia> @unpackall_baz baz
("one", '2', 3)

julia> p, q, r
("one", '2', 3)

```
"""
macro defunpack(name::Symbol, expr)
    macroname = Symbol(:unpack, name)
    atmacroname = Symbol('@', macroname)
    val = Core.eval(__module__, expr)
    names = val isa Tuple{Vararg{Symbol}} ? val :
            val isa Type ? fieldnames(val) : propertynames(val)
    docstr = """`$atmacroname(obj)` unpacks the properties `$names` of `obj`."""
    quote
        macro $macroname(obj)
            Expr(:(=),
                Expr(:tuple, $names...),
                Expr(:tuple, (:($obj.$name) for name in $names)...)
            ) |> esc
        end
        @doc $docstr $(:($atmacroname))
        $atmacroname
    end |> esc
end

end

# %%
@doc DefUnPack

# %%
using .DefUnPack
@doc @defunpack

# %%
@defunpack _cat_and_dog (:cat, :dog)

# %%
?@unpack_cat_and_dog

# %%
@macroexpand @unpack_cat_and_dog x

# %%
@unpack_cat_and_dog (dog = "bowwow", mouse = "squeak", cat = "meow")

# %%
cat, dog

# %% tags=[]
struct Foo{A, B, C} a::A; b::B; c::C end

# %%
@defunpack all_Foo Foo

# %%
?@unpackall_Foo

# %%
@macroexpand @unpackall_Foo foo

# %%
@unpackall_Foo Foo(1, 2.0, "three")

# %%
a, b, c

# %%
function f(foo::Foo)
    @unpackall_Foo foo
    @show a b c
    a, b, c
end

# %%
f(Foo(1, 2.0, "three"))

# %%
baz = (p = "one", q = 2.0, r = 3)

# %%
@defunpack all_baz baz

# %%
?@unpackall_baz

# %%
@macroexpand @unpackall_baz baz

# %%
@unpackall_baz baz

# %%
p, q, r

# %%
module A
using ..DefUnPack
struct Bar{X, Y} x::X; y::Y end
@defunpack Bar Bar
end

# %%
?A.@unpackBar

# %%
@macroexpand A.@unpackBar bar

# %%
A.@unpackBar A.Bar(1, 2)

# %%
x, y

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# ### parametric function plot
# http://docs.juliaplots.org/latest/examples/pyplot/

#md # ![FILE_NAME.png](images/FILE_NAME.png)

xₜ(t) = sin(t)
yₜ(t) = sin(2t)

plot(xₜ, yₜ, 0, 2π, leg=false, fill=(0,:orange))


text:
using MinFEM

mesh = import_mesh("../meshes/stokes.msh");
L = asmVectorLaplacian(mesh, qdim=2);
B = asmIncompressibility(mesh);
C = MinFEM.asmStokesStabilization(mesh);
alpha = 0.05;

A = [L B; B' -alpha*C];

boundingBox = MinFEM.getMeshBoundingBox(mesh)
f(x) = [0.5*cos(2*pi*x[2]/(boundingBox[1][2]-boundingBox[2][2]))+0.5,0.0]
bc = evaluateMeshFunction(mesh, f, region=mesh.Boundaries[1003].Nodes,
                                   qdim=2)

rhs = zeros(3*mesh.nnodes)
asmDirichletCondition(A, union(mesh.Boundaries[1001].Nodes,
                               mesh.Boundaries[1002].Nodes,
                               mesh.Boundaries[1003].Nodes),
                      rhs=rhs, bc=bc, qdim=2)

x = A\rhs;

vtkfile = open_vtk_file(mesh, "stokes.vtu")
write_point_data(vtkfile, [reshape(x[1:2*mesh.nnodes], 2, mesh.nnodes); zeros(1,mesh.nnodes)], "v")
write_point_data(vtkfile, x[2*mesh.nnodes+1:end], "p")
save_vtk_file(vtkfile)


text:
"""
Generate figure 3 which shows embedding of the lensed grid from Section 6.
===================================================== """

"""
Load packages. Set prelim grids etc.
----------------------------------------------------- """
srcpath  = "/Users/ethananderes/Dropbox/BayesLense/src/"
savepath = "/Users/ethananderes/Dropbox/BayesLense/paper/"

include(srcpath*"Interp.jl")
include(srcpath*"cmb.jl")
include(srcpath*"fft.jl")

using PyPlot
using Interp

seed = 12
srand(seed)

par = setpar(0.35, 2^11, 0, 0, 3000, 1000, srcpath)
x,y = meshgrid(1:7,1:7)

"""
Generate figure3a.pdf
----------------------------------------------------- """
phik   = fft2(randn(size(par.grd.x))./par.grd.deltx, par) .* sqrt(par.cPP)
phidx1 = ifft2r(im .* par.grd.k1 .* phik, par)
phidx2 = ifft2r(im .* par.grd.k2 .* phik, par)

grng = linspace(1, 100, size(x,1)) |> int
fac = 1.9e3
xnew = x + fac*phidx1[grng,grng] .- mean(fac*phidx1[grng,grng])
ynew = y + fac*phidx2[grng,grng] .- mean(fac*phidx2[grng,grng])
plt.scatter(xnew, ynew, s = 58, facecolors="k")
plt.annotate(L"$(x+\nabla\phi(x),\, data(x))$",
            xy=(xnew[6,6], ynew[6,6]),
            xytext=(4.5, 7.2),
            fontsize=18,
            arrowprops={:arrowstyle=>"->"}
)
axis("tight")
axis("off")
savefig(joinpath(savepath, "figure3a.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)

"""
Generate figure3b.pdf
----------------------------------------------------- """
a = linspace(2,8, 25) - ones(1,25) - .2
at = a' +.1
xnewround = zeros(size(xnew))
ynewround = zeros(size(ynew))
for k in 1:length(xnewround)
    indd = indmin( (xnew[k] .- a).^2 + (ynew[k] .- at).^2)
    xnewround[k] = a[indd]
    ynewround[k] = at[indd]
end
plt.scatter(xnewround, ynewround, s = 58, facecolors="k")
plt.scatter(a, at, facecolors="none", edgecolors="k", s = 58, lw = .95) #, alpha = .95)
plt.annotate(L"$(y,\, T(y) + \tilde n(y))$",
  xy=(xnewround[6,6], ynewround[6,6]),
  xytext= (4.5, 7.2),
  fontsize = 18,
  arrowprops={:arrowstyle => "->"}
)
plt.annotate(L"$(y,\, T(y) + \tilde n(y))$",
  xy=(a[22,23], at[22,24]),
  xytext = (4.5, 7.2),
  fontsize = 18,
  arrowprops={:arrowstyle => "->"}
)
axis("tight")
axis("off")
savefig(joinpath(savepath, "figure3b.pdf"), dpi=300, bbox_inches="tight")

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
## Types
@concrete struct Saturated <: Function
    f
    lb
    ub
end
(f::Saturated)(args...; kwargs...) = clamp(f.f(args...; kwargs...), f.lb, f.ub)

# Proportional-derivative controller
function PD_controller(;kp, kd)
    gains = ComponentArray(;kp, kd)
    f = (x, p, t) -> -p.kp*x.θ - p.kd*x.ω
    return gains, f
end

# Single-layer dense neural controller
function dense_neural_controller(layer_size, act_fun=identity, out_fun=only)
    layers = ComponentArray(;
        L1 = dense_layer(2, layer_size),
        L2 = dense_layer(layer_size, 1),
    )

    f = function (x, p, t)
        @unpack L1, L2 = p
        out = act_fun.(L2.W * act_fun.(L1.W * x + L1.b) + L2.b)
        return out_fun(out)
    end

    return layers, f
end


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# define Markov model's parameters (λ)
Ns = 3          # number of states
a = [0.1; 0.3; 0.6]  # initial state probability
A = [0.1 0.9 0.0; 0.0 0.2 0.8; 0.5 0.0 0.5] # transition matrix
B = [0.8 0.1 0.1; 0.5 0.2 0.3; 0.7 0.1 0.2] # observation prob matrix
# Must be row stotastic
@test HMMGradients.is_transition_matrix(A)

Nt = 5                     # time window length
x = rand(1:Ns,5)           # observations

using Combinatorics
s_all = collect(multiset_permutations(1:Ns,Nt.*ones(Int,Ns),Nt))
@test length(unique(s_all)) == Ns^Nt
# this has all possible state sequences e.g. [1;1;1;1;1], [1;1;1;1;2] ...
# these are Ns^Nt permutations

# Probability of sequence of states
Ps = ones(Ns^Nt)
for z = 1:length(s_all)
  for t = 1:Nt
    if t == 1
      Ps[z] *= a[ s_all[z][t] ] # initial prob
    else
      Ps[z] *= A[ s_all[z][t-1], s_all[z][t] ] # trans prob
    end
  end
end
@test sum(Ps) ≈ 1 # check it's a probability

# Likelihood of x for a given sequence
Px_s = ones(length(s_all))
for z = 1:length(s_all)
  for t = 1:Nt
    Px_s[z] *= B[ x[t], s_all[z][t] ]
  end
end

P = sum(Ps .* Px_s)
y = [B[x[t],i] for t = 1:Nt, i = 1:Ns]

# testing Baum forward backward
alpha_b = Baum_forward(a,A,y)
alpha_m = Baum_forward_matrix(a,A,y)
beta_b  = Baum_backward(A,y)
beta_m = Baum_backward_matrix(A,y)
P2 = sum(alpha_m[end,:])
gamma_b = alpha_b .* beta_b ./ P2
@test norm(alpha_m-alpha_b) < 1e-7
@test norm(beta_m-beta_b) < 1e-7
@test P ≈ P2
@test all(sum(alpha_b.*beta_b,dims=2) .≈ P)

## testing forward backward with normalization
alpha, c = forward(Nt,a,A,y)
logalpha, logML = logforward(Nt,log.(a),log.(A),log.(y))
logbeta = logbackward(Nt,log.(A),log.(y))
beta = backward(Nt,A,c,y)
gamma = alpha.*beta
gamma2 = posterior(Nt,a,A,y)
loggamma = logalpha .+ logbeta .- logML*Nt
loggamma2 = logposterior(Nt,log.(a),log.(A),log.(y))
@test P ≈ sum(exp.(logalpha[end,:]))
@test nlogML(Nt,a,A,y) ≈ -1/Nt * log(P)
@test nlogML(Nt,a,A,y) ≈ nlogMLlog(Nt,log.(a),log.(A),log.(y))
@test all(sum(exp.(logalpha .+ logbeta),dims=2) .≈ P)
@test sum(log.(c)) ≈ log(sum(alpha_b[end,:]) )
@test all( sum(gamma,dims=2) .≈ 1.0 )
@test norm( gamma - gamma_b ) < 1e-7
@test norm( gamma2 - gamma_b ) < 1e-7
@test norm( gamma - exp.(loggamma) ) < 1e-7
@test norm( gamma - exp.(loggamma2) ) < 1e-7
# equivalence between baum alpha/beta and normalized ones
for t = 1:Nt
  @test norm( prod(1 ./ c[1:t]).*alpha_b[t,:] - alpha[t,:] ) < 1e-8
  @test norm( exp.(logalpha[t,:]) - (prod(c[1:t])*alpha[t,:] ) ) < 1e-8
end
for t = 1:Nt
  @test norm( prod(1 ./c[t+1:end]).*beta_b[t,:] - ( beta[t,:] ) )  < 1e-8
  @test norm( exp.(logbeta[t,:]) - (prod(c[t+1:end])*beta[t,:] ) )  < 1e-8
end

# testing with long sequence
T = Float64
Nt = 4000
Ns = 100
A = rand(T,Ns,Ns)
a = rand(T,Ns)
a ./= norm(a,1)
for i in 1:Ns A[i,:] ./= norm(A[i,:],1) end
@assert HMMGradients.is_transition_matrix(A)
y = rand(T,Nt,Ns)

alpha, c = forward(Nt,a,A,y)
logalpha, logML = logforward(Nt,log.(a),log.(A),log.(y))
logbeta = logbackward(Nt,log.(A),log.(y))
beta = backward(Nt,A,c,y)
gamma = alpha.*beta
loggamma = logalpha .+ logbeta .- logML*Nt
@test norm( gamma - exp.(loggamma) ) < 1e-5
@test nlogML(Nt,a,A,y) ≈ nlogMLlog(Nt,log.(a),log.(A),log.(y))


text:
include(joinpath(DFC.DEPS_DIR, "wannier90flags.jl"))
const WAN_FLAGS = _WAN_FLAGS()
issoc(c::Calculation{Wannier90}) = get(c, :spinors, false)
flagtype(::Type{Wannier90}, flag) = haskey(WAN_FLAGS, flag) ? WAN_FLAGS[flag] : Nothing
flagtype(::Calculation{Wannier90}, flag) = flagtype(Wannier90, flag)

"""
    gencalc_wan(nscf::Calculation{QE}, structure::Structure, Emin, wanflags...;
                Epad     = 5.0,
                wanexec  = Exec(exec="wannier90.x", dir=""))

Generates a Wannier90 calculation to follow on the supplied `nscf` calculation. It uses the projections defined in the `structure`, and starts counting the required amount of bands from `Emin`.
The `nscf` needs to have a valid output since it will be used in conjunction with `Emin` to find the required amount of bands and energy window for the Wannier90 calculation.
"""
function gencalc_wan(nscf::Calculation{QE}, structure::Structure, bands, Emin, wanflags...; Epad = 5.0,
                     wanexec = Exec(; exec = "wannier90.x", dir = ""))
    projs = vcat(map(structure.atoms) do x
                     ps = x.projections
                     @assert !isempty(ps) "Please first set projections for all atoms in the Structure."
                     return ps
                 end...)

    Structures.sanitize!(projs, Calculations.issoc(nscf))

    if Structures.iscolin(structure)
        wannames = ["wanup", "wandn"]
        @info "Spin polarized calculation found (inferred from nscf calculation)."
    else
        wannames = ["wan"]
    end

    wanflags = wanflags !== nothing ? Dict{Symbol,Any}(wanflags) : Dict{Symbol,Any}()
    if issoc(nscf)
        wanflags[:spinors] = true
    end

    nwann = issoc(nscf) ? 2*sum(length, projs) : sum(length, projs)
    @info "num_wann=$nwann (inferred from provided projections)."
    wanflags[:num_wann] = nwann
    kpoints = data(nscf, :k_points).data
    wanflags[:mp_grid] = length.(unique.([[n[i] for n in kpoints] for i in 1:3]))
    @info "mp_grid=$(join(wanflags[:mp_grid]," ")) (inferred from nscf calculation)."
    wanflags[:preprocess] = true
    Structures.isnoncolin(structure) && (wanflags[:spinors] = true)
    kdata = InputData(:kpoints, :none, [k[1:3] for k in kpoints])

    wancalculations = Calculation{Wannier90}[]
    for wanfil in wannames
        push!(wancalculations,
              Calculation{Wannier90}(; name = wanfil,
                                     flags = copy(wanflags), data = [kdata],
                                     exec  = wanexec, run = true))
    end

    if length(wancalculations) > 1
        set_flags!(wancalculations[1], :spin => "up")
        set_flags!(wancalculations[2], :spin => "down")
    end

    map(x -> set_wanenergies!(x, structure, bands, Emin; Epad = Epad), wancalculations)
    return wancalculations
end

## Wannier90
function kgrid(na, nb, nc, ::Type{Wannier90})
    return reshape([(a, b, c)
                    for a in collect(range(0; stop = 1, length = na + 1))[1:end-1],
                        b in collect(range(0; stop = 1, length = nb + 1))[1:end-1],
                        c in collect(range(0; stop = 1, length = nc + 1))[1:end-1]],
                   (na * nb * nc))
end

function set_kpoints!(calculation::Calculation{Wannier90}, k_grid::NTuple{3,Int};
                      print = true)
    set_flags!(calculation, :mp_grid => [k_grid...]; print = print)
    d = data(calculation, :kpoints)
    if d !== nothing
        d.data = kgrid(k_grid..., calculation)
    else
        push!(calculation.data, InputData(:kpoints, :none, kgrid(k_grid..., calculation)))
    end
    return calculation
end

"""
    set_wanenergies!(wancalculation::Calculation{Wannier90}, structure::Structure, bands::Vector{Band}, Emin::Real; Epad=5.0)

Automatically calculates and sets the wannier energies. This uses the projections,
`Emin` and the output of the nscf calculation to infer the other limits.
`Epad` allows one to specify the padding around the inner and outer energy windows
"""
function set_wanenergies!(calculation::Calculation{Wannier90}, structure::Structure, bands,
                          Emin::Real; Epad = 5.0)
    nprojs = sum(x -> sum(length.(x.projections)), structure.atoms)
    nwann = Calculations.issoc(calculation) ? 2*nprojs : nprojs
    @assert nwann != 0 "Please specify projections first."
    @info "num_wann=$nwann (inferred from provided projections)."

    if length(bands) == 2
        num_bands = length(bands[1])
        if hasflag(calculation, :spin) && calculation[:spin] == "up"
            winmin, frozmin, frozmax, winmax = wanenergyranges(Emin, nwann, bands.up, Epad)
        elseif hasflag(calculation, :spin) && calculation[:spin] == "down"
            winmin, frozmin, frozmax, winmax = wanenergyranges(Emin, nwann, bands.down,
                                                               Epad)
        end
    else
        num_bands = length(bands)
        winmin, frozmin, frozmax, winmax = wanenergyranges(Emin, nwann, bands, Epad)
    end

    set_flags!(calculation, :dis_win_min => winmin, :dis_froz_min => frozmin,
               :dis_froz_max => frozmax, :dis_win_max => winmax, :num_wann => nwann,
               :num_bands => num_bands; print = false)
    return calculation
end

#TODO: there is probably a slightly more optimal way for the frozen window, by checking max and min
#      such that every k-point has nbands inbetween.
function Emax(Emin, nbnd, bands)
    nbndfound = 0
    max = 0
    for b in bands
        if maximum(b.eigvals) >= Emin && nbndfound <= nbnd
            nbndfound += 1
            #maximum of allowed frozen window is the minimum of the first band>nbnd
            max = minimum(b.eigvals) - 0.005
        end
    end

    nbndfound < nbnd &&
        error("Number of needed bands for the projections ($nbnd) exceeds the amount of bands starting from \nEmin=$Emin ($nbndfound).\nRerun nscf with nbnd=$(length(bands) + nbnd - nbndfound).")
    return max
end

function wanenergyranges(Emin, nbnd, bands, Epad = 5)
    max = Emax(Emin, nbnd, bands)
    return (Emin - Epad, Emin, max, max + Epad)
end

function Emin_from_projwfc(structure::Structure, states, bands::Vector{DFC.Band},
                           threshold::Number)
    mask = zeros(length(states))
    for (atid, at) in enumerate(structure.atoms)
        projs = at.projections

        if isempty(projs)
            continue
        end

        stateids = Int[]
        for proj in projs
            orb = proj.orbital
            push!.((stateids,), findall(x -> x.atom_id == atid && x.l == orb.l, states))
        end
        mask[stateids] .= 1.0
    end
    Emin = -10000.0
    for b in bands
        ψ = sum(b.extra[:ψ]) / length(b.extra[:ψ])
        tot_relevant_occupation = dot(mask, ψ)

        if tot_relevant_occupation > threshold
            Emin = minimum(b.eigvals)
            break
        end
    end
    if Emin == -10000.0
        error("Couldn't find any band with occupation of relevant projections above $threshold, were any set in the structure?")
    end
    return Emin
end

# for f in (:cp, :mv)
#     @eval function Base.$f(i::Calculation{Wannier90}, dest::String; kwargs...)
#         for glob in ("$(i.name)", "UNK") # seedname should also cover generated pw2wannier90 files
#             for file in searchdir(i, glob)
#                 $f(file, joinpath(dest, splitdir(f)[end]); kwargs...)
#             end
#         end
#     end
# end


text:
# Autogenerated wrapper script for xrootdgo_jll for powerpc64le-linux-gnu
export xrootdgo

JLLWrappers.@generate_wrapper_header("xrootdgo")
JLLWrappers.@declare_library_product(xrootdgo, "xrootdgo.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        xrootdgo,
        "lib/xrootdgo.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testsuite "constructors" AT->begin
    @testset "constructors + similar" begin
        for T in supported_eltypes()
            B = AT{T}(undef, 10)
            @test B isa AT{T,1}
            @test size(B) == (10,)
            @test eltype(B) == T

            B = AT{T}(undef, 10, 10)
            @test B isa AT{T,2}
            @test size(B) == (10, 10)
            @test eltype(B) == T

            B = AT{T}(undef, (10, 10))
            @test B isa AT{T,2}
            @test size(B) == (10, 10)
            @test eltype(B) == T

            B = similar(B, Int32, 11, 15)
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(B, T)
            @test B isa AT{T,2}
            @test size(B) == (11, 15)
            @test eltype(B) == T

            B = similar(B, (5,))
            @test B isa AT{T,1}
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(B, 7)
            @test B isa AT{T,1}
            @test size(B) == (7,)
            @test eltype(B) == T

            B = similar(AT{Int32}, (11, 15))
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(AT{T}, (5,))
            @test B isa AT{T,1}
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(AT{T}, 7)
            @test B isa AT{T,1}
            @test size(B) == (7,)
            @test eltype(B) == T

            B = similar(Broadcast.Broadcasted(*, (B, B)), Int32, 11, 15)
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(Broadcast.Broadcasted(*, (B, B)), T)
            @test B isa AT{T,2}
            @test size(B) == (11, 15)
            @test eltype(B) == T
        end
    end
    @testset "comparison against Array" begin
        for typs in [(), (Int,), (Int,1), (Int,2), (Float32,), (Float32,1), (Float32,2)],
            args in [(), (1,), (1,2), ((1,),), ((1,2),),
                     (undef,), (undef, 1,), (undef, 1,2), (undef, (1,),), (undef, (1,2),),
                     (Int,), (Int, 1,), (Int, 1,2), (Int, (1,),), (Int, (1,2),),
                     ([1,2],), ([1 2],)]
            cpu = try
                Array{typs...}(args...)
            catch ex
                isa(ex, MethodError) || rethrow()
                nothing
            end

            gpu = try
                AT{typs...}(args...)
            catch ex
                isa(ex, MethodError) || rethrow()
                cpu == nothing || rethrow()
                nothing
            end

            if cpu == nothing
                @test gpu == nothing
            else
                @test typeof(cpu) == typeof(convert(Array, gpu))
            end
        end
    end
end

@testsuite "conversions" AT->begin
    for T in supported_eltypes()
        Bc = round.(rand(10, 10) .* 10.0)
        B = AT{T}(Bc)
        @test size(B) == (10, 10)
        @test eltype(B) == T
        @test Array(B) ≈ Bc

        Bc = rand(T, 10)
        B = AT(Bc)
        @test size(B) == (10,)
        @test eltype(B) == T
        @test Array(B) ≈ Bc

        Bc = rand(T, 10, 10)
        B = AT{T, 2}(Bc)
        @test size(B) == (10, 10)
        @test eltype(B) == T
        @test Array(B) ≈ Bc

        Bc = rand(Int32, 3, 3, 3)
        B = convert(AT{T, 3}, Bc)
        @test size(B) == (3, 3, 3)
        @test eltype(B) == T
        @test Array(B) ≈ Bc
    end
end

@testsuite "value constructors" AT->begin
    for T in supported_eltypes()
        x = fill(zero(T), (2, 2))
        x1 = fill(AT{T}, T(0), (2, 2))
        x2 = fill(AT{T}, T(0), (2, 2))
        x3 = fill(AT{T, 2}, T(0), (2, 2))
        @test Array(x1) ≈ x
        @test Array(x2) ≈ x
        @test Array(x3) ≈ x

        x = fill(T(1), (2, 2))
        x1 = fill(AT{T}, T(1), (2, 2))
        x2 = fill(AT{T}, T(1), (2, 2))
        x3 = fill(AT{T, 2}, T(1), (2, 2))
        @test Array(x1) ≈ x
        @test Array(x2) ≈ x
        @test Array(x3) ≈ x

        fill!(x, 0)

        x1 = fill(AT{T}, 0f0, (2, 2))
        x2 = fill(AT{T}, T(0), (2, 2))
        x3 = fill(AT{T, 2}, T(0), (2, 2))
        @test Array(x1) ≈ x
        @test Array(x2) ≈ x
        @test Array(x3) ≈ x

        fill!(x1, 2f0)
        x2 = fill!(AT{Int32}(undef, (4, 4, 4)), 77f0)
        @test all(x-> x == 2f0, Array(x1))
        @test all(x-> x == Int32(77), Array(x2))

        x = fill(zero(T), (0, 2))
        x1 = fill(AT{T}, T(0), (0, 2))
        x2 = fill(AT{T}, T(0), (0, 2))
        x3 = fill(AT{T, 2}, T(0), (0, 2))
        @test Array(x1) ≈ x
        @test Array(x2) ≈ x
        @test Array(x3) ≈ x

        x = Matrix{T}(I, 4, 2)

        x1 = AT{T, 2}(I, 4, 2)
        x2 = AT{T}(I, (4, 2))
        x3 = AT{T, 2}(I, (4, 2))

        @test Array(x1) ≈ x
        @test Array(x2) ≈ x
        @test Array(x3) ≈ x

        x = Matrix(T(3) * I, 2, 4)
        x1 = AT(T(3) * I, 2, 4)
        @test eltype(x1) == T
        @test Array(x1) ≈ x

        x = fill(T(3), (2, 4))
        x1 = fill(AT{T}, T(3), (2, 4))
        copyto!(x, 2I)
        copyto!(x1, 2I)
        @test Array(x1) ≈ x
    end
end

@testsuite "iterator constructors" AT->begin
    for T in supported_eltypes()
        @test AT(Fill(T(0), (10,))) == fill(AT{T}, T(0), (10,))
        @test AT(Fill(T(0), (10, 10))) == fill(AT{T}, T(0), (10, 10))
        if T <: Real
            x = AT{Float32}(Fill(T(0), (10, 10)))
            @test eltype(x) == Float32
            @test AT(Eye{T}((10))) == AT{T}(I, 10, 10)
            x = AT{Float32}(Eye{T}(10))
            @test eltype(x) == Float32
        end
    end
end


text:
"""
An aggregator interface for SSA-like algorithms.

### Required Fields
- `next_jump`
- `next_jump_time`
- `end_time`
- `cur_rates`
- `sum_rate`
- `ma_jumps`
- `rates`
- `affects!`
- `save_positions`
- `rng`
"""
abstract type AbstractSSAJumpAggregator <: AbstractJumpAggregator end

DiscreteCallback(c::AbstractSSAJumpAggregator) = DiscreteCallback(c, c, initialize = c, save_positions = c.save_positions)

########### The following routines should be templates for all SSAs ###########

# # condition for jump to occur
# @inline function (p::SSAJumpAggregator)(u, t, integrator)
#   p.next_jump_time == t
# end

# # executing jump at the next jump time
# function (p::SSAJumpAggregator)(integrator)
#   execute_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
#   generate_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
#   register_next_jump_time!(integrator, p, integrator.t)
#   nothing
# end

# # setting up a new simulation
# function (p::SSAJumpAggregator)(dj, u, t, integrator) # initialize
#   initialize!(p, integrator, u, integrator.p, t)
#   register_next_jump_time!(integrator, p, integrator.t)
#   nothing
# end

############################## Generic Routines ###############################

@inline function register_next_jump_time!(integrator, p::AbstractSSAJumpAggregator, t)
  if p.next_jump_time < p.end_time
    add_tstop!(integrator, p.next_jump_time)
  end
  nothing
end

# helper routine for setting up standard fields of SSA jump aggregations
function build_jump_aggregation(jump_agg_type, u, p, t, end_time, ma_jumps, rates,
                                affects!, save_positions, rng; kwargs...)

  # mass action jumps
  majumps = ma_jumps
  if majumps === nothing
    majumps = MassActionJump(Vector{typeof(t)}(),
                             Vector{Vector{Pair{Int,eltype(u)}}}(),
                             Vector{Vector{Pair{Int,eltype(u)}}}())
  end

  # current jump rates, allows mass action rates and constant jumps
  cur_rates = Vector{typeof(t)}(undef,get_num_majumps(majumps) + length(rates))

  sum_rate = zero(typeof(t))
  next_jump = 0
  next_jump_time = typemax(typeof(t))
  jump_agg_type(next_jump, next_jump_time, end_time, cur_rates, sum_rate,
                majumps, rates, affects!, save_positions, rng; kwargs...)
end

# reevaluate all rates and total rate
function fill_rates_and_sum!(p::AbstractSSAJumpAggregator, u, params, t)
  sum_rate = zero(typeof(p.sum_rate))

  # mass action jumps
  majumps   = p.ma_jumps
  cur_rates = p.cur_rates
  @inbounds for i in 1:get_num_majumps(majumps)
      cur_rates[i] = evalrxrate(u, i, majumps)
      sum_rate    += cur_rates[i]
  end

  # constant rates
  rates = p.rates
  idx   = get_num_majumps(majumps) + 1
  @inbounds for rate in rates
      cur_rates[idx] = rate(u, params, t)
      sum_rate += cur_rates[idx]
      idx += 1
  end

  p.sum_rate = sum_rate
  nothing
end

# recalculate jump rates for jumps that depend on the just executed jump
# requires dependency graph
function update_dependent_rates!(p::AbstractSSAJumpAggregator, u, params, t)
  @inbounds dep_rxs = p.dep_gr[p.next_jump]
  num_majumps = get_num_majumps(p.ma_jumps)
  cur_rates   = p.cur_rates
  sum_rate    = p.sum_rate
  majumps     = p.ma_jumps
  @inbounds for rx in dep_rxs
      sum_rate -= cur_rates[rx]
      if rx <= num_majumps
          @inbounds cur_rates[rx] = evalrxrate(u, rx, majumps)
      else
          @inbounds cur_rates[rx] = p.rates[rx-num_majumps](u, params, t)
      end
      sum_rate += cur_rates[rx]
  end

  p.sum_rate = sum_rate
  nothing
end


text:
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random

function R_sir!(rate,xc,xd,parms,t::Float64,issum::Bool)
	(S,I,R,~) = xd
	(beta,mu) = parms
	infection = beta*S*I
	recovery = mu*I
	rate_display = 0.01
	if issum == false
			rate[1] = infection
			rate[2] = recovery
			rate[3] = rate_display
			return 0.
	else
		return infection+recovery + rate_display
	end
end

function F_sir!(xdot,xc,xd,parms,t::Float64)
	# vector field used for the continuous variable
	xdot[1] = 0.0
	nothing
end

xc0 = [0.0]
xd0 = [99,10,0,0]
nu = [[-1 1 0 0];[0 -1 1 0];[0 0 0 1]]
parms = [0.1/100.0,0.01]
tf = 150.0

Random.seed!(1234)
	problem = PDMP.PDMPProblem(F_sir!,R_sir!,nu, xc0, xd0, parms, (0.0, tf))
	result = PDMP.solve(problem, CHV(Tsit5()); n_jumps = 1000)


text:
function rejection(m)
    r = makeRand(m)
    function f(data)
        prop = r(;data...)
        while NamedTupleTools.select(prop, data) != data
            prop = r(;data...)
        end
        return prop
    end
end


text:
INTEGER_MOD(s::State) = if length(s.integer) >= 2
  divisor = pop!(s.integer)
  if divisor != zero(Int32)
    # FILL IN.
  end
end

INTEGER_MUL(s::State) = if length(s.integer) >= 2
  push!(s.integer, pop!(s.integer) * pop!(s.integer))
end

INTEGER_PLUS(s::State) = if length(s.integer) >= 2
  push!(s.integer, pop!(s.integer) + pop!(s.integer))
end

INTEGER_SUB(s::State) = if length(s.integer) >= 2
  push!(s.integer, - pop!(s.integer) + pop!(s.integer))
end

INTEGER_DIV(s::State) = if length(s.integer) >= 2
  divisor = peek(s.integer)
  if divisor != zero(Int32)
    pop!(s.integer)
    push!(s.integer, div(pop!(s.integer), divisor))
  end
end

INTEGER_LT(s::State) = if length(s.integer) >= 2
  push!(s.boolean, pop!(s.integer) > pop!(s.integer))
end

INTEGER_EQ(s::State) = if length(s.integer) >= 2
  push!(s.boolean, pop!(s.integer) == pop!(s.integer))
end

INTEGER_GT(s::State) = if length(s.integer) >= 2
  push!(s.boolean, pop!(s.integer) < pop!(s.integer))
end

INTEGER_DEFINE(s::State) = if !isempty(s.integer) && !isempty(s.name)
  s.definitions[pop!(s.name)] = pop!(s.integer)
end

INTEGER_DUP(s::State) = if !isempty(s.integer)
  push!(s.integer, peek(s.integer))
end

INTEGER_FLUSH(s::State) = empty!(s.integer)

INTEGER_FROM_BOOLEAN(s::State) = if !isempty(s.boolean)
  push!(s.integer, pop!(s.boolean) ? one(Int32) : zero(Int32))
end

#
# TODO:
# Check Push3 reference integer to floating point conversion semantics.
#
INTEGER_FROM_FLOAT(s::State) = if !isempty(s.float)
  push!(s.integer, pop!(s.float))
end

INTEGER_MAX(s::State) = if length(s.integer) >= 2
  push!(s.integer, max(pop!(s.integer), pop!(s.integer)))
end

INTEGER_MIN(s::State) = if length(s.integer) >= 2
  push!(s.integer, min(pop!(s.integer), pop!(s.integer)))
end

INTEGER_POP(s::State) = !isempty(s.integer) && pop!(s.integer)

INTEGER_RAND(s::State) = push!(s.integer, random_integer(s))

INTEGER_ROT(s::State) = if length(s.integer) > 2
  s.integer[end], s.integer[end-1], s.integer[end-2] =
    s.integer[end-1], s.integer[end-2], s.integer[end]
end

INTEGER_SHOVE(s::State) = if length(s.integer) >= 2
  shove!(s.integer, pop!(s.integer))
end

INTEGER_STACK_DEPTH(s::State) = push!(s.integer, length(s.integer))

INTEGER_SWAP(s::State) = if length(s.integer) >= 2
  s.integer[end], s.integer[end-1] = s.integer[end-1], s.integer[end]
end

INTEGER_YANK(s::State) = if length(s.integer) >= 2
  yank!(s.integer, pop!(s.integer))
end

INTEGER_YANKDUP(s::State) = if length(s.integer) >= 2
  yankdup!(s.integer, pop!(s.integer))
end

INTEGER_ERC(s::State) = return

Push.register("INTEGER.%",           INTEGER_MOD)
Push.register("INTEGER.*",           INTEGER_MUL)
Push.register("INTEGER.+",           INTEGER_PLUS)
Push.register("INTEGER.-",           INTEGER_SUB)
Push.register("INTEGER./",           INTEGER_DIV)
Push.register("INTEGER.<",           INTEGER_LT)
Push.register("INTEGER.=",           INTEGER_EQ)
Push.register("INTEGER.>",           INTEGER_GT)
Push.register("INTEGER.DEFINE",      INTEGER_DEFINE)
Push.register("INTEGER.DUP",         INTEGER_DUP)
Push.register("INTEGER.FLUSH",       INTEGER_FLUSH)
Push.register("INTEGER.FROMBOOLEAN", INTEGER_FROM_BOOLEAN)
Push.register("INTEGER.FROMFLOAT",   INTEGER_FROM_FLOAT)
Push.register("INTEGER.MAX",         INTEGER_MAX)
Push.register("INTEGER.MIN",         INTEGER_MIN)
Push.register("INTEGER.POP",         INTEGER_POP)
Push.register("INTEGER.RAND",        INTEGER_RAND)
Push.register("INTEGER.ROT",         INTEGER_ROT)
Push.register("INTEGER.SHOVE",       INTEGER_SHOVE)
Push.register("INTEGER.STACKDEPTH",  INTEGER_STACK_DEPTH)
Push.register("INTEGER.SWAP",        INTEGER_SWAP)
Push.register("INTEGER.YANK",        INTEGER_YANK)
Push.register("INTEGER.YANKDUP",     INTEGER_YANKDUP)
Push.register("INTEGER.ERC",         INTEGER_ERC)


text:
let
    v = collect(1:10)
    t = 5.5
    idx,Δt = linear_interp(v,t)
    @test idx == 5

end
let
    v = collect(1:10)
    t1 = 2.5
    t2 = 5.5
    idx1,Δt1 = linear_interp(v,t1)
    idx2,Δt2 = linear_interp(v,t2)
    vp = prune_sorted_array(v,idx1,idx2,Δt1,Δt2)
    @test length(vp) == 5
end


text:
module Util
    export human_readable_size, parse_size

    const ORDERS_OF_MAGNITUDE = "KMGTPEZY"

    """
    Convert `bytes` to a human-readable string.

    # Examples
    - 1024 → 1 KiB
    - 8589934592 → 8 GiB
    """
    function human_readable_size(bytes:: Integer) :: String
        order_of_magnitude = floor(Int, log(1024, bytes)) # Order of magnitude of bytes in base 1024 (0 for B, 1 for KiB, ...)
        order_of_magnitude = min(order_of_magnitude, length(ORDERS_OF_MAGNITUDE))
        bytes ÷= 1024^order_of_magnitude

        prefix = if order_of_magnitude == 0; ""
                 else ORDERS_OF_MAGNITUDE[order_of_magnitude] * "i"
                 end

        return "$bytes $(prefix)B"
    end

    """
    Convert human-readable string representing a number of bytes to that number of bytes.

    # Examples
    - 1KiB → 1024
    - 8GiB → 8589934592
    """
    function parse_size(size:: String) :: Integer
        if isnothing(match(r"^[0-9]+[KMGTPEZY]$", size))
            throw(ArgumentError("The given size string '$size' is invalid. Valid examples: 8G, 512M"))
        end

        magnitude_letter = size[end]
        magnitude = parse(Int, size[1:end-1])
        order_of_magnitude = findfirst(magnitude_letter, ORDERS_OF_MAGNITUDE)

        return magnitude * 1024^order_of_magnitude
    end
end

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
module testUtils

using Base.Test
using DataFrames
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running strategies tests\n")

tm = TimeData.readTimedata("/home/chris/.julia/v0.3/TimeData/data/logRet.csv")

####################
## random weights ##
####################

## randInvestments
wgts = AssetMgmt.randInvestments(tm)
@test isa(wgts, AssetMgmt.Investments)
@test AssetMgmt.idx(wgts) == idx(tm)
@test AssetMgmt.names(wgts) == names(tm)

######################
## equally weighted ##
######################

## equWgtInvestments
wgts = AssetMgmt.equWgtInvestments(tm)
@test isa(wgts, AssetMgmt.Investments)
@test AssetMgmt.idx(wgts) == idx(tm)
@test AssetMgmt.names(wgts) == names(tm)

wgts1 = AssetMgmt.core(wgts)[1, :]
wgts2 = AssetMgmt.core(wgts)[2, :]
@test (wgts1 == wgts2)

end


text:
export locally_free_class_group

################################################################################
#
#  Locally free class group
#
################################################################################

# Bley, Boltje "Computation of Locally Free Class Groups"
# If the left and right conductor of O in a maximal order coincide (which is the
# case if O is the integral group ring of a group algebra), the computation can
# be speeded up be setting cond = :left.
@doc Markdown.doc"""
    locally_free_class_group(O::AlgAssAbsOrd) -> GrpAbFinGen

> Given an order $O$ in a semisimple algebra over $\mathbb Q$, this function
> returns the locally free class group of $O$.
"""
function locally_free_class_group(O::AlgAssAbsOrd, cond::Symbol = :center)
  A = algebra(O)
  OA = maximal_order(O)
  Z, ZtoA = center(A)
  Fl = conductor(O, OA, :left)
  if cond == :left
    F = Fl
    FinZ = _as_ideal_of_smaller_algebra(ZtoA, F, OA)
  elseif cond == :center
    FinZ = _as_ideal_of_smaller_algebra(ZtoA, Fl, OA)
    # Compute FinZ*OA but as an ideal of O
    bOA = basis(OA, copy = false)
    bFinZ = basis(FinZ, copy = false)
    basis_F = Vector{elem_type(O)}()
    t = one(A)
    for x in bOA
      for y in bFinZ
        yy = ZtoA(elem_in_algebra(y, copy = false))
        t = mul!(t, yy, elem_in_algebra(x, copy = false))
        push!(basis_F, O(t))
      end
    end
    F = ideal_from_z_gens(O, basis_F)
  elseif cond == :product
    Fr = conductor(O, OA, :right)
    F = Fr*Fl
    FinZ = _as_ideal_of_smaller_algebra(ZtoA, F, OA)
  else
    error("Option :$(cond) for cond not implemented")
  end

  Adec = decompose(A)
  fields_and_maps = as_number_fields(Z)

  # Find the infinite places we need for the ray class group of FinZ
  inf_plc = Vector{Vector{InfPlc}}(undef, length(fields_and_maps))
  for i = 1:length(fields_and_maps)
    inf_plc[i] = Vector{InfPlc}()
  end
  for i = 1:length(Adec)
    B, BtoA = Adec[i]
    C, BtoC, CtoB = _as_algebra_over_center(B)
    K = base_ring(C)
    @assert K === fields_and_maps[i][1]

    places = real_places(K)
    for p in places
      if !issplit(C, p)
        push!(inf_plc[i], p)
      end
    end
  end

  R, mR = ray_class_group(FinZ, inf_plc)

  # Compute K_1(O/F) and the subgroup of R generated by nr(a)*OZ for a in k1 where
  # nr is the reduced norm and OZ the maximal order in Z
  k1 = K1_order_mod_conductor(O, OA, F, FinZ)

  k1_as_subgroup = Vector{elem_type(R)}()
  for x in k1
    # It is possible that x is not invertible in A
    t = isinvertible(elem_in_algebra(x, copy = false))[1]
    while !t
      r = rand(F, 100)
      x += r
      t = isinvertible(elem_in_algebra(x, copy = false))[1]
    end
    s = _reduced_norms(elem_in_algebra(x, copy = false), mR)
    push!(k1_as_subgroup, s)
  end

  Cl, CltoR = quo(R, k1_as_subgroup)

  return snf(Cl)[1]
end

# Helper function for locally_free_class_group
# Computes the representative in the ray class group (domain(mR)) for the ideal
# nr(a)*O_Z, where nr is the reduced norm and O_Z the maximal order of the centre
# of A.
function _reduced_norms(a::AlgAssElem, mR::MapRayClassGroupAlg)
  A = parent(a)
  Adec = decompose(A)
  r = zero_matrix(FlintZZ, 1, 0)

  for i = 1:length(Adec)
    B, BtoA = Adec[i]
    C, BtoC, CtoB = _as_algebra_over_center(B)
    c = BtoC(BtoA\a)
    G, GtoIdl = mR.groups_in_number_fields[i]
    K = number_field(order(codomain(GtoIdl)))
    OK = maximal_order(K)
    @assert K === base_ring(C)
    nc = norm(c)
    I = OK(nc)*OK
    m = isqrt(dim(C))
    @assert m^2 == dim(C)
    b, J = ispower(I, m)
    @assert b
    g = GtoIdl\J
    r = hcat(r, g.coeff)
  end
  G = codomain(mR.into_product_of_groups)
  return mR.into_product_of_groups\G(r)
end

################################################################################
#
#  K1
#
################################################################################

# Computes generators for K_1(O/F) where F is the product of the left and right
# conductor of O in the maximal order.
# FinZ should be F intersected with the centre of algebra(O).
# See Bley, Boltje "Computation of Locally Free Class Groups"
function K1_order_mod_conductor(O::AlgAssAbsOrd, OA::AlgAssAbsOrd, F::AlgAssAbsOrdIdl, FinZ::AlgAssAbsOrdIdl)
  A = algebra(O)
  Z, ZtoA = center(A)
  OZ = maximal_order(Z)
  OinZ = _as_order_of_smaller_algebra(ZtoA, O, OA)

  facFinZ = factor(FinZ)
  prime_ideals = Dict{ideal_type(OinZ), Vector{ideal_type(OZ)}}()
  for (p, e) in facFinZ
    q = contract(p, OinZ)
    if haskey(prime_ideals, q)
      push!(prime_ideals[q], p)
    else
      prime_ideals[q] = [ p ]
    end
  end

  primary_ideals = Vector{Tuple{ideal_type(O), ideal_type(O)}}()
  for p in keys(prime_ideals)
    primes_above = prime_ideals[p]
    q = primes_above[1]^facFinZ[primes_above[1]]
    for i = 2:length(primes_above)
      q = q*primes_above[i]^facFinZ[primes_above[i]]
    end
    pO = _as_ideal_of_larger_algebra(ZtoA, p, O)
    qO = _as_ideal_of_larger_algebra(ZtoA, contract(q, OinZ), O)
    # The qO are primary ideals such that F = \prod (qO + F)
    push!(primary_ideals, (pO, qO))
  end

  moduli = Vector{ideal_type(O)}()
  for i = 1:length(primary_ideals)
    push!(moduli, primary_ideals[i][2] + F)
  end

  # Compute generators of K_1(O/q + F) for each q and put them together with the CRT
  elements_for_crt = Vector{Vector{elem_type(O)}}(undef, length(primary_ideals))
  for i = 1:length(primary_ideals)
    # We use the exact sequence
    # (1 + p + F)/(1 + q + F) -> K_1(O/q + F) -> K_1(O/p + F) -> 1
    (p, q) = primary_ideals[i]
    pF = p + F
    qF = moduli[i]
    char = basis_matrix(p, copy = false)[1, 1]
    B, OtoB = AlgAss(O, pF, char)
    k1_B = K1(B)
    k1_O = [ OtoB\x for x in k1_B ]
    if pF != qF
      append!(k1_O, _1_plus_p_mod_1_plus_q_generators(pF, qF))
    end
    elements_for_crt[i] = k1_O
  end
  # Make the generators coprime to the other ideals
  if length(moduli) != 0 # maybe O is maximal
    k1 = make_coprime(elements_for_crt, moduli)
  else
    k1 = elem_type(O)[]
  end

  return k1
end

@doc Markdown.doc"""
    K1(A::AlgAss{T}) where { T <: Union{gfp_elem, Generic.ResF{fmpz}, fq, fq_nmod } }
      -> Vector{AbsAlgAssElem}

> Given an algebra over a finite field, this function returns generators for $K_1(A)$.
"""
function K1(A::AlgAss{T}) where { T <: Union{gfp_elem, Generic.ResF{fmpz}, fq, fq_nmod } }
  # We use the exact sequence 1 + J -> K_1(A) -> K_1(B/J) -> 1
  J = radical(A)
  onePlusJ = _1_plus_j(A, J)

  B, AtoB = quo(A, J)
  k1B = K1_semisimple(B)
  k1 = append!(onePlusJ, [ AtoB\x for x in k1B ])
  return k1
end

# Computes generators for K_1(A) with A semisimple as described in
# Bley, Boltje "Computation of Locally Free Class Groups", p. 84.
function K1_semisimple(A::AlgAss{T}) where { T <: Union{ gfp_elem, Generic.ResF{fmpz}, fq, fq_nmod } }

  Adec = decompose(A)
  k1 = Vector{elem_type(A)}()
  idems = [ BtoA(one(B)) for (B, BtoA) in Adec ]
  sum_idems = sum(idems)
  minus_idems = map(x -> -one(A)*x, idems)
  for i = 1:length(Adec)
    B, BtoA = Adec[i]
    C, BtoC, CtoB = _as_algebra_over_center(B)
    F = base_ring(C)
    # Consider C as a matrix algebra over F. Then the matrices with a one somewhere
    # on the diagonal are given by primitive idempotents (see also _as_matrix_algebra).
    prim_idems = _primitive_idempotents(C)
    a = _primitive_element(F)
    # aC is the identity matrix with a at position (1, 1)
    aC = a*prim_idems[1]
    if dim(C) > 1
      for j = 2:length(prim_idems)
        aC = add!(aC, aC, prim_idems[j])
      end
    end
    aA = BtoA(CtoB(aC))
    # In the other components aA should be 1 (this is not mentioned in the Bley/Boltje-Paper)
    aA = add!(aA, aA, sum_idems)
    aA = add!(aA, aA, minus_idems[i])
    push!(k1, aA)
  end
  return k1
end

# Computes generators for 1 + J where J is the Jacobson Radical of A
function _1_plus_j(A::AlgAss{T}, jacobson_radical::AbsAlgAssIdl...) where { T <: Union{ gfp_elem, Generic.ResF{fmpz}, fq_nmod, fq } }
  F = base_ring(A)

  if length(jacobson_radical) == 1
    J = jacobson_radical[1]
  else
    J = radical(A)
  end

  onePlusJ = Vector{elem_type(A)}()

  if iszero(J)
    return onePlusJ
  end

  # We use the filtration 1 + J \supseteq 1 + J^2 \subseteq ... \subseteq 1
  oneA = one(A)
  while !iszero(J)
    J2 = J^2
    Q, AtoQ = quo(J, J2)
    for i = 1:dim(Q)
      push!(onePlusJ, one(A) + AtoQ\Q[i])
    end
    J = J2
  end
  return onePlusJ
end


text:
using Revise, Parameters, Logging
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\kernelEvolutions.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\structs\\BasicStructs.jl")

includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\kernels\\TpfpfnKernel.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\overLap\\MainOverlap.jl")
using ..BasicPreds, ..CUDAGpuUtils,Cthulhu,BenchmarkTools , CUDA
using ..MainOverlap, ..TpfpfnKernel

goldBoolGPU,segmBoolGPU,tp,tn,fp,fn, tpArr,tnArr,fpArr, fnArr, blockNum , nx,ny,nz ,tpTotalTrue,tnTotalTrue,fpTotalTrue, fnTotalTrue ,tpPerSliceTrue,  tnPerSliceTrue,fpPerSliceTrue,fnPerSliceTrue ,flattG, flattSeg ,FlattGoldGPU,FlattSegGPU,intermediateResTp,intermediateResFp,intermediateResFn = getSmallTestBools();

sizz= size(goldBoolGPU)

TpfpfnKernel.getTpfpfnData!(FlattGoldGPU,FlattSegGPU,tp,tn,fp,fn, intermediateResTp,intermediateResFp,intermediateResFn,sizz[1],sizz[1]*sizz[1],sizz[3],Float32(1))
tp[1]
tp[1] ==tpTotalTrue && fp[1] ==fpTotalTrue && fn[1] ==fnTotalTrue #tn[1] == tnTotalTrue &&

using CUDA

#@benchmark CUDA.@sync
# we are adding false in the end to make indexing easier
blockss = Int64(round((length(FlattGoldGPU)/512)/16))

#@benchmark CUDA.@sync
args = (FlattGoldGPU,FlattSegGPU,tp,tn,fp,fn, intermediateResTp,intermediateResFp,intermediateResFn, loopNumb, indexCorr)
 @cuda threads=512 blocks=blockNum-1 getBlockTpFpFn(args...)
tp[1]
fp[1]

tp[1] ==tpTotalTrue && fp[1] ==fpTotalTrue && fn[1] ==fnTotalTrue #tn[1] == tnTotalTrue &&

tp[1]
fp[1]
fn[1]

sum(intermediateResTp)
sum(intermediateResFp)
sum(intermediateResFn)

sum(FlattGoldGPU)

tpNew = (9*9*9)*2 +18

  @device_code_warntype interactive=true @cuda getBlockTpFpFn(FlattGoldGPU,FlattSegGPU,tp,tn,fp,fn, intermediateResTp,intermediateResFp,intermediateResFn)

## using occupancy API to calculate  threads number, block number etc...
#args = (FlattGoldGPU,FlattSegGPU,tp,tn,fp,fn, intermediateResTp,intermediateResFp,intermediateResFn)
calcBlocks,valcThreads,maxBlocks = computeBlocksFromOccupancy(args,3 )