text:
# # Bukdu v0.4.5
using Bukdu # ApplicationController Conn HTML Changeset render routes get post
using Bukdu.HTML5.Form # form_for file_input submit
using Documenter.Utilities.DOM: @tags
using HTTP: Multipart
using Base64

struct FormController <: ApplicationController
    conn::Conn
end

function layout(body)
    """

  <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />

Multipart example

$body

"""
end

function preview(multipart::Multipart)
    if eof(multipart)
        string("empty")
    else
        @tags div img
        mark(multipart)
        data = read(multipart)
        reset(multipart)
        contents = [
            div(multipart.filename),
            div(repr(length(data)), " bytes"),
        ]
        if startswith(multipart.contenttype, "image")
            base64 = base64encode(data)
            push!(contents, img[:src => string("data:", multipart.contenttype, ";base64, ", base64),
                                :alt => multipart.filename])
        end
        div(contents)
    end
end

function preview(::Nothing)
    ""
end

struct User
end

function post_result(c::FormController)
    @tags div h3
    user_photo1 = c.params["user_photo1"]
    user_photo2 = c.params["user_photo2"]
    body = div([
        h3("user_photo1"),
        preview(user_photo1),
        h3("user_photo2"),
        preview(user_photo2),
    ])
    render(HTML, layout(body))
end

function index(c::FormController)
    @tags div
    form1 = form_for(Changeset(User), (FormController, post_result), method=post, multipart=true) do f
        div(
            div.([
                label_for(file_input(f, :photo1), "Photo1"),
                label_for(file_input(f, :photo2), "Photo2"),
                submit("Submit"),
            ])
        )
    end
    body = div(form1)
    render(HTML, layout(body))
end

routes() do
    get("/", FormController, index)
    post("/post_result", FormController, post_result)
end

Bukdu.start(8080)
Router.call(get, "/") #
Base.JLOptions().isinteractive==0 && wait()


text:
function gnomesort!(arr::Vector)
    i, j = 1, 2
    while i < length(arr)
        if arr[i] ≤ arr[i+1]
            i = j
            j += 1
        else
            arr[i], arr[i+1] = arr[i+1], arr[i]
            i -= 1
            if i == 0
                i = j
                j += 1
            end
        end
    end
    return arr
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", bogosort!(v))

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
macro c(ex)
    Meta.isexpr(ex, :call) || throw(ArgumentError("not a function call expression."))
    prologue = Expr(:block)
    epilogue = Expr(:block)
    func_expr = Expr(:call, first(ex.args))
    for arg in ex.args[2:end]
        if Meta.isexpr(arg, :&)
            refee = arg.args[]
            if Meta.isexpr(refee, :ref) && length(refee.args) == 2
                # &a[n] => pointer(a) + n * Core.sizeof(eltype(a))
                array_name, n = refee.args
                push!(func_expr.args, :(pointer($array_name) + $n * Core.sizeof(eltype($array_name))))
            else
                ref_sym = gensym("cref")
                push!(prologue.args, Expr(:(=), ref_sym, Expr(:call, :Ref, refee)))
                push!(func_expr.args, ref_sym)
                push!(epilogue.args, Expr(:(=), refee, Expr(:ref, ref_sym)))
            end
        else
            push!(func_expr.args, arg)
        end
    end
    func_ret = gensym("cref_ret")
    push!(prologue.args, Expr(:(=), func_ret, func_expr))
    append!(prologue.args, epilogue.args)
    push!(prologue.args, func_ret)
    return esc(prologue)
end

function glfwGetMonitors(count)
    ccall((:glfwGetMonitors, GLFW.libglfw), Ptr{Ptr{Cvoid}}, (Ptr{Cint},), count)
end

function glfwGetMonitorPos(monitor, xpos, ypos)
    ccall((:glfwGetMonitorPos, GLFW.libglfw), Cvoid, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}), monitor, xpos, ypos)
end

struct GLFWvidmode
    width::Cint
    height::Cint
    redBits::Cint
    greenBits::Cint
    blueBits::Cint
    refreshRate::Cint
end

function glfwGetVideoMode(monitor)
    ccall((:glfwGetVideoMode, GLFW.libglfw), Ptr{GLFWvidmode}, (Ptr{Cvoid},), monitor)
end

function glfwGetMonitorWorkarea(monitor, xpos, ypos, width, height)
    ccall((:glfwGetMonitorWorkarea, GLFW.libglfw), Cvoid, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), monitor, xpos, ypos, width, height)
end

function glfwGetMonitorContentScale(monitor, xscale, yscale)
    ccall((:glfwGetMonitorContentScale, GLFW.libglfw), Cvoid, (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}), monitor, xscale, yscale)
end

function glfwGetWindowSize(window, width, height)
    ccall((:glfwGetWindowSize, GLFW.libglfw), Cvoid, (GLFW.Window, Ptr{Cint}, Ptr{Cint}), window, width, height)
end

function glfwGetFramebufferSize(window, width, height)
    ccall((:glfwGetFramebufferSize, GLFW.libglfw), Cvoid, (GLFW.Window, Ptr{Cint}, Ptr{Cint}), window, width, height)
end

function glfwGetTime()
    ccall((:glfwGetTime, GLFW.libglfw), Cdouble, ())
end

function glfwSetWindowTitle(window, title)
    ccall((:glfwSetWindowTitle, GLFW.libglfw), Cvoid, (GLFW.Window, Ptr{Cchar}), window, title)
end

function glfwSetMouseButtonCallback(window, callback)
    ccall((:glfwSetMouseButtonCallback, GLFW.libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), window, callback)
end

function glfwSetScrollCallback(window, callback)
    ccall((:glfwSetScrollCallback, GLFW.libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), window, callback)
end

function glfwSetKeyCallback(window, callback)
    ccall((:glfwSetKeyCallback, GLFW.libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), window, callback)
end

function glfwSetCharCallback(window, callback)
    ccall((:glfwSetCharCallback, GLFW.libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), window, callback)
end

function glfwSetMonitorCallback(callback)
    ccall((:glfwSetMonitorCallback, GLFW.libglfw), Ptr{Cvoid}, (Ptr{Cvoid},), callback)
end


text:
#!/usr/bin/env julia

# This file is a part of SIS3316.jl, licensed under the MIT License (MIT).

using ArgParse
using SIS3316
using ROOTFramework

function sis3316_to_root end
export sis3316_to_root

mktemp_custom(parent=tempdir(), format="tmpXXXXXX") = begin
    b = joinpath(parent, format)
    p = ccall(:mkstemp, Int32, (Ptr{UInt8},), b) # modifies b
    systemerror(:mktemp, p == -1)
    return (b, fdio(p, true))
end

sis3316_to_root(input_io::IO, output_tfile::TFile; evt_merge_window::AbstractFloat = 100e-9, waveform_format = :none) = begin
    _time(x::Pair{Int64, SIS3316.RawChEvent}) = time(x.second)

    bindings = TTreeBindings()

    const info_idx = bindings[:info_idx] = Ref{Int32}(0)
    const info_time = bindings[:info_time] = Ref{Float64}(0)

    const raw_pp_ch = bindings[:raw_pp_ch] = Int32[]
    const raw_pp_mca = bindings[:raw_pp_mca] = Int32[]
    const raw_pp_trig_max = bindings[:raw_pp_trig_max] = Int32[]
    const raw_pp_peak_pos = bindings[:raw_pp_peak_pos] = Int32[]
    const raw_pp_peak_height = bindings[:raw_pp_peak_height] = Int32[]
    const raw_pp_acc = map(i -> bindings[Symbol("raw_pp_acc_$i")] = Int32[], 1:8)

    const raw_trig_ch = bindings[:raw_trig_ch] = Int32[]
    const raw_trig_trel = bindings[:raw_trig_trel] = Float64[]
    const raw_trig_pileup = bindings[:raw_trig_pileup] = Int32[]
    const raw_trig_overflow = bindings[:raw_trig_overflow] = Int32[]

    const raw_wf_ch = Int32[]
    const raw_wf_smpl_n = Int32[]
    const raw_wf_smpl_v = Int32[]

    if waveform_format == :none
    elseif waveform_format == :integers
        bindings[:raw_wf_ch] = raw_wf_ch
        bindings[:raw_wf_smpl_n] = raw_wf_smpl_n
        bindings[:raw_wf_smpl_v] = raw_wf_smpl_v
    else
        error("Unknown waveform format $waveform_format.")
    end

    const ch_sized_vecs = Vector[raw_pp_ch, raw_pp_mca, raw_pp_trig_max,
        raw_pp_peak_pos, raw_pp_peak_height,
        raw_trig_ch, raw_trig_trel, raw_trig_pileup, raw_trig_overflow
    ]
    append!(ch_sized_vecs, raw_pp_acc)

    for v in ch_sized_vecs sizehint!(v, 16) end

    ttree = create_ttree!(output_tfile, "data", "Data")
    ttoutput = TTreeOutput(ttree, bindings)

    const reader = eachchunk(input_io, SIS3316.UnsortedEvents)

    local evtno = 0

    for unsorted in reader
        const sorted = sortevents(unsorted, merge_window = evt_merge_window)
        const evtv = Vector{Pair{Int64, SIS3316.RawChEvent}}()
        const timestamps = Vector{Float64}()

        const energynull = SIS3316.EnergyValues(0, 0)
        const mawnull = SIS3316.MAWValues(0, 0, 0)
        const psanull = SIS3316.PSAValue(0, 0)
        const flagsnull = SIS3316.EvtFlags(false,false,false,false)

        for evt in sorted
            evtno += 1

            resize!(evtv, length(evt))
            copy!(evtv, evt)
            sort!(evtv, by = first)
            resize!(timestamps, length(evtv))
            map!(_time, timestamps, evtv)
            const starttime = isempty(timestamps) ? zero(Float64) : minimum(timestamps)

            info_idx.x = evtno
            info_time.x = starttime

            for v in ch_sized_vecs empty!(v) end

            empty!(raw_wf_ch)
            empty!(raw_wf_smpl_n)
            empty!(raw_wf_smpl_v)

            for (ch, chevt) in evtv
                push!(raw_pp_ch, ch)
                push!(raw_pp_mca, get(chevt.energy, energynull).maximum)
                push!(raw_pp_trig_max, get(chevt.trig_maw, mawnull).maximum)
                push!(raw_pp_peak_pos, get(chevt.peak_height, psanull).index)
                push!(raw_pp_peak_height, get(chevt.peak_height, psanull).value)

                for i in eachindex(raw_pp_acc)
                    push!(raw_pp_acc[i], get(chevt.accsums, i, 0))
                end

                push!(raw_trig_ch, ch)
                push!(raw_trig_trel, time(chevt) - starttime)
                push!(raw_trig_pileup, chevt.pileup_flag + 2 * get(chevt.flags, flagsnull).pileup +  4 * get(chevt.flags, flagsnull).repileup)
                push!(raw_trig_overflow, 1 * get(chevt.flags, flagsnull).overflow +  2 * get(chevt.flags, flagsnull).underflow)

                if !isempty(chevt.samples)
                    if waveform_format == :integers
                        push!(raw_wf_ch, ch)
                        push!(raw_wf_smpl_n, length(chevt.samples))
                        append!(raw_wf_smpl_v, chevt.samples)
                    end
                end
            end

            push!(ttoutput)
        end
    end
end

sis3316_to_root(input_fname::AbstractString; evt_merge_window::AbstractFloat = 100e-9, waveform_format = :none) = begin
    const fnexpr = r"(.*)\.dat(\.[^.]+)?"
    const fnbase = match(fnexpr, basename(input_fname))[1]
    const output_fname = "$(fnbase).root"

    if !ispath(output_fname)
        output_tmpname, tmpio = mktemp_custom(pwd(), "$(output_fname).tmp-XXXXXX")
        close(tmpio)

        open(SIS3316.DecompessIO, input_fname) do input_io
            open(TFile, output_tmpname, "recreate") do output_tfile
                sis3316_to_root(input_io, output_tfile, evt_merge_window = evt_merge_window, waveform_format = waveform_format)
            end
        end

        mv(output_tmpname, output_fname, remove_destination = false)
    else
        info("Output file \"$(output_fname)\" already exists, skipping \"$(input_fname)\".")
    end
end

sis3316_to_root() = begin
    argsettings = ArgParseSettings()

    @add_arg_table argsettings begin
        "--merge-within", "-m"
            help = "time window for merging events (no merging for negative values)"
            arg_type = Float64
            default = 100e-9
        "--waveform-format", "-w"
            help = "waveform format (\"none\" or \"integers\")"
            arg_type = AbstractString
            default = "integers"
        "inputs"
            help = "input files"
            nargs = '*'
    end

    const parsed_args = parse_args(argsettings)

    for input_fname in parsed_args["inputs"]
        try
            info("Converting \"$(input_fname)\"")
            @time sis3316_to_root(
                input_fname,
                evt_merge_window = parsed_args["merge-within"],
                waveform_format = Symbol(parsed_args["waveform-format"])
            )
        catch err
            print_with_color(:red, STDERR, "ERROR: $err\n")
        end
    end
end


text:
mutable struct Topology
    topology::Union{Nothing,Hwloc.Object}
end
Topology() = Topology(nothing)
function safe_topology_load!()
    try
        TOPOLOGY.topology = Hwloc.topology_load();
    catch e
        @warn e
        @warn """
            Using Hwloc failed. Please file an issue with the above warning at: https://github.com/JuliaParallel/Hwloc.jl
            Proceeding with generic topology assumptions. This may result in reduced performance.
        """
    end
end

const TOPOLOGY = Topology()
safe_topology_load!()

function count_attr(topology::Hwloc.Object, attr::Symbol)
    count = 0
    for t ∈ topology
        count += t.type_ == attr
    end
    count
end

function count_attr(attr::Symbol)
    topology = TOPOLOGY.topology
    topology === nothing && return nothing
    count_attr(topology, attr)
end

function define_attr_count(fname::Symbol, v)
    if v === nothing
        @eval $fname() = nothing
        return
    elseif v isa Integer
        @eval $fname() = StaticInt{$(convert(Int,v))}()
    elseif v isa Bool
        if v
            @eval $fname() = True()
        else
            @eval $fname() = False()
        end
    else
        @eval $fname() = $v
    end
    nothing
end
# define_attr(fname::Symbol, attr::Symbol) = define_attr_count(fname, )

for (f, attr) ∈ [
    (:num_l1cache, :L1Cache),
    (:num_l2cache, :L2Cache),
    (:num_l3cache, :L3Cache),
    (:num_l4cache, :L4Cache),
    (:num_machines, :Machine),
    (:num_sockets, :Package),
    (:num_cores, :Core),
    (:sys_threads, :PU),
    (:num_threads, :PU)
]
    define_attr_count(f, count_attr(attr))
end

function redefine_attr_count()
    sys_thread = sys_threads()
    iter = [
        (num_l1cache(), :num_l1cache, :L1Cache),
        (num_l2cache(), :num_l2cache, :L2Cache),
        (num_l3cache(), :num_l3cache, :L3Cache),
        (num_l4cache(), :num_l4cache, :L4Cache),
        (num_machines(), :num_machines, :Machine),
        (num_sockets(), :num_sockets, :Package),
        (num_cores(), :num_cores, :Core),
        (sys_thread, :sys_threads, :PU)
    ]
    for (v, f, attr) ∈ iter
        ref = count_attr(attr)
        if ref ≠ v
            @info "Redefining attr count $f = $ref."
            define_attr_count(f, ref)
        end
    end
    if sys_thread != Threads.nthreads()
        @eval num_threads() = StaticInt{$(Threads.nthreads())}()
    end
    nothing
end

num_cache(::Union{Val{1},StaticInt{1}}) = num_l1cache()
num_cache(::Union{Val{2},StaticInt{2}}) = num_l2cache()
num_cache(::Union{Val{3},StaticInt{3}}) = num_l3cache()
num_cache(::Union{Val{4},StaticInt{4}}) = num_l4cache()

function num_cache_levels()
    numl4 = num_l4cache()
    numl4 === nothing && return nothing
    ifelse(
        eq(numl4, Zero()),
        ifelse(
            eq(num_l3cache(), Zero()),
            ifelse(
                eq(num_l2cache(), Zero()),
                ifelse(
                    eq(num_l1cache(), Zero()),
                    Zero(),
                    One()
                ),
                StaticInt{2}()
            ),
            StaticInt{3}()
        ),
        StaticInt{4}()
    )
end

function dynamic_cache_inclusivity()::NTuple{4,Bool}
    if !((Sys.ARCH === :x86_64) || (Sys.ARCH === :i686))
        return (false,false,false,false)
    end
    function get_cache_edx(subleaf)
        # source: https://github.com/m-j-w/CpuId.jl/blob/401b638cb5a020557bce7daaf130963fb9c915f0/src/CpuInstructions.jl#L38
        # credit Markus J. Weber, copyright: https://github.com/m-j-w/CpuId.jl/blob/master/LICENSE.md
        Base.llvmcall(
            """
            ; leaf = %0, subleaf = %1, %2 is some label
            ; call 'cpuid' with arguments loaded into registers EAX = leaf, ECX = subleaf
            %2 = tail call { i32, i32, i32, i32 } asm sideeffect "cpuid",
                "={ax},={bx},={cx},={dx},{ax},{cx},~{dirflag},~{fpsr},~{flags}"
                (i32 4, i32 %0) #2
            ; retrieve the result values and return eax and edx contents
            %3 = extractvalue { i32, i32, i32, i32 } %2, 0
            %4 = extractvalue { i32, i32, i32, i32 } %2, 3
            %5  = insertvalue [2 x i32] undef, i32 %3, 0
            %6  = insertvalue [2 x i32]   %5 , i32 %4, 1
            ; return the value
            ret [2 x i32] %6
            """
            # llvmcall requires actual types, rather than the usual (...) tuple
            , Tuple{UInt32,UInt32}, Tuple{UInt32}, subleaf % UInt32
        )
    end
    # eax0, edx1 = get_cache_edx(0x00000000)
    t = (false,false,false,false)
    i = zero(UInt32)
    j = 0
    while (j < 4)
        eax, edx = get_cache_edx(i)
        i += one(UInt32)
        iszero(eax & 0x1f) && break
        iszero(eax & 0x01) && continue
        ci = ((edx & 0x00000002) != 0x00000000) & (eax & 0x1f != 0x00000000)
        t = Base.setindex(t, ci, (j += 1))
    end
    t
end

nothing_cache_summary() = (size = 0, linesize = 64, associativity = nothing, type = nothing, inclusive = nothing)
function dynamic_cache_summary(N)
    topology = TOPOLOGY.topology
    cache_name = (:L1Cache, :L2Cache, :L3Cache, :L4Cache)[N]
    if (topology === nothing) || (iszero(count_attr(cache_name)))
        return nothing_cache_summary()
    end
    c = first(t for t in topology if t.type_ == cache_name && t.attr.depth == N).attr
    (
        size = c.size,
        linesize = c.linesize,
        associativity = c.associativity,
        type = c.type_,
        inclusive = dynamic_cache_inclusivity()[N]
    )
end
cache_size(_) = StaticInt{0}()
cache_linesize(_) = StaticInt{64}() # assume...
cache_associativity(_) = nothing
cache_type(_) = nothing
cache_inclusive(_) = nothing
function define_cache(N, c = dynamic_cache_summary(N))
    c === nothing_cache_summary() || _define_cache(N, c)
end
function _define_cache(N, c)
    @eval begin
        cache_size(::Union{Val{$N},StaticInt{$N}}) = StaticInt{$(c.size)}()
        cache_linesize(::Union{Val{$N},StaticInt{$N}}) = StaticInt{$(c.linesize)}()
        cache_associativity(::Union{Val{$N},StaticInt{$N}}) = StaticInt{$(c.associativity)}()
        cache_type(::Union{Val{$N},StaticInt{$N}}) = Val{$(c.type === nothing ? nothing : QuoteNode(c.type))}()
        cache_inclusive(::Union{Val{$N},StaticInt{$N}}) = $(c.inclusive ? :True : :False)()
    end
    nothing
end

function redefine_cache(N)
    s = cache_size(StaticInt(N))
    l = cache_linesize(StaticInt(N))
    a = cache_associativity(StaticInt(N))
    t = cache_type(StaticInt(N))
    i = cache_inclusive(StaticInt(N))
    c = (
        size = Int(s)::Int,
        linesize = l === nothing ? nothing : Int(l)::Int,
        associativity = a === nothing ? nothing : Int(a)::Int,
        type = t === nothing ? nothing : unwrap(t)::Symbol,
        inclusive = i === nothing ? nothing : Bool(i)::Bool
    )
    correct = dynamic_cache_summary(N)
    if c !== correct
        @info "Redefining cache $N."
        _define_cache(N, correct)
    end
    nothing
end
foreach(define_cache, 1:4)

cache_linesize() = cache_linesize(Val(1))


text:
using st16SseqJuliaTools
using Test
using BioTools.BLAST
using BioSequences
using IntervalSets
@testset "Blast results manipulation" begin
    playxml = "test/playdata/chimrem/blast_results.xml"
    data = split_by_query(readblastXML(read(playxml, String)))
    @test [10..16, 17..19] == add_new_interval([10..12,15..16,17..19],12..15)
    @test [10..12, 13..14, 15..16, 17..19] == add_new_interval([10..12,15..16,17..19],13..14)
    @test covered_length([10..16,17..19]) == 10
    @test 0.0 == unique_length_fraction([10..12, 13..14, 15..16, 17..19],12..13)
    @test 0.5 == unique_length_fraction([10..12, 13..14, 15..16, 17..19],19..20)
    @test [("gi|1384345818|gb|CP029032.1|", 556866, 556195), ("gi|1384345818|gb|CP029032.1|", 1997488, 1997716)] == map(x -> (x.hitid,x.hitfrom,x.hitto),get_coverage_maximising_hits("1_1828;size=1",data))
end


text:
"""
    arrays(f::ROOTFile, treename)

Reads all branches from a tree.
"""
function arrays(f::ROOTFile, treename)
    names = keys(f[treename])
    res = Vector{Any}(undef, length(names))
    Threads.@threads for i in eachindex(names)
        res[i] = array(f, "$treename/$(names[i])")
    end
    return res
end

"""
    array(f::ROOTFile, path; raw=false)

Reads an array from a branch. Set `raw=true` to return raw data and correct offsets.
"""
function array(f::ROOTFile, path::AbstractString; raw=false)
    return array(f::ROOTFile, _getindex(f, path); raw=raw)
end

function array(f::ROOTFile, branch; raw=false)
    ismissing(branch) && error("No branch found at $path")
    (!raw && length(branch.fLeaves.elements) > 1) && error(
        "Branches with multiple leaves are not supported yet. Try reading with `array(...; raw=true)`.",
    )

    rawdata, rawoffsets = readbranchraw(f, branch)
    if raw
        return rawdata, rawoffsets
    end
    T, J = auto_T_JaggT(f, branch; customstructs=f.customstructs)
    return interped_data(rawdata, rawoffsets, T, J)
end

"""
    basketarray(f::ROOTFile, path::AbstractString, ith)
    basketarray(f::ROOTFile, branch::Union{TBranch, TBranchElement}, ith)

Reads actual data from ith basket of a branch. This function first calls [`readbasket`](@ref)
to obtain raw bytes and offsets of a basket, then calls [`auto_T_JaggT`](@ref) followed
by [`interped_data`](@ref) to translate raw bytes into actual data.
"""
function basketarray(f::ROOTFile, path::AbstractString, ithbasket)
    return basketarray(f, f[path], ithbasket)
end
function basketarray(f::ROOTFile, branch, ithbasket)
    ismissing(branch) && error("No branch found at $path")
    length(branch.fLeaves.elements) > 1 && error(
        "Branches with multiple leaves are not supported yet. Try reading with `array(...; raw=true)`.",
    )

    rawdata, rawoffsets = readbasket(f, branch, ithbasket)
    T, J = auto_T_JaggT(f, branch; customstructs=f.customstructs)
    return interped_data(rawdata, rawoffsets, T, J)
end

# function barrior to make getting individual index faster
# TODO upstream some types into parametric types for Branch/BranchElement
"""
    LazyBranch(f::ROOTFile, branch)

Construct an accessor for a given branch such that `BA[idx]` and or `BA[1:20]` is
type-stable. And memory footprint is a single basket (<1MB usually). You can also
iterate or map over it. If you want a concrete `Vector`, simply `collect()` the
LazyBranch.

# Example
```julia
julia> rf = ROOTFile("./test/samples/tree_with_large_array.root");

julia> b = rf["t1/int32_array"];

julia> ab = UnROOT.LazyBranch(rf, b);

julia> for entry in ab
           @show entry
           break
       end
entry = 0

julia> ab[begin:end]
0
1
...
```
"""
mutable struct LazyBranch{T,J,B} <: AbstractVector{T}
    f::ROOTFile
    b::Union{TBranch,TBranchElement}
    L::Int64
    fEntry::Vector{Int64}
    buffer::Vector{B}
    buffer_range::Vector{UnitRange{Int64}}

    function LazyBranch(f::ROOTFile, b::Union{TBranch,TBranchElement})
        T, J = auto_T_JaggT(f, b; customstructs=f.customstructs)
        _buffer = J === Nojagg ? T[] : VectorOfVectors{eltype(T)}()
        return new{T,J,typeof(_buffer)}(f, b, length(b),
                                        b.fBasketEntry,
                                        [_buffer for _ in 1:Threads.nthreads()],
                                        [0:-1 for _ in 1:Threads.nthreads()])
    end
end

function Base.hash(lb::LazyBranch, h::UInt)
    h = hash(lb.f, h)
    h = hash(lb.b.fClassName, h)
    h = hash(lb.L, h)
    for br in lb.buffer_range
        h = hash(br, h)
    end
    return h
end
Base.size(ba::LazyBranch) = (ba.L,)
Base.length(ba::LazyBranch) = ba.L
Base.firstindex(ba::LazyBranch) = 1
Base.lastindex(ba::LazyBranch) = ba.L
Base.eltype(ba::LazyBranch{T,J,B}) where {T,J,B} = T

function Base.show(io::IO, lb::LazyBranch)
    summary(io, lb)
    println(":")
    println("  File: $(lb.f.filename)")
    println("  Branch: $(lb.b.fName)")
    println("  Description: $(lb.b.fTitle)")
    println("  NumEntry: $(lb.L)")
    print("  Entry Type: $(eltype(lb))")
    nothing
end

"""
    Base.getindex(ba::LazyBranch{T, J}, idx::Integer) where {T, J}

Get the `idx`-th element of a `LazyBranch`, starting at `1`. If `idx` is
within the range of `ba.buffer_range`, it will directly return from `ba.buffer`.
If not within buffer, it will fetch the correct basket by calling [`basketarray`](@ref)
and update buffer and buffer range accordingly.

!!! warning
    Because currently we only cache a single basket inside `LazyBranch` at any given
    moment, access a `LazyBranch` from different threads at the same time can cause
    performance issue and incorrect event result.
"""
function Base.getindex(ba::LazyBranch{T,J,B}, idx::Integer) where {T,J,B}
    tid = Threads.threadid()
    br = ba.buffer_range[tid]
    if idx ∉ br
        seek_idx = findfirst(x -> x > (idx - 1), ba.fEntry) - 1 #support 1.0 syntax
        bb = basketarray(ba.f, ba.b, seek_idx)
        if typeof(bb) !== B
            error("Expected type of interpreted data: $(B), got: $(typeof(bb))")
        end
        ba.buffer[tid] = bb
        br = (ba.fEntry[seek_idx] + 1):(ba.fEntry[seek_idx + 1])
        ba.buffer_range[tid] = br
    end
    localidx = idx - br.start + 1
    return ba.buffer[tid][localidx]
end

function Base.iterate(ba::LazyBranch{T,J,B}, idx=1) where {T,J,B}
    idx > ba.L && return nothing
    return (ba[idx], idx + 1)
end

const _LazyTreeType =
    TypedTables.Table{<:NamedTuple,1,NamedTuple{S,N}} where {S,N<:Tuple{Vararg{LazyBranch}}}

struct LazyTree{T} <: DataFrames.AbstractDataFrame
    treetable::T
    colidx::DataFrames.Index
end
@inline innertable(t::LazyTree) = Core.getfield(t, :treetable)

# a specific branch
Base.getindex(lt::LazyTree, row::Int) = innertable(lt)[row]
function Base.getindex(lt::LazyTree, rang::UnitRange)
    return LazyTree(innertable(lt)[rang], Core.getfield(lt, :colidx))
end
Base.getindex(lt::LazyTree, ::typeof(!), s::Symbol) = lt[:, s]
Base.getindex(lt::LazyTree, ::Colon, i::Int) = lt[:, propertynames(lt)[i]]
Base.getindex(lt::LazyTree, ::typeof(!), i::Int) = lt[:, propertynames(lt)[i]]
Base.getindex(lt::LazyTree, ::Colon, s::Symbol) = getproperty(innertable(lt), s) # the real deal

# a specific event
Base.getindex(lt::LazyTree, row::Int, col::Int) = lt[:, col][row]
Base.getindex(lt::LazyTree, row::Int, col::Symbol) = lt[:, col][row]
Base.getindex(lt::LazyTree, rows::UnitRange, col::Symbol) = lt[:, col][rows]
Base.getindex(lt::LazyTree, ::Colon) = lt[1:end]
Base.firstindex(lt::LazyTree) = 1
Base.lastindex(lt::LazyTree) = length(lt)

# allow enumerate() to be chunkable (eg with Threads.@threads)
Base.firstindex(e::Iterators.Enumerate{LazyTree{T}}) where T = firstindex(e.itr)
Base.lastindex(e::Iterators.Enumerate{LazyTree{T}}) where T = lastindex(e.itr)
Base.getindex(e::Iterators.Enumerate{LazyTree{T}}, row::Int) where T = (row, first(iterate(e.itr, row)))

# interfacing AbstractDataFrame
DataFrames._check_consistency(lt::LazyTree) = nothing #we're read-only
Base.names(lt::LazyTree) = collect(String.(propertynames(innertable(lt))))
DataFrames.index(lt::LazyTree) = Core.getfield(lt, :colidx)
DataFrames.ncol(lt::LazyTree) = length(DataFrames.index(lt))
Base.length(lt::LazyTree) = length(innertable(lt))
DataFrames.nrow(lt::LazyTree) = length(lt)

function getbranchnamesrecursive(obj)
    out = Vector{String}()
    for b in obj.fBranches.elements
        push!(out, b.fName)
        for subname in getbranchnamesrecursive(b)
            push!(out, "$(b.fName)/$(subname)")
        end
    end
    return out
end

"""
    LazyTree(f::ROOTFile, s::AbstractString, branche::Union{AbstractString, Regex})
    LazyTree(f::ROOTFile, s::AbstractString, branches::Vector{Union{AbstractString, Regex}})

Constructor for `LazyTree`, which is close to an `AbstractDataFrame` (interface wise),
and a lazy `TypedTables.Table` (speed wise). Looping over a `LazyTree` is fast and type
stable. Internally, `LazyTree` contains a typed table whose branch are [`LazyBranch`](@ref).
This means that at any given time only `N` baskets are cached, where `N` is the number of branches.

!!! note
    Accessing with `[start:stop]` will return a `LazyTree` with concrete internal table.

# Example
```julia
julia> mytree = LazyTree(f, "Events", ["Electron_dxy", "nMuon", r"Muon_(pt|eta)\$"])
 Row │ Electron_dxy     nMuon   Muon_eta         Muon_pt
     │ Vector{Float32}  UInt32  Vector{Float32}  Vector{Float32}
─────┼───────────────────────────────────────────────────────────
 1   │ [0.000371]       0       []               []
 2   │ [-0.00982]       2       [0.53, 0.229]    [19.9, 15.3]
 3   │ []               0       []               []
 4   │ [-0.00157]       0       []               []
 ⋮   │     ⋮            ⋮             ⋮                ⋮
```
"""
function LazyTree(f::ROOTFile, s::AbstractString, branches)
    tree = f[s]
    tree isa TTree || error("$s is not a tree name.")
    if length(branches) > 30
        @warn "Your tree is quite wide, with $(length(branches)) branches, this will take compiler a moment."
    end
    d = Dict{Symbol,LazyBranch}()
    d_colidx = Dict{Symbol,Int}()
    _m(s::AbstractString) = isequal(s)
    _m(r::Regex) = Base.Fix1(occursin, r)
    branches = mapreduce(b -> filter(_m(b), getbranchnamesrecursive(tree)), ∪, branches)
    SB = Symbol.(branches)
    for (i, b) in enumerate(SB)
        d[b] = f["$s/$b"]
        d_colidx[b] = i
    end
    return LazyTree(TypedTables.Table(d), DataFrames.Index(d_colidx, SB))
end

function LazyTree(f::ROOTFile, s::AbstractString)
    return LazyTree(f, s, keys(f[s]))
end

function LazyTree(f::ROOTFile, s::AbstractString, branch::Union{AbstractString,Regex})
    return LazyTree(f, s, [branch])
end

struct LazyEvent{T<:TypedTables.Table}
    tree::T
    idx::Int64
end
Base.show(io::IO, evt::LazyEvent) = show(io, "LazyEvent with: $(propertynames(evt))")
function Base.getproperty(evt::LazyEvent, s::Symbol)
    @inbounds getproperty(Core.getfield(evt, :tree), s)[Core.getfield(evt, :idx)]
end
Base.collect(evt::LazyEvent) = Core.getfield(evt, :tree)[Core.getfield(evt, :idx)]

function Base.iterate(tree::T, idx=1) where {T<:LazyTree}
    idx > length(tree) && return nothing
    return LazyEvent(Core.getfield(tree, :treetable), idx), idx + 1
end

# TODO this is not terribly slow, but we can get faster implementation still ;)
function Base.getindex(ba::LazyBranch{T,J,B}, rang::UnitRange) where {T,J,B}
    return [ba[i] for i in rang]
end


text:
using Distributions
using Base.Test

mu = [1.0, 2.0, 3.0]
C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]
h = mu
J = C
mv_test_distributions = [
	                 Dirichlet(3, 2.0),
	                 Dirichlet([2.0, 1.0, 3.0]),
	                 IsoNormal(mu, 2.0),
	                 DiagNormal(mu, [1.5, 2.0, 2.5]),
	                 MvNormal(mu, C),
	                 IsoNormalCanon(h, 2.0),
	                 DiagNormalCanon(h, [1.5, 2.0, 1.2]),
	                 MvNormalCanon(h, J)]

for d in mv_test_distributions
    @test size(d) == size(rand(d))
    @test length(d) == length(rand(d))
end


text:
#!/usr/bin/env julia
# Author: Joseph Sahyoun
# This node uses MPC algorithms to calculate the optimal inputs that minimizes
# the distance between the BARC and the reference trajectory. For this demonstration
# a dummy reference trajectory is fed into this node which then publishes the optimal
# state trajectory.

println("Initializing Packages...")

using RobotOS
@rosimport std_msgs.msg: Float64MultiArray, Float64
@rosimport barc.msg: Input, barc_state
@rosimport data_service.msg: TimeData
@rosimport geometry_msgs.msg: Vector3
rostypegen()
using barc.msg
using data_service.msg
using geometry_msgs.msg
using std_msgs.msg
using JuMP
using Ipopt

println("DONE")

# -------------------------------------------------------------------
#                           MPC PARAMETERS
# -------------------------------------------------------------------

println("Defining Parameters...")

# Number of states / Number of inputs
nz = get_param("/numStates")
nu = get_param("/numInputs")

#### MAKE MPC TERMS PARAMETERS
# Horizons
n = get_param("/numStepsToLookAhead")           # MPC Horizon (CFTOC Horizon)

# Cost matrices
Q = get_param("/Q")
Q = reshape(Q,(nz,nz))
R = get_param("/R")
R = reshape(R,(nu,nu))

# Hardware parameters
lr = get_param("/lr")
lf = get_param("/lf")

#Box constraints on states
x_min = get_param("/x_min")
x_max = get_param("/x_max")
y_min = get_param("/y_min")
y_max = get_param("/y_max")

# Box constraints on input, A*u <= B
v_max = get_param("/v_max")
v_min = get_param("/v_min")
steer_min = get_param("/steer_min")
steer_max = get_param("/steer_max")

vref = get_param("/reference_velocity")
numMovingAveragePoints = get_param("/numMovingAveragePoints")
velSum = []
deltaSum = []
count = 0
version = 1

println("DONE")

# -------------------------------------------------------------------
#                         SET UP CFTOC MODEL
# -------------------------------------------------------------------

println("Creating CFTOC Model...")

m = Model(solver = IpoptSolver(print_level=0))

# Variables
@variable(m,z[1:nz,1:n+1])              # State vectors (nz x n+1)
@variable(m,u[1:nu,1:n])                # Input vectors (nu x n)

# Reference Trajectory
@NLparameter(m, xref[i=1:n] == 0)
@NLparameter(m, yref[i=1:n] == 0)
@NLparameter(m, pref[i=1:n] == 0)

# Objective Function

if version == 1
    @NLobjective(m, :Min, sum(Q[1,1]*(z[1,i+1]-xref[i])^2 + Q[2,2]*(z[2,i+1]-yref[i])^2 for i in 1:n) + sum(R[1,1]*(u[1,i]-vref)^2 for i in 1:n))
end

# State Constraints

for i in 1:n+1
    @constraint(m, x_min <= z[1,i] <= x_max)
    @constraint(m, y_min <= z[2,i] <= y_max)
    @constraint(m, -3.1415926/2 <= z[3,i] <= 3.1415926/2)
end

for i in 1:n
    @constraint(m, v_min <= u[1,i] <= v_max)
    @constraint(m, steer_min <= u[2,i] <= steer_max)
end

# Initial Condition and Constraints
@NLparameter(m, x0 == 0); @NLconstraint(m, z[1,1] == x0);
@NLparameter(m, y0 == 0); @NLconstraint(m, z[2,1] == y0);
@NLparameter(m, p0 == 0); @NLconstraint(m, z[3,1] == p0);
@NLparameter(m, dt == 0.05);

# Dynamics Constraints
for i in 1:n
    @NLconstraint(m, z[1,i+1] == z[1,i] + dt*u[1,i]*cos(z[3,i]+atan(lr*tan(u[2,i])/(lf+lr))))
    @NLconstraint(m, z[2,i+1] == z[2,i] + dt*u[1,i]*sin(z[3,i]+atan(lr*tan(u[2,i])/(lf+lr))))
    @NLconstraint(m, z[3,i+1] == z[3,i] + dt*(u[1,i]/lr)*sin(atan(lr*tan(u[2,i])/(lf+lr))))
end

println("DONE")

# -------------------------------------------------------------------
#                             INITIAL SOLVE
# -------------------------------------------------------------------

println("Initial solve...")
solve(m)
println("DONE")

# -------------------------------------------------------------------
#                             ROS FUNCTIONS
# -------------------------------------------------------------------

println("Running MPC...")

function dt_callback(msg::Float64Msg)
    setvalue(dt, msg.data)
end

function callback2(msg::barc_state)
    # Assign new reference trajectory
    for i in 1:n
        setvalue(xref[i], msg.x[i])
        setvalue(yref[i], msg.y[i])
        #setvalue(pref[i], msg.psi[i])
    end
end

function loop(pub_one, pub_two)
    global count,velSum,deltaSum
    while ! is_shutdown()
	loop_rate = Rate(get_param("/loop_rate"))
	count = count+1
        # -----------------------------------------------------------
        #                      RUNNING MPC
        # -----------------------------------------------------------
		if count>1
	    toc()
		end
        tic()
        solve(m)                                # Solve CFTOC
        zOpt = getvalue(z)
        uOpt = getvalue(u)
        JOpt = getobjectivevalue(m)
		if count<=numMovingAveragePoints
			velSum = [velSum;uOpt[1,1]]
			deltaSum = [deltaSum;uOpt[2,1]]
		else
			velSum = [velSum[2:numMovingAveragePoints];uOpt[1,1]]
			deltaSum = [deltaSum[2:numMovingAveragePoints];uOpt[2,1]]
		end
        # Publish optimal inputs and optimal state trajectory
        # Only need to publish state traj for this demonstration

        cmdInt = Input(mean(velSum), mean(deltaSum))    # Command to publish
        optTraj = barc_state(zOpt[1,1:end],zOpt[2,1:end],zOpt[3,1:end])
        publish(pub_one,cmdInt)
        publish(pub_two,optTraj)

        rossleep(loop_rate)

    end
end

function main()
    init_node("MPC_node")
	global count,velSum,deltaSum
    pub = Publisher("uOpt", Input, queue_size=10)
    pub2 = Publisher("optimal_state_trajectory", barc_state, queue_size=1)
    sub1 = Subscriber("dt", Float64Msg, dt_callback, queue_size=1)
    sub2 = Subscriber("reference_trajectory", barc_state, callback2, queue_size=1)
    loop(pub, pub2)
end

if ! isinteractive()
    main()
end

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc

max_forks_repo_name: treyfortmuller/barc


text:
using Test
using NumberIntervals

const a = NumberInterval(-1, 0)
const b = NumberInterval(-0.5, 0.5)
const c = NumberInterval(0.5, 2)
const d = NumberInterval(0.25, 0.8)
const z = zero(NumberInterval)
const e = NumberInterval(Inf, -Inf)

@testset "number comparison" begin
    @test a < c
    @test c > a
    @test (a < b) |> ismissing
    @test (c > b) |> ismissing
    @test !(c < a)
    @test !(a > c)
    @test z == z
    @test z != c
    @test (a == b) |> ismissing
    @test (b != c) |> ismissing
    @test b <= c
end
@testset "testing for zero" begin
    @test !iszero(c)
    @test iszero(z)
    @test iszero(a) |> ismissing
    @test iszero(b) |> ismissing
end
@testset "test sign" begin
    @test signbit(c) == false
    @test signbit(-a) == false
    @test signbit(-c) == true
    @test signbit(a) |> ismissing
    @test sign(c) == 1
    @test sign(z) == 0
    @test sign(-c) == -1
    @test sign(b) |> ismissing
end
@testset "isinteger" begin
    @test isinteger(z)
    @test isinteger(NumberInterval(4))
    @test !isinteger(NumberInterval(4.5))
    @test isinteger(c) |> ismissing
    @test !isinteger(d)
end
@testset "isfinite" begin
    @test isfinite(a)
    @test isfinite(b)
    @test isfinite(c)
    @test isfinite(z)
    @test isfinite(NumberInterval(0., Inf))
    @test isfinite(e) |> ismissing
end
@testset "IndeterminateException" begin
    @test_throws IndeterminateException throw(IndeterminateException())
end
@testset "constructor" begin
    @test NumberInterval(a) === a
    @test NumberInterval{Float64}(a) === a
    @test NumberInterval{Float32}(4) isa NumberInterval{Float32}
    @test real(a) === a
    @test_throws ErrorException NumberInterval(2., 1.)
    @test isnan(NumberInterval(NaN))
    @test_throws ErrorException NumberInterval(Inf)
    @test_throws ErrorException NumberInterval(-Inf)
end
@testset "promotion" begin
    @test promote_rule(NumberInterval{Float32}, Float64) ==
        NumberInterval{Float64}
end


text:
using Test
using MARY_fRG.Basics
using MARY_fRG.Fermi, MARY_fRG.Triangulated
using MARY_fRG.Fermi.Surface, MARY_fRG.Fermi.Patch
using MARY_fRG.FlowEquation

@testset "ec测试结果" begin
    qs = common_square_lattice(0.)
    pnum = 44
    ltris, ladjs = split_square(qs.brillouin, 40)
    lpats = group_ltris_into_patches_mt(ltris, qs.brillouin, pnum)
    sf = const_energy_line(ltris, ladjs, 0., qs.dispersion[1])
    res = pi_αβ_plus_ec(
        sf, sf, 4pi^2, 0.2, Point2D(0., 0.2), qs.dispersion[1], qs.kadd
    )
    @test isapprox(res, 0.116678333187, atol=1e-12)
    res = pi_αβ_minus_ec(
        sf, sf, 4pi^2, 0.2, Point2D(-0.15, 0.), qs.dispersion[1], qs.kadd
    )
    @test isapprox(res, 0.107264645635, atol=1e-12)
    #
    pnum = 8
    lpats = group_ltris_into_patches_mt(ltris, qs.brillouin, pnum)
    patcs = patches_under_vonhove(qs.brillouin, qs.dispersion[1], pnum)
    mpatcs = [patcs]
    #Γ4 = Gamma4(
    #    qs, 8.0, zeros(1, 1, 1, 1, pnum, pnum, pnum),
    #    zeros(Int64, 1, 1, 1, 1, pnum, pnum, pnum),
    #    pnum, mpatcs,
    #    ltris, lpats, ladjs, nothing
    #)
    Γ4 = ECGamma4(qs, 8.0, pnum, 40)
    @time all_bubble_ec_mt(Γ4, 1.0)
end

@testset "tf测试结果" begin
    qs = common_square_lattice(0.)
    pnum = 44
    ltris, ladjs = split_square(qs.brillouin, 40)
    lpats = group_ltris_into_patches_mt(ltris, qs.brillouin, pnum)
    res = pi_αβ_plus_tf(
        ltris, 4pi^2, area(ltris[1]), 0.2, Point2D(0., 0.2),
        qs.dispersion[1], qs.dispersion[1], qs.kadd
    )
    @test isapprox(res, 0.0503171, atol=1e-6)
    res = pi_αβ_minus_tf(
        ltris, 4pi^2, area(ltris[1]), 0.2, Point2D(-0.15, 0.),
        qs.dispersion[1], qs.dispersion[1], qs.kadd
    )
    @test isapprox(res, 0.2063251, atol=1e-6)
    #
    pnum = 8
    lpats = group_ltris_into_patches_mt(ltris, qs.brillouin, pnum)
    patcs = patches_under_vonhove(qs.brillouin, qs.dispersion[1], pnum)
    mpatcs = [patcs]
    #Γ4 = Gamma4(
    #    qs, 8.0, zeros(1, 1, 1, 1, pnum, pnum, pnum),
    #    zeros(Int64, 1, 1, 1, 1, pnum, pnum, pnum),
    #    pnum, mpatcs,
    #    ltris, lpats, ladjs, nothing
    #)
    Γ4 = TFGamma4(qs, 8.0, pnum, 40)
    @time all_bubble_tf_mt(Γ4, 1.0)
    #println(Γ4.ltris_pat)
    @time all_bubble_tf_mt(Γ4, 1.0)
    #println(Γ4.ltris_pat)
end


text:
function mkelemblocks(ne,ns::Int=512)

ns = min(ns,ne);
nb = Int(floor(ne/ns));  # number of blocks
na = Int(round(ne/nb)); # number of elements per block
nk = collect(1:na:ne);
nk = reshape(nk,1,length(nk));
if length(nk)==1
    nm = [1; ne];
else
    tm = 0*nk;
    tm[1:end-1] = nk[2:end].-1;
    tm[end] = ne;
    nm = [nk; tm];
    if (nm[2,end]-nm[1,end]) < (na/2)
        nm[2,end-1] = nm[2,end];
        nm = nm[:,1:(end-1)];
    end
end

nb = size(nm,2);
while minimum(nm[2,:]-nm[1,:])<(ns/2-1) || maximum(nm[2,:]-nm[1,:])>ns
    nb = nb+1;
    na = round(ne/nb); # number of elements per block
    nk = collect(1:na:ne);
    nk = reshape(nk,1,length(nk));
    tm = 0*nk;
    tm[1:end-1] = nk[2:end].-1;
    tm[end] = ne;
    nm = [nk; tm];
    #nm = [nk[1:end]; [nk[2:end]-1 ne]];
    if (nm[2,end]-nm[1,end]) < (na/2)
        nm[2,end-1] = nm[2,end];
        nm = nm[:,1:(end-1)];
    end
end
nm = Int.(nm);
nb = size(nm,2);

if nm[2,end]!=ne
    error("something wrong");
end
if maximum(nm[2,:]-nm[1,:])>ns
    error("something wrong");
end
if minimum(nm[2,:]-nm[1,:])<(ns/2-1)
    error("something wrong");
end

return nm,nb

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
function get_lists_list(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/list.json", options))
end

#Requires list_id and slug
function get_lists_statuses(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/statuses.json", options))
end

function get_lists_memberships(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/memberships.json", options))
end

#Requires list_id and slug
function get_lists_subscribers(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/subscribers.json", options))
end

#Requires list_id and slug
function get_lists_subscribers_show(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/subscribers/show.json", options))
end

#Requires list_id and slug
function get_lists_members_show(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/members/show.json", options))
end

#Requires list_id and slug
function get_lists_members(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/members.json", options))
end

#Requires list_id and slug
function get_lists_show(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/show.json", options))
end

function get_lists_subscriptions(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/subscriptions.json", options))
end

function get_lists_ownerships(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/lists/ownerships.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
# Define abstract interface
"""
$(TYPEDEF)

Abstract supertype for excitation fields.
"""
abstract type ExcitationFields{T,N} end

"""
$(SIGNATURES)

Initialize an excitation field given a `selectionField` and multiple
`excitationFields`.
"""
function ExcitationField(selectionField::MagneticFields{T}, excitationFields::U) where {T,U<:Tuple{Vararg{W} where W<:MagneticFields{T}}}
    error("ExcitationField is not implemented for selectionField::$(typeof(selectionField)) and excitationFields::$(typeof(excitationFields)).")
end

function ExcitationField(selectionField::MagneticFields, excitationFields...)
    return ExcitationField(selectionField,tuple(excitationFields...))
end

"""
$(SIGNATURES)

Return the magnetic field strength H [T/μ₀] (μ₀ being the vacuum permeability)
as `SVector{3,T}` for a given location r [m] and time t [s].
"""
function magneticFieldStrength(excitation::ExcitationFields{T,N},r::SVector{3,T},t::T) where {T,N}
    error("magneticFieldStrength is not implemented for $(typeof(excitation)).")
end

# Define specific subtypes
## Selection field and N-dimensional sinusoidal excitation field
struct SinusoidalExcitationField{T,N,U<:StaticMagneticFields{T}} <: ExcitationFields{T,N}
    selectionField::U
    sinusoidalFields::NTuple{N,SinusoidalField{T}}
end

function ExcitationField(selectionField::U, magneticFields::NTuple{N,SinusoidalField{T}}) where {T,N,U<:StaticMagneticFields{T}}
    return SinusoidalExcitationField{T,N,U}(selectionField,magneticFields)
end

function magneticFieldStrength(excitation::SinusoidalExcitationField{T,N},r::SVector{3,T},t::T) where {T,N}
    H = magneticFieldStrength(excitation.selectionField,r,t)
    @inbounds for i =1:N
        H += magneticFieldStrength(excitation.sinusoidalFields[i],r,t)
    end
    return H
end


text:
using Test
using CuYao
using CuYao: tri2ij
using LinearAlgebra
using BitBasis
using Statistics: mean
using StaticArrays
using CUDA
CUDA.allowscalar(false)

@testset "basics" begin
    a = randn(ComplexF64, 50, 20)
    ca = a|>cu
    batch_normalize!(ca)
    batch_normalize!(a)
    @test ca |> Matrix ≈ a

    for l = 1:100
        i, j = tri2ij(l)
        @test (i-1)*(i-2)÷2+j == l
    end
end

@testset "constructor an measure" begin
    reg = rand_state(10)
    greg = reg |> cu
    @test greg isa GPUReg
    @test eltype(greg.state) == ComplexF64
    myvec(x) = Vector(x)
    myvec(x::Number) = [x]
    for reg in [rand_state(10, nbatch=333), rand_state(10)]
        greg = reg |> cu
        @test size(measure(greg |> copy, nshots=10)) == size(measure(reg, nshots=10))
        @test size(measure!(greg |> copy)) == size(measure!(reg |> copy))
        @test size(measure!(ResetTo(0), greg |> copy)) == size(measure!(ResetTo(0), reg |> copy))
        @test size(measure!(RemoveMeasured(), greg |> copy)) == size(measure!(RemoveMeasured(), reg |> copy))
        @test select(greg |> copy, 12) |> cpu ≈ select(reg, 12)
        @test size(measure!(greg |> copy |> focus!(3,4,1))) == size(measure!(reg |> copy |> focus!(3,4,1)))
        @test greg |> copy |> focus!(3,4,1) |> relax!(3,4,1) |> cpu ≈ reg

        if nbatch(greg) == 1
            greg1 = greg |> copy |> focus!(1,4,3)
            greg0 = copy(greg1)
            res = measure!(RemoveMeasured(), greg1)
            @test select(greg0, res |> myvec) |> normalize! |> cpu ≈ greg1 |> cpu
        end

        greg1 = greg |> copy |> focus!(1,4,3)
        greg0 = copy(greg1)
        res = measure!(ResetTo(3), greg1)
        @test all(measure(greg1, nshots=10) .== 3)
        @test greg1 |> isnormalized
        @test all(select.(greg0 |> cpu, res |> myvec) .|> normalize! .≈ select.(greg1 |> cpu, 3))

        greg1 = greg |> copy |> focus!(1,4,3)
        greg0 = copy(greg1)
        res = measure!(greg1)
        @test all(select.(greg0 |> cpu, res |> myvec) .|> normalize! .≈ select.(greg1 |> cpu, res|>myvec))
    end

    @test join(rand_state(3) |> cu, rand_state(3) |> cu) |> nactive == 6
    @test join(rand_state(3, nbatch=10) |> cu, rand_state(3, nbatch=10) |> cu) |> nactive == 6
end

@testset "insert_qubits!" begin
    reg = rand_state(5; nbatch=10)
    res = insert_qubits!(reg |> cu, 3; nqubits=2) |> cpu
    @test insert_qubits!(reg, 3; nqubits=2) ≈ res

    reg = rand_state(5, nbatch=10) |>focus!(2,3)
    res = insert_qubits!(reg |> cu, 3; nqubits=2) |> cpu
    @test insert_qubits!(reg, 3; nqubits=2) ≈ res
end

@testset "cuda-op-measures" begin
    reg = rand_state(8; nbatch=32) |> cu
    op = repeat(5, X, 1:5)

    # measure!
    reg2 = reg |> copy
    res = measure!(op, reg2, 2:6)
    res2 = measure!(op, reg2, 2:6)
    @test size(res) == (32,)
    @test res2 == res

    # measure! and reset
    reg2 = reg |> copy
    res = measure!(ResetTo(0), op, reg2, 2:6)
    reg2 |> repeat(8, H, 2:6)
    res2 = measure!(ResetTo(0), op, reg2, 2:6)
    @test size(res) == (32,) == size(res2)
    @test all(res2 .== 1)

    # measure! and remove
    reg2 = reg |> copy
    res = measure!(RemoveMeasured(), op, reg2, 2:6)
    @test size(res) == (32,)

    reg = repeat(ArrayReg([1,-1+0im]/sqrt(2.0)), 10) |> cu
    @test measure!(X, reg) |> mean ≈ -1
    reg = repeat(ArrayReg([1.0,0+0im]), 1000)
    @test abs(measure!(X, reg) |> mean) < 0.1
end

@testset "cuda kron getindex" begin
    a = randn(3,4)
    b = randn(4,2)
    c = zeros(12,8)
    ca, cb, cc = cu(a), cu(b), cu(c)
    @test kron(ca, cb) |> Array ≈ kron(a, b)
    @test kron!(cc, ca, cb) |> Array ≈ kron(a,b)

    kron!(c,a,b)
    @test cc |> Array ≈ c

    v = randn(100) |> cu
    inds = [3,5,2,1,7,1]
    @test v[inds] ≈ v[inds |> CuVector]
end

@testset "cuda batched_kron" begin
    a = randn(3,4,5)
    b = randn(4,2,5)
    c = zeros(12,8,5)
    ca, cb, cc = cu(a), cu(b), cu(c)

    @test batched_kron(ca, cb) |> Array ≈ batched_kron(a, b)
    @test batched_kron!(cc, ca, cb) |> Array ≈ batched_kron(a, b)

    batched_kron!(c, a, b)
    @test cc |> Array ≈ c
end


text:
################################################################################
#  Copyright 2021, Tom Van Acker, Frederik Geth                                #
################################################################################
# StochasticPowerModels.jl                                                     #
# An extention package of PowerModels.jl for Stochastic (Optimal) Power Flow   #
# See http://github.com/timmyfaraday/StochasticPowerModels.jl                  #
################################################################################

# variables
""
function variable_branch_current(pm::AbstractIVRModel; nw::Int=nw_id_default, aux::Bool=true, bounded::Bool=true, report::Bool=true, aux_fix::Bool=false, kwargs...)
    _PM.variable_branch_series_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    _PM.variable_branch_series_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    _PM.variable_branch_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    _PM.variable_branch_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    if aux
        variable_branch_series_current_squared(pm, nw=nw, bounded=bounded, report=report, aux_fix=aux_fix; kwargs...)
    else
        if nw == nw_id_default
            variable_branch_series_current_expectation(pm, nw=nw, bounded=bounded, report=report, aux_fix=aux_fix; kwargs...)
            variable_branch_series_current_variance(pm, nw=nw, bounded=bounded, report=report, aux_fix=aux_fix; kwargs...)
    end end
end

""
function variable_branch_current_reduced(pm::AbstractIVRModel; nw::Int=nw_id_default, aux::Bool=true, bounded::Bool=true, report::Bool=true, aux_fix::Bool=false, kwargs...)
    _PM.variable_branch_series_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    _PM.variable_branch_series_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    expression_variable_branch_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    expression_variable_branch_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    if aux
        variable_branch_series_current_squared(pm, nw=nw, bounded=bounded, report=report, aux_fix=aux_fix; kwargs...)
    else
        if nw == nw_id_default
            variable_branch_series_current_expectation(pm, nw=nw, bounded=bounded, report=report, aux_fix=aux_fix; kwargs...)
            variable_branch_series_current_variance(pm, nw=nw, bounded=bounded, report=report, aux_fix=aux_fix; kwargs...)
    end end
end

""
function variable_load_current(pm::AbstractIVRModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true, kwargs...)
    variable_load_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    variable_load_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)
end

""
function variable_gen_current(pm::AbstractIVRModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true, kwargs...)
    _PM.variable_gen_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    _PM.variable_gen_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)
end

# general constraints
""
function constraint_current_balance(pm::AbstractIVRModel, n::Int, i, bus_arcs, bus_gens, bus_loads, bus_gs, bus_bs)
    vr = _PM.var(pm, n, :vr, i)
    vi = _PM.var(pm, n, :vi, i)

    cr = _PM.var(pm, n, :cr)
    ci = _PM.var(pm, n, :ci)

    crd = _PM.var(pm, n, :crd)
    cid = _PM.var(pm, n, :cid)
    crg = _PM.var(pm, n, :crg)
    cig = _PM.var(pm, n, :cig)

    JuMP.@constraint(pm.model,  sum(cr[a] for a in bus_arcs)
                                ==
                                sum(crg[g] for g in bus_gens)
                                - sum(crd[d] for d in bus_loads)
                                - sum(gs for gs in values(bus_gs))*vr + sum(bs for bs in values(bus_bs))*vi
                                )
    JuMP.@constraint(pm.model,  sum(ci[a] for a in bus_arcs)
                                ==
                                sum(cig[g] for g in bus_gens)
                                - sum(cid[d] for d in bus_loads)
                                - sum(gs for gs in values(bus_gs))*vi - sum(bs for bs in values(bus_bs))*vr
                                )
end

# galerkin projection
""
function constraint_gp_branch_series_current_squared(pm::AbstractIVRModel, n::Int, i, T2, T3)
    css  = _PM.var(pm, n, :css, i)
    csr = Dict(nw => _PM.var(pm, nw, :csr, i) for nw in _PM.nw_ids(pm))
    csi = Dict(nw => _PM.var(pm, nw, :csi, i) for nw in _PM.nw_ids(pm))

    JuMP.@constraint(pm.model,  T2.get([n-1,n-1]) * css
                                ==
                                sum(T3.get([n1-1,n2-1,n-1]) *
                                    (csr[n1] * csr[n2] + csi[n1] * csi[n2])
                                    for n1 in _PM.nw_ids(pm), n2 in _PM.nw_ids(pm))
                    )
end

""
function constraint_gp_gen_power_real(pm::AbstractIVRModel, n::Int, i, g, T2, T3)
    vr  = Dict(nw => _PM.var(pm, nw, :vr, i) for nw in _PM.nw_ids(pm))
    vi  = Dict(nw => _PM.var(pm, nw, :vi, i) for nw in _PM.nw_ids(pm))

    crg = Dict(nw => _PM.var(pm, nw, :crg, g) for nw in _PM.nw_ids(pm))
    cig = Dict(nw => _PM.var(pm, nw, :cig, g) for nw in _PM.nw_ids(pm))

    pg  = _PM.var(pm, n, :pg, g)

    JuMP.@constraint(pm.model,  T2.get([n-1,n-1]) * pg
                                ==
                                sum(T3.get([n1-1,n2-1,n-1]) *
                                    (vr[n1] * crg[n2] + vi[n1] * cig[n2])
                                    for n1 in _PM.nw_ids(pm), n2 in _PM.nw_ids(pm))
                    )
end

""
function constraint_gp_gen_power_imaginary(pm::AbstractIVRModel, n::Int, i, g, T2, T3)
    vr  = Dict(nw => _PM.var(pm, nw, :vr, i) for nw in _PM.nw_ids(pm))
    vi  = Dict(nw => _PM.var(pm, nw, :vi, i) for nw in _PM.nw_ids(pm))

    crg = Dict(nw => _PM.var(pm, nw, :crg, g) for nw in _PM.nw_ids(pm))
    cig = Dict(nw => _PM.var(pm, nw, :cig, g) for nw in _PM.nw_ids(pm))

    qg  = _PM.var(pm, n, :qg, g)

    JuMP.@constraint(pm.model,  T2.get([n-1,n-1]) * qg
                                ==
                                sum(T3.get([n1-1,n2-1,n-1]) *
                                    (vi[n1] * crg[n2] - vr[n1] * cig[n2])
                                    for n1 in _PM.nw_ids(pm), n2 in _PM.nw_ids(pm))
                    )
end

""
function constraint_gp_load_power_real(pm::AbstractIVRModel, n::Int, i, l, pd, T2, T3)
    vr  = Dict(nw => _PM.var(pm, nw, :vr, i) for nw in _PM.nw_ids(pm))
    vi  = Dict(nw => _PM.var(pm, nw, :vi, i) for nw in _PM.nw_ids(pm))

    crd = Dict(nw => _PM.var(pm, nw, :crd, l) for nw in _PM.nw_ids(pm))
    cid = Dict(nw => _PM.var(pm, nw, :cid, l) for nw in _PM.nw_ids(pm))

    JuMP.@constraint(pm.model,  T2.get([n-1,n-1]) * pd
                                ==
                                sum(T3.get([n1-1,n2-1,n-1]) *
                                    (vr[n1] * crd[n2] + vi[n1] * cid[n2])
                                    for n1 in _PM.nw_ids(pm), n2 in _PM.nw_ids(pm))
                    )
end

""
function constraint_gp_load_power_imaginary(pm::AbstractIVRModel, n::Int, i, l, qd, T2, T3)
    vr  = Dict(n => _PM.var(pm, n, :vr, i) for n in _PM.nw_ids(pm))
    vi  = Dict(n => _PM.var(pm, n, :vi, i) for n in _PM.nw_ids(pm))

    crd = Dict(n => _PM.var(pm, n, :crd, l) for n in _PM.nw_ids(pm))
    cid = Dict(n => _PM.var(pm, n, :cid, l) for n in _PM.nw_ids(pm))

    JuMP.@constraint(pm.model,  T2.get([n-1,n-1]) * qd
                                ==
                                sum(T3.get([n1-1,n2-1,n-1]) *
                                    (vi[n1] * crd[n2] - vr[n1] * cid[n2])
                                    for n1 in _PM.nw_ids(pm), n2 in _PM.nw_ids(pm))
                    )
end

# chance constraints
""
function constraint_branch_series_current_cc_limit(pm::AbstractIVRModel, b, cmax, λmax, T2, T4, gs, bs)
    ntws = _PM.nw_ids(pm)

    cse  = _PM.var(pm, nw_id_default, :cse, b)
    csv  = _PM.var(pm, nw_id_default, :csv, b)

    csr  = Dict(n => _PM.var(pm, n, :csr, b) for n in ntws)
    csi  = Dict(n => _PM.var(pm, n, :csi, b) for n in ntws)

    T44 = Dict((n1,n2,n3,n4) => T4.get([n1-1,n2-1,n3-1,n4-1]) for n1 in ntws, n2 in ntws, n3 in ntws, n4 in ntws)

    # expectation
    JuMP.@constraint(pm.model, cse == sum((csr[n]^2 + csi[n]^2) * T2.get([n-1,n-1]) for n in ntws))
    # 'variance'
    JuMP.@NLconstraint(pm.model, cse^2 + csv^2
                                 ==
                                 sum(
                                    (csr[n1] * csr[n2] * csr[n3] * csr[n4] +
                                     2 * csr[n1] * csr[n2] * csi[n3] * csi[n4] +
                                     csi[n1] * csi[n2] * csi[n3] * csi[n4]
                                    ) *
                                    T44[(n1,n2,n3,n4)]
                                    for n1 in ntws, n2 in ntws, n3 in ntws, n4 in ntws
                                 )
                    )
    # chance constraint bounds
    JuMP.@constraint(pm.model,  cse + λmax * csv
                                <=
                                cmax^2
                    )
end

""
function sol_data_model!(pm::AbstractIVRModel, solution::Dict)
    _PM.apply_pm!(_sol_data_model_ivr!, solution)
end

""
function _sol_data_model_ivr!(solution::Dict)
    if haskey(solution, "bus")
        for (i, bus) in solution["bus"]
            if haskey(bus, "vr") && haskey(bus, "vi")
                bus["vm"] = hypot(bus["vr"], bus["vi"])
                bus["va"] = atan(bus["vi"], bus["vr"])
            end
        end
    end
end


text:
"""
Kalman filters, log-likelihood and reasonable priors for a "Celerite"
process, as defined in [Foreman-Mackey, et
al. (2017)[https://arxiv.org/abs/1703.09710].

A Celerite process is a sum of AR(1) and CARMA(2,0) processes; an
alternative description is as a sum of damped random walks and SHOs
driven by white noise.  It is, in fact, a re-parameterization of a
CARMA process, but a convenient one for stably finding mode
frequencies and damping rates and correlation timescales because it is
a sum of low-order terms.

"""

module Celerite

using Ensemble

mutable struct CeleriteKalmanFilter
    mu::Float64 # Mean
    x::Array{ComplexF64, 1} # State mean
    Vx::Array{ComplexF64, 2} # State variance
    K::Array{ComplexF64, 1} # Kalman Gain
    lambda::Array{ComplexF64, 1} # Evolution factors exp(roots*dt)
    b::Array{ComplexF64, 2} # Rotated observation vector
    roots::Array{ComplexF64, 1} # Eigenvalues of the ODEs
    V::Array{ComplexF64, 2} # Stationary covariance
    Vtemp::Array{ComplexF64, 2} # Storage for matrix ops
    drw_rms::Array{Float64, 1}
    drw_rates::Array{Float64, 1}
    osc_rms::Array{Float64, 1}
    osc_freqs::Array{Float64, 1}
    osc_Qs::Array{Float64, 1}
end

function reset!(filt::CeleriteKalmanFilter)
    p = size(filt.x, 1)
    filt.x = zeros(ComplexF64, p)
    filt.Vx = copy(filt.V)
    filt.K = zeros(ComplexF64, p)
    filt
end

function osc_roots(osc_freqs::Array{Float64, 1}, osc_Qs::Array{Float64, 1})
    nosc = size(osc_freqs, 1)

    oscroots = zeros(ComplexF64, nosc)
    for i in 1:nosc
        omega0 = 2.0*pi*osc_freqs[i]
        alpha = omega0/(2.0*osc_Qs[i])
        oscroots[i] = -alpha + omega0*1im
    end
    oscroots
end

function CeleriteKalmanFilter(mu::Float64, drw_rms::Array{Float64, 1}, drw_rates::Array{Float64, 1}, osc_rms::Array{Float64, 1}, osc_freqs::Array{Float64, 1}, osc_Qs::Array{Float64, 1})
    ndrw = size(drw_rms,1)
    nosc = size(osc_rms,1)
    dim = ndrw + 2*nosc

    oscroots = osc_roots(osc_freqs, osc_Qs)

    roots = zeros(ComplexF64, dim)
    ii = 1
    for i in 1:ndrw
        roots[ii] = -drw_rates[i]
        ii+=1
    end

    for i in 1:nosc
        roots[ii] = oscroots[i]
        roots[ii+1] = conj(oscroots[i])
        ii += 2
    end

    b = zeros(ComplexF64, dim)

    iobs = 1
    for i in 1:ndrw
        b[iobs] = 1.0
        iobs += 1
    end

    for i in 1:nosc
        b[iobs] = 1.0
        iobs += 2 # skip the derivative term
    end

    b = b'

    U = zeros(ComplexF64, (dim,dim))

    ii = 1
    for i in 1:ndrw
        U[ii,ii] = 1.0
        ii += 1
    end

    for i in 1:nosc
        U[ii,ii] = 1.0
        U[ii, ii+1] = 1.0
        U[ii+1,ii] = oscroots[i]
        U[ii+1,ii+1] = conj(oscroots[i])
        ii += 2
    end

    b = b*U # Rotated observation vector

    e = zeros(ComplexF64, dim)
    ii = 1
    for i in 1:ndrw
        e[ii] = 1.0
        ii += 1
    end

    for i in 1:nosc
        e[ii+1] = 1.0
        ii += 2
    end

    J = U \ e

    V = zeros(ComplexF64, (dim,dim))

    ii = 1
    for i in 1:ndrw
        V[ii,ii] = -J[ii]*conj(J[ii])/(roots[ii] + conj(roots[ii]))
        ii += 1
    end

    for i in 1:nosc
        for j in 0:1
            for k in 0:1
                V[ii+j, ii+k] = -J[ii+j]*conj(J[ii+k])/(roots[ii+j] + conj(roots[ii+k]))
            end
        end
        ii += 2
    end

    ii = 1
    for i in 1:ndrw
        s2 = b[1,ii]*V[ii,ii]*conj(b[1,ii])
        s2 = s2[1]
        V[ii,ii] *= drw_rms[i]*drw_rms[i]/s2
        ii += 1
    end

    for i in 1:nosc
        s2 = b[:,ii:ii+1]*V[ii:ii+1,ii:ii+1]*b[:,ii:ii+1]'
        s2 = s2[1]
        V[ii:ii+1, ii:ii+1] *= osc_rms[i]*osc_rms[i]/s2
        ii += 2
    end

    CeleriteKalmanFilter(mu, zeros(ComplexF64, dim), V, zeros(ComplexF64, dim), zeros(ComplexF64, dim), b, roots, copy(V), zeros(ComplexF64, (dim,dim)), drw_rms, drw_rates, osc_rms, osc_freqs, osc_Qs)
end

@inbounds function advance!(filt::CeleriteKalmanFilter, dt::Float64)
    p = size(filt.x, 1)

    for i in 1:p
        x = filt.roots[i]*dt
        filt.lambda[i] = exp(x)
    end
    lam = filt.lambda

    for i in 1:p
        filt.x[i] = lam[i]*filt.x[i]
    end

    for j in 1:p
        for i in 1:p
            a::ComplexF64 = lam[i]*conj(lam[j])
            b::ComplexF64 = a*(filt.Vx[i,j] - filt.V[i,j])
            filt.Vx[i,j] = b + filt.V[i,j]
        end
    end

    filt
end

@inbounds function observe!(filt::CeleriteKalmanFilter, y::Float64, dy::Float64)
    p = size(filt.x, 1)

    ey, vy = predict(filt)
    vy += dy*dy

    for i in 1:p
        filt.K[i] = zero(filt.K[i])
        for j in 1:p
            filt.K[i] += filt.Vx[i,j]*conj(filt.b[1,j])/vy
        end
    end

    for i in 1:p
        filt.x[i] = filt.x[i] + (y - ey)*filt.K[i]
    end

    for j in 1:p
        for i in 1:p
            a::ComplexF64 = vy*filt.K[i]
            b::ComplexF64 = a*conj(filt.K[j])
            filt.Vx[i,j] = filt.Vx[i,j] - b
        end
    end

    filt
end

@inbounds function predict(filt::CeleriteKalmanFilter)
    p = size(filt.x,1)

    yp = filt.mu
    for i in 1:p
        yp += real(filt.b[1,i]*filt.x[i])
    end

    vyp = 0.0
    for i in 1:p
        for j in 1:p
            a::ComplexF64 = filt.b[1,i]*filt.Vx[i,j]
            b::ComplexF64 = a*conj(filt.b[1,j])
            vyp = vyp + real(b)
        end
    end

    yp, vyp
end

function whiten(filt::CeleriteKalmanFilter, ts, ys, dys)
    n = size(ts, 1)

    reset!(filt)
    zs = zeros(n)

    for i in 1:n
        yp, vyp = predict(filt)

        zs[i] = (ys[i] - yp)/sqrt(vyp + dys[i]*dys[i])

        observe!(filt, ys[i], dys[i])

        if i < n
            advance!(filt, ts[i+1]-ts[i])
        end
    end

    zs
end

function draw_and_collapse!(filt::CeleriteKalmanFilter)
    nd = size(filt.x, 1)
    try
        for i in 1:nd
            filt.Vx[i,i] = real(filt.Vx[i,i]) # Fix a roundoff error problem?
        end
        L = ctranspose(chol(Hermitian(filt.Vx)))
        filt.x = filt.x + L*randn(nd)
        filt.Vx = zeros(ComplexF64, (nd, nd))
    catch e
        if isa(e, Base.LinearAlgebra.PosDefException)
            warn("Current variance matrix not pos. def.---may be roundoff problem in generation.")
            F = eigfact(filt.Vx)
            for i in eachindex(F[:values])
                l = real(F[:values][i])
                v = F[:vectors][:,i]

                if l < 0.0
                    l = 0.0
                end

                filt.x = filt.x + sqrt(l)*randn()*v
            end
            filt.Vx = zeros(ComplexF64, (nd, nd))
        else
            rethrow()
        end
    end
end

function generate(filt::CeleriteKalmanFilter, ts, dys)
    n = size(ts, 1)
    nd = size(filt.x, 1)

    ys = zeros(n)

    reset!(filt)

    for i in 1:n
        # Draw a new state
        draw_and_collapse!(filt)

        y, _ = predict(filt)

        ys[i] = y + dys[i]*randn()

        if i < n
            advance!(filt, ts[i+1]-ts[i])
        end
    end

    ys
end

function log_likelihood(filt, ts, ys, dys)
    n = size(ts, 1)

    ll = -0.5*n*log(2.0*pi)

    reset!(filt)
    for i in 1:n
        yp, vyp = predict(filt)

        if vyp < 0.0
            warn("Kalman filter has gone unstable!")
            return -Inf
        end

        dy = ys[i] - yp
        vy = vyp + dys[i]*dys[i]

        ll -= 0.5*log(vy)
        ll -= 0.5*dy*dy/vy

        observe!(filt, ys[i], dys[i])

        if i < n
            advance!(filt, ts[i+1]-ts[i])
        end
    end

    ll
end

function residuals(filt, ts, ys, dys)
    n = size(ts, 1)

    resid = zeros(n)
    dresid = zeros(n)

    reset!(filt)
    for i in 1:n
        yp, vyp = predict(filt)

        if vyp < 0.0
            warn("Kalman filter has gone unstable!")
            return resid, dresid
        end

        resid[i] = ys[i] - yp
        dresid[i] = sqrt(vyp + dys[i]*dys[i])

        observe!(filt, ys[i], dys[i])

        if i < n
            advance!(filt, ts[i+1]-ts[i])
        end
    end

    resid, dresid
end

function raw_covariance(ts::Array{Float64, 1}, dys::Array{Float64, 1}, drw_rms::Array{Float64, 1}, drw_rates::Array{Float64, 1}, osc_rms::Array{Float64, 1}, osc_freqs::Array{Float64, 1}, osc_Qs::Array{Float64,1})
    N = size(ts, 1)
    ndrw = size(drw_rms,1)
    nosc = size(osc_rms,1)

    cov = zeros((N,N))
    dts = zeros((N,N))

    oscroots = osc_roots(osc_freqs, osc_Qs)

    for j in 1:N
        for i in 1:N
            dts[i,j] = abs(ts[i] - ts[j])
        end
    end

    for i in 1:ndrw
        cov += drw_rms[i]*drw_rms[i]*exp.(-drw_rates[i].*dts)
    end

    for i in 1:nosc
        A = 1.0 / (-4.0*real(oscroots[i])*(conj(oscroots[i]) - oscroots[i])*oscroots[i])
        B = 1.0 / (-4.0*real(oscroots[i])*(oscroots[i] - conj(oscroots[i]))*conj(oscroots[i]))

        s2 = osc_rms[i]*osc_rms[i] / (A+B)
        cov = cov .+ real.(s2.*(A.*exp.(oscroots[i].*dts) .+ B.*exp.(conj(oscroots[i]).*dts)))
    end

    for i in 1:N
        cov[i,i] += dys[i]*dys[i]
    end

    cov
end

function psd_drw(rms_amp, damp_rate, fs)
    4.0*damp_rate*rms_amp*rms_amp./abs2.(2.0*pi*1im.*fs .+ damp_rate)
end

function psd_osc(rms_amp, freq, Q, fs)
    r1 = osc_roots([freq], [Q])[1]
    r2 = conj(r1)
    norm = 1.0/real(2.0*r1*(r1-r2)*(r1+r2))

    rms_amp*rms_amp/norm./abs2.((2.0*pi*1im.*fs .- r1).*(2.0*pi*1im.*fs .- r2))
end

function psd(filt::CeleriteKalmanFilter, fs::Array{Float64, 1})
    Pfs = zeros(size(fs,1))

    for i in 1:size(filt.drw_rms,1)
        Pfs .= Pfs .+ psd_drw(filt.drw_rms[i], filt.drw_rates[i], fs)
    end

    for i in 1:size(filt.osc_rms,1)
        Pfs .= Pfs .+ psd_osc(filt.osc_rms[i], filt.osc_freqs[i], filt.osc_Qs[i], fs)
    end

    Pfs
end

function predict(filt::CeleriteKalmanFilter, ts, ys, dys, tsp)
    # The algorithm here is to run the filter twice over the union of
    # the data and prediction times.  The first run is "forward" in
    # time, and predicts at each time *before* observing the output at
    # data times, and the second run is "backward" in time, and
    # predicts at each time *after* observing the output at data
    # times.  Weighted-averaging of the two runs together gives the
    # full internal state prediction incorporating all the data at the
    # given times.

    allts = vcat(ts, tsp)
    obsflag = convert(Array{Bool, 1}, vcat(trues(size(ts, 1)), falses(size(tsp, 1))))
    inds = sortperm(allts)
    rinds = reverse(inds)

    yspforward = zeros(size(allts, 1))
    vyspforward = zeros(size(allts, 1))

    reset!(filt)
    for i in eachindex(allts)
        yspforward[inds[i]], vyspforward[inds[i]] = predict(filt)

        if obsflag[inds[i]]
            observe!(filt, ys[inds[i]], dys[inds[i]])
        end

        if i < size(allts, 1)
            advance!(filt, allts[inds[i+1]]-allts[inds[i]])
        end
    end

    yspbackward = zeros(size(allts, 1))
    vyspbackward = zeros(size(allts, 1))

    reset!(filt)
    for i in eachindex(allts)
        if obsflag[rinds[i]]
            observe!(filt, ys[rinds[i]], dys[rinds[i]])
        end

        yspbackward[rinds[i]], vyspbackward[rinds[i]] = predict(filt)

        if i < size(allts, 1)
            advance!(filt, allts[rinds[i]] - allts[rinds[i+1]])
        end
    end

    ysp = (yspforward.*vyspbackward .+ yspbackward.*vyspforward) ./ (vyspbackward .+ vyspforward)
    vysp = 1.0./(1.0./vyspforward .+ 1.0./vyspbackward)

    ysp[size(ts,1)+1:end], vysp[size(ts, 1)+1:end]
end

struct CeleriteKalmanPosterior
    ts::Array{Float64, 1}
    ys::Array{Float64, 1}
    dys::Array{Float64, 1}
    ndrw::Int
    nosc::Int
    mumin::Float64
    mumax::Float64
    fmin::Float64
    fmax::Float64
    ratemin::Float64
    ratemax::Float64
    Qmin::Float64
    Qmax::Float64
    drw_rms_min::Float64
    drw_rms_max::Float64
    osc_rms_min::Float64
    osc_rms_max::Float64
end

struct CeleriteKalmanParams
    mu::Float64
    nu::Float64
    drw_rms::Array{Float64, 1}
    drw_rates::Array{Float64, 1}
    osc_rms::Array{Float64, 1}
    osc_freqs::Array{Float64, 1}
    osc_Qs::Array{Float64, 1}
end

function ndim(post::CeleriteKalmanPosterior)
    2 + 2*post.ndrw + 3*post.nosc
end

function to_params(post::CeleriteKalmanPosterior, x::Array{Float64, 1})
    nd = post.ndrw
    no = post.nosc

    mu = Parameterizations.bounded_value(x[1], post.mumin, post.mumax)
    nu = Parameterizations.bounded_value(x[2], 0.1, 10.0)
    drw_rms = Parameterizations.bounded_value.(x[3:3+nd-1], post.drw_rms_min, post.drw_rms_max)
    drw_rates = Parameterizations.bounded_value.(x[3+nd:3+2*nd-1], post.ratemin, post.ratemax)
    osc_rms = Parameterizations.bounded_value.(x[3+2*nd:3+2*nd+no-1], post.osc_rms_min, post.osc_rms_max)
    osc_freqs = Parameterizations.bounded_value.(x[3+2*nd+no:3+2*nd+2*no-1], post.fmin, post.fmax)
    osc_Qs = Parameterizations.bounded_value.(x[3+2*nd+2*no:3+2*nd+3*no-1], post.Qmin, post.Qmax)

    CeleriteKalmanParams(mu, nu, drw_rms, drw_rates, osc_rms, osc_freqs, osc_Qs)
end

function to_array(post::CeleriteKalmanPosterior, p::CeleriteKalmanParams)
    x = zeros(ndim(post))

    nd = post.ndrw
    no = post.nosc

    x[1] = Parameterizations.bounded_param(p.mu, post.mumin, post.mumax)
    x[2] = Parameterizations.bounded_param(p.nu, 0.1, 10)
    x[3:3+nd-1] = Parameterizations.bounded_param.(p.drw_rms, post.drw_rms_min, post.drw_rms_max)
    x[3+nd:3+2*nd-1] = Parameterizations.bounded_param.(p.drw_rates, post.ratemin, post.ratemax)
    x[3+2*nd:3+2*nd+no-1] = Parameterizations.bounded_param.(p.osc_rms, post.osc_rms_min, post.osc_rms_max)
    x[3+2*nd+no:3+2*nd+2*no-1] = Parameterizations.bounded_param.(p.osc_freqs, post.fmin, post.fmax)
    x[3+2*nd+2*no:3+2*nd+3*no-1] = Parameterizations.bounded_param.(p.osc_Qs, post.Qmin, post.Qmax)

    x
end

function log_prior(post::CeleriteKalmanPosterior, x::Array{Float64, 1})
    log_prior(post, to_params(post, x), x)
end

function log_prior(post::CeleriteKalmanPosterior, p::CeleriteKalmanParams)
    log_prior(post, p, to_array(post, p))
end

function log_prior(post::CeleriteKalmanPosterior, p::CeleriteKalmanParams, x::Array{Float64, 1})
    nd = post.ndrw
    no = post.nosc

    lp = 0.0

    # Flat prior on mu
    lp += Parameterizations.bounded_logjac(p.mu, x[1], post.mumin, post.mumax)

    # Flat-in-log prior on nu
    lp += -log(p.nu) + Parameterizations.bounded_logjac(p.nu, x[2], 0.1, 10.0)

    # Flat in log prior on rms
    lp += -sum(log.(p.drw_rms)) + sum(Parameterizations.bounded_logjac(p.drw_rms, x[3:3+nd-1], post.drw_rms_min, post.drw_rms_max))

    # Flat in log prior on rates
    lp += -sum(log.(p.drw_rates)) + sum(Parameterizations.bounded_logjac(p.drw_rates, x[3+nd:3+2*nd-1], post.ratemin, post.ratemax))

    # Flat in log prior on osc rms
    lp += -sum(log.(p.osc_rms)) + sum(Parameterizations.bounded_logjac(p.osc_rms, x[3+2*nd:3+2*nd+no-1], post.osc_rms_min, post.osc_rms_max))

    # Flat in log prior on osc freqs
    lp += -sum(log.(p.osc_freqs)) + sum(Parameterizations.bounded_logjac(p.osc_freqs, x[3+2*nd+no:3+2*nd+2*no-1], post.fmin, post.fmax))

    # Flat in log prior on Q
    lp += -sum(log.(p.osc_Qs)) + sum(Parameterizations.bounded_logjac(p.osc_Qs, x[3+2*nd+2*no:3+2*nd+3*no-1], post.Qmin, post.Qmax))

    lp
end

function log_likelihood(post::CeleriteKalmanPosterior, x::Array{Float64, 1})
    log_likelihood(post, to_params(post, x))
end

function log_likelihood(post::CeleriteKalmanPosterior, p::CeleriteKalmanParams)
    filt = CeleriteKalmanFilter(p.mu, p.drw_rms, p.drw_rates, p.osc_rms, p.osc_freqs, p.osc_Qs)

    log_likelihood(filt, post.ts, post.ys, p.nu*post.dys)
end

function init(post::CeleriteKalmanPosterior, n::Int)
    nd = post.ndrw
    no = post.nosc

    xs = zeros(ndim(post), n)

    for i in 1:n
        mu = post.mumin + (post.mumax-post.mumin)*rand()
        nu = exp(log(0.1) + (log(10.0) - log(0.1))*rand())

        drw_rms = exp.(log(post.drw_rms_min) + log(post.drw_rms_max/post.drw_rms_min)*rand(nd))
        drw_rates = exp.(log(post.ratemin) + log(post.ratemax/post.ratemin)*rand(nd))

        osc_rms = exp.(log(post.osc_rms_min) + log(post.rms_osc_max/post.osc_rms_min)*rand(no))
        osc_freqs = exp.(log(post.fmin) + log(post.fmax/post.fmin)*rand(no))
        osc_Qs = exp.(log(post.Qmin) + log(post.Qmax/post.Qmin)*rand(no))

        p = CeleriteKalmanParams(mu, nu, drw_rms, drw_rates, osc_rms, osc_freqs, osc_Qs)

        xs[:,i] = to_array(post, p)
    end

    xs
end

end


text:
module IndexSet

# WE USE MATPOWER'S DATA FORMAT

function idx_bus()
    BUS_I       = 1    # bus number (1 to 29997)
    BUS_TYPE    = 2    # bus type (1 - PQ bus, 2 - PV bus, 3 - reference bus, 4 - isolated bus)
    PD          = 3    # Pd, real power demand (MW)
    QD          = 4    # Qd, reactive power demand (MVAr)
    GS          = 5    # Gs, shunt conductance (MW at V = 1.0 p.u.)
    BS          = 6    # Bs, shunt susceptance (MVAr at V = 1.0 p.u.)
    BUS_AREA    = 7    # area number, 1-100
    VM          = 8    # Vm, voltage magnitude (p.u.)
    VA          = 9    # Va, voltage angle (degrees)
    BASE_KV     = 10   # baseKV, base voltage (kV)
    ZONE        = 11   # zone, loss zone (1-999)
    VMAX        = 12   # maxVm, maximum voltage magnitude (p.u.)      (not in PTI format)
    VMIN        = 13   # minVm, minimum voltage magnitude (p.u.)      (not in PTI format)

    LAM_P       = 14   # Lagrange multiplier on real power mismatch (u/MW)
    LAM_Q       = 15   # Lagrange multiplier on reactive power mismatch (u/MVAr)
    MU_VMAX     = 16   # Kuhn-Tucker multiplier on upper voltage limit (u/p.u.)
    MU_VMIN     = 17   # Kuhn-Tucker multiplier on lower voltage limit (u/p.u.)

    return BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, VA, BASE_KV, ZONE, VMAX, VMIN,
    LAM_P, LAM_Q, MU_VMAX, MU_VMIN
end

function idx_branch()
    F_BUS       = 1    # f, from bus number
    T_BUS       = 2    # t, to bus number
    BR_R        = 3    # r, resistance (p.u.)
    BR_X        = 4    # x, reactance (p.u.)
    BR_B        = 5    # b, total line charging susceptance (p.u.)
    RATE_A      = 6    # rateA, MVA rating A (long term rating)
    RATE_B      = 7    # rateB, MVA rating B (short term rating)
    RATE_C      = 8    # rateC, MVA rating C (emergency rating)
    TAP         = 9    # ratio, transformer off nominal turns ratio
    SHIFT       = 10   # angle, transformer phase shift angle (degrees)
    BR_STATUS   = 11   # initial branch status, 1 - in service, 0 - out of service
    ANGMIN      = 12   # minimum angle difference, angle(Vf) - angle(Vt) (degrees)
    ANGMAX      = 13   # maximum angle difference, angle(Vf) - angle(Vt) (degrees)

    PF          = 14   # real power injected at "from" bus end (MW)       (not in PTI format)
    QF          = 15   # reactive power injected at "from" bus end (MVAr) (not in PTI format)
    PT          = 16   # real power injected at "to" bus end (MW)         (not in PTI format)
    QT          = 17   # reactive power injected at "to" bus end (MVAr)   (not in PTI format)

    MU_SF       = 18   # Kuhn-Tucker multiplier on MVA limit at "from" bus (u/MVA)
    MU_ST       = 19   # Kuhn-Tucker multiplier on MVA limit at "to" bus (u/MVA)
    MU_ANGMIN   = 20   # Kuhn-Tucker multiplier lower angle difference limit (u/degree)
    MU_ANGMAX   = 21   # Kuhn-Tucker multiplier upper angle difference limit (u/degree)

    return F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, TAP, SHIFT, BR_STATUS,
    ANGMIN, ANGMAX, PF, QF, PT, QT, MU_SF, MU_ST, MU_ANGMIN, MU_ANGMAX
end

function idx_gen()
    GEN_BUS     = 1    # bus number
    PG          = 2    # Pg, real power output (MW)
    QG          = 3    # Qg, reactive power output (MVAr)
    QMAX        = 4    # Qmax, maximum reactive power output at Pmin (MVAr)
    QMIN        = 5    # Qmin, minimum reactive power output at Pmin (MVAr)
    VG          = 6    # Vg, voltage magnitude setpoint (p.u.)
    MBASE       = 7    # mBase, total MVA base of this machine, defaults to baseMVA
    GEN_STATUS  = 8    # status, 1 - machine in service, 0 - machine out of service
    PMAX        = 9    # Pmax, maximum real power output (MW)
    PMIN        = 10   # Pmin, minimum real power output (MW)
    PC1         = 11   # Pc1, lower real power output of PQ capability curve (MW)
    PC2         = 12   # Pc2, upper real power output of PQ capability curve (MW)
    QC1MIN      = 13   # Qc1min, minimum reactive power output at Pc1 (MVAr)
    QC1MAX      = 14   # Qc1max, maximum reactive power output at Pc1 (MVAr)
    QC2MIN      = 15   # Qc2min, minimum reactive power output at Pc2 (MVAr)
    QC2MAX      = 16   # Qc2max, maximum reactive power output at Pc2 (MVAr)
    RAMP_AGC    = 17   # ramp rate for load following/AGC (MW/min)
    RAMP_10     = 18   # ramp rate for 10 minute reserves (MW)
    RAMP_30     = 19   # ramp rate for 30 minute reserves (MW)
    RAMP_Q      = 20   # ramp rate for reactive power (2 sec timescale) (MVAr/min)
    APF         = 21   # area participation factor

    MU_PMAX     = 22   # Kuhn-Tucker multiplier on upper Pg limit (u/MW)
    MU_PMIN     = 23   # Kuhn-Tucker multiplier on lower Pg limit (u/MW)
    MU_QMAX     = 24   # Kuhn-Tucker multiplier on upper Qg limit (u/MVAr)
    MU_QMIN     = 25   # Kuhn-Tucker multiplier on lower Qg limit (u/MVAr)

    return GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, PC1, PC2, QC1MIN,
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF, MU_PMAX, MU_PMIN, MU_QMAX,
    MU_QMIN
end

function idx_cost()
    MODEL       = 1    # cost model, 1 = piecewise linear, 2 = polynomial
    STARTUP     = 2    # startup cost in US dollars
    SHUTDOWN    = 3    # shutdown cost in US dollars
    NCOST       = 4    # number breakpoints in piecewise linear cost function,
                       # or number of coefficients in polynomial cost function
    COST        = 5    # parameters defining total cost function begin in this col

    return MODEL, STARTUP, SHUTDOWN, NCOST, COST
end

end


text:
"""
    read_dlm_line_after(filename, startpos = 0) -> (Vector{String}(), lastpos)
Read vector of strings from commented text file.
The result is from the first encountered line which doesn't start with '#'.
Also returns the last read position.
"""
function read_dlm_line_after(filename; startpos = 0)
    strconfig, lastpos = open(filename, read = true) do f
        seek(f, startpos)
        st = ""
        if !eof(f)
            st = readline(f)
            while !eof(f) && startswith(st, "#") || st == ""
                st = readline(f)
            end
        end
        st, position(f)
    end
    strip.(split(strconfig, " ", keepempty = false)), lastpos
end

"""
    specbyname(configfi) -> Dict{String, NamedTuple}

Make a channel name => [bytepos, bitpad, bitlength] dictionary.
Argument 'configfi' points to a user configuration file.
"""
function specbyname(configfi)
    # Names of channels
    chnam, lp = read_dlm_line_after(configfi)
    # address field values
    l2, lp = read_dlm_line_after(configfi, startpos = lp)
    l3, lp = read_dlm_line_after(configfi, startpos = lp)
    l4, lp = read_dlm_line_after(configfi, startpos = lp)
    v2 = parse.(Int64, l2)
    v3 = parse.(Int64, l3)
    v4 = parse.(Int64, l4)
    pairs = map(chnam, v2, v3, v4) do nam, va2, va3, va4
        (string(nam) => (; :bytepos => va2, :bitpad => va3, :bitlength => va4))
    end
    Dict{String, NamedTuple}(pairs)
end

"Read DevState given filename"
function devstate(filename)
    st, tist = open(filename, read = true) do f
        st = ""
        tist = ""
        if !eof(f)
            st = readline(f)
            while !eof(f) && startswith(st, "#") || st == ""
                st = readline(f)
            end
        end
        tist = ""
        if !eof(f)
            tist = readline(f)
            while !eof(f) && startswith(tist, "#") || tist == ""
                tist = readline(f)
            end
        end
        st, tist
    end
    if st == "" || tist == ""
        return DevState([], 0, false)
    end
    strstatevec =  split(st, " ", keepempty = false )
    statevec = parse.(UInt8, strstatevec)
    timestamp = parse(Float64, split(tist, " ", keepempty = false)[1])
    DevState(statevec, timestamp, true)
end


text:
#TODO: Fix offset in visualization for periodic bounds
#TODO: Fix visualization for different bounds -- add simple bound

using FFTW
using Plots
using Plots.PlotMeasures
using LinearAlgebra

import Plots: default
default(show=false, reuse=true)

ENV["GKSwstype"]="nul"

abstract type Bounds end

abstract type Periodic <: Bounds end

function signal_plot(i, signal1, signal2, out; filebase = "out",
                     out_name = "convolution", scale = 1.0,
                     bounds_type = :periodic)

    println(i)
    n = length(out)
    temp = -ones(n)

    signal1 = vcat(signal1, zeros(n - length(signal1)))
    signal2 = reverse(signal2)

    if bounds_type == :periodic
        if i < length(signal2)
            shift_signal = real([signal2[end-i+1:end]; zeros(n-length(signal2));
                                signal2[1:end-i]])
        else
            shift_signal = real([zeros(i-length(signal2)); signal2; zeros(n-i)])
        end
    else
        if i < length(signal2)
            shift_signal = real([signal2[end-i+1:end]; zeros(n-i)])
        else
            shift_signal = real([zeros(i-length(signal2)); signal2; zeros(n-i)])
        end
    end

    conv_sum = real(signal1).*shift_signal

#=
    n1 = length(signal1)
    n2 = length(signal2)
    if n2 < n1
        signal2 = [reverse(signal2);zeros(n1-n2)]
    elseif n1 < n2
        signal1 = [signal1;zeros(n2-n1)]
        signal2 = reverse(signal2)
    else
        signal2 = reverse(signal2)
    end
    if bounds_type == :periodic
        shift_signal = real([signal2[n1-i:n1-1];signal2[1:n1-i]])
    elseif bounds_type == :simple
        shift_signal = vcat(zeros(i), signal2[1:end-i])
    else
        println("bound type not found")
    end
=#

    conv_sum = real(signal1).*shift_signal

    plt = plot(real(signal1); label = "signal (f(x))",
               ylabel = "Signal and Filter values", dpi=200,
               ylims=(0,1.35), size=(600,400), right_margin=12mm,
               xlabel = "x", legend=:topright)
    plt = plot!(plt, shift_signal; label = "filter (g(x-t))")
    plt = plot!(plt, zeros(n); label = "area under f(t)g(x-t)",
                ribbon=(zeros(n),conv_sum))
    plt = plot!(plt, temp; label = out_name, linecolor=:purple,
                linestyle=:dash)

    plt = twinx()
    plt = plot!(plt, real(out[1:i]); label = "",
                linewidth=2, linestyle=:dash, ylabel = "Convolution output",
                grid=:off, ylims = (0,scale*1.35), linecolor=:purple)
    plt = plot!(plt, ([i,i],[0,real(out[i])/scale]); label = "",
                linewidth=2, linecolor=:purple)
#=
    plt = plot!(plt, vcat(zeros(i-1), real(out[i])/scale); label = "",
                linewidth=2, linecolor=:purple)
=#
    savefig(filebase*lpad(string(i),5,string(0))*".png")
end

function conv_fft(signal1::Array{Float64,1},
                  signal2::Array{Float64,1})
    return ifft(fft(signal1).*fft(signal2))
end

function conv_periodic(signal1::Array, signal2::Array; norm_factor = 1,
                  plot = false, frame_output = 1)

    n = length(signal1)

    out = Array{Float64,1}(undef,n)
    sum = 0

    for i = 1:n
        for j = 1:n
            if j < length(signal1) && mod1(i-j, n)+1 < length(signal2)
                sum += signal1[j] * signal2[mod1(i-j, n)+1]
            end
        end
        out[i] = sum

        if plot && i % frame_output == 0
            signal_plot(i, signal1, signal2, out; scale = norm_factor,
                        bounds_type = :periodic)
        end

        sum = 0
    end
    return out
end

function conv_lin(signal1::Array, signal2::Array;
                  norm_factor = 1, plot = false, frame_output = 1)

    n = length(signal1) + length(signal2)

    out = Array{Float64,1}(undef,n)
    sum = 0

    for i = 1:n
        for j = 1:i
            if j < length(signal1) && i-j+1 < length(signal2)
                sum += signal1[j] * signal2[i-j+1]
            end
        end
        out[i] = sum

        if plot && i % frame_output == 0
            signal_plot(i, signal1, signal2, out; scale = norm_factor,
                        bounds_type = :none)
        end

        sum = 0
    end

#=
    n = length(signal1)
    out = zeros(n)
    rsum = 0

    for i = 1:n-1
        if i+length(signal2) <= length(signal1)
            rsum = sum(signal1[i:i+length(signal2)-1].*reverse(signal2))
        elseif bounds_type == :periodic
            offset = i+length(signal2)-length(signal1)
            println("offset is ", offset)
            rsum = sum(vcat(signal1[i:end], signal1[1:offset-1])
                       .*reverse(signal2))
        elseif bounds_type == :simple
            diff = i+length(signal2)-length(signal1)
            rsum = sum(signal1[i:i+length(signal2)-1-diff].*
                       reverse(signal2)[1:end-diff])
        else
            println("bound type not found")
        end
        out[i] = rsum
        if plot && i % frame_output == 0
            signal_plot(i, signal1, signal2, out; scale = norm_factor,
                        bounds_type = bounds_type)
        end
        rsum = 0
    end
=#

    return out
end

function conv_check(signal1::Array, signal2::Array)
    n = length(signal1) + length(signal2)
    out = Array{Complex{Float64},1}(undef,n)
    sum = 0

    for i = 1:n-1
        for j = 1:i-1
            if i-j < length(signal2) && j < length(signal1)
                sum += signal1[j] * signal2[i-j]
            end
        end
        out[i] = sum
        sum = 0
    end

    return out
end

function main()
    # Gaussian signals
    #x = [exp(-((i-500)/1000)^2/.01) for i = 1:1000]
    y = [exp(-((i-50)/100)^2/.01) for i = 1:100]
    #y = vcat(y[div(length(y),2)+1:end], y[1:div(length(y),2)])

    # Square wave
    x = zeros(1000)
    x[400:600] .= 1.0 + 0im

    # Random
    x = rand(100)

    # Saw tooth
    #y = [(float(i))/200 for i = 200:-1:0]

    # Normalize
    #normalize!(x)
    #normalize!(y)
    out_conv = conv_periodic(x,y, norm_factor=sum(x.*x); plot=true,
                             frame_output = 1)
    #out_conv = conv_lin(x,y, norm_factor=sum(x.*x); plot=true,
    #                   frame_output = 1)

    return out_conv
end

max_stars_repo_name: leios/simuleios

max_issues_repo_name: leios/simuleios

max_forks_repo_name: leios/simuleios