text:
#This module works as a server for the calls from Java

module Ju4ja

using Sockets
using JSON

export startServer, closeServer, JavaCallInst, JavaCallResult

_server = Nothing
include("Ju4jaParser.jl")

#
function startServer(;port::Int=6996,
    funcProcessor::Function = executeFunction,
    scriptProcessor::Function = executeScript)
    _server = Sockets.listen(port)
    while true
        conn=accept(_server)
        @async begin
            try #执行一次客户端的请求过程，并返回结果
                line = readline(conn)
                println(line)
                exeResult=parseAndExecute(line,funcProcessor,scriptProcessor)
                write(conn,JSON.json(exeResult))
            catch ex
                println("Connection ended with error $ex")
            finally
                close(conn)
            end
        end
    end
end

#解析调用中的参数并执行相应的操作，返回结果
function parseAndExecute(line::String, funcProcessor::Function = executeFunction, scriptProcessor::Function = executeScript)
    result=JavaCallResult("","","","")
    jcall=parseLine(line)

    if jcall.operation == "FUNCTION"
        func=jcall.func
        if(func != "")
            try
                args=jcall.args #将参数试图解析为支持的类型
                res=funcProcessor(func,args)
                result.resultStr=res
                result.resultType = typeof(res)
                result.status="success"
            catch ex
                println(ex)
                result.resultStr=""
                result.status="error"
                #throw(ex)
            end
        end
    elseif jcall.operation == "SCRIPT"
        script=jcall.script
        #执行脚本
        try
            if (script != null)
                scriptProcessor(script)
                result.status= "success"
            else
                result.status= "null"
            end
        catch
            #println(ex)
            stacktrace()
            result.status="error"
        end
        result.resultStr=""
    end
    return result
end

function executeFunction(func::String,args::Array{Any})
    nargs = parseAllParams(args)
    ass = Tuple(nargs)
    println("Prepare to invoke function $func  with args: $ass")
    result=getfield(Main, Symbol(func))(ass...)
    return result
end

function executeScript(script::String)
    include(script)
end

function parseLine(ln::String)
    dict=JSON.parse(ln)
    jcall = JavaCallInst("","","","","","","")
    jcall.operation = get(dict,"operation","")
    jcall.script = get(dict,"script","")
    jcall.modn = get(dict,"modn","")
    jcall.func = get(dict,"func","")
    jcall.args = get(dict,"args",Any[])
    jcall.resultType = get(dict,"resultType","")
    if jcall.modn == ""
        jcall.modn = "Main"
    end
    println("JSON parsing finished successfuly.")
    return jcall
end

function closeServer()
    try
        if _server != Nothing
            close(_server)
        end
    catch ex
        print(ex)
    end
end

#对Java调用的参数解析
mutable struct JavaCallInst
    id#id
    operation#操作类型
    script#脚本名
    modn#模块名
    func#调用的函数名
    args#提供的参数
    resultType#返回值类型

end

#对Java调用的结果返回
mutable struct JavaCallResult
    id#
    resultStr#结果（返回值）
    resultType
    status#函数执行结果状态
end

end


text:
using Setfield

@testset "code2lens & lens2code" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2,"c":"oh"}]}""")
	j2 = JSON.parse("""{"a": [{"a":1,"b":3,"c":"hi"},{"b":2,"a":1,"c":"Mark"}]}""")
	j3 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")
	j4 = JSON.parse("""{"a": []}""")
	j5 = JSON.parse("""{}""")

	sch = JsonGrinder.schema([j1,j2,j3,j4, j5])
	ext = suggestextractor(sch)

    for t in list_traversal(sch)
        @test all(l -> all(t .== lens2code(sch, l)), code2lens(sch, t))
    end
    for t in list_traversal(ext)
        @test all(l -> all(t .== lens2code(ext, l)), code2lens(ext, t))
    end

	l_sch = code2lens(sch, "w") |> only
	l_ext = code2lens(ext, "w") |> only
	@test l_sch == (@lens _.childs[:a].items.childs[:c])
	@test l_ext == (@lens _.dict[:a].item.dict[:c])
	@test get(sch, l_sch) == sch["w"]
	@test get(ext, l_ext) == ext["w"]
	@test lens2code(sch, l_sch) |> only == "w"
	@test lens2code(ext, l_ext) |> only == "w"
end

ext: jl


text:
# TODO add support for reverse indexing with -1
#       e.g. a[0:-1:5]
# MOON add support for Julia matrix notation
#       e.g. a[5:0, 3:1, 2:0]

# Prerequitisite: Must be passed a preprocessed Juliog function block
function JULIOGtoJulia(ex::Expr, ad::Dict{Symbol, Arrow})
    global arrowDict = ad
    return JJfunction(ex)
end

# TODO add symbolic link for nameDict for direct assignments of one wire to another
function JJfunction(ex::Expr)
    println(ex)
    global nameDict = Dict{Symbol, Array{Symbol}}()
    jexpr = Expr(:Function)
    jargs = Array{Any}(2)
    jargs[1] = ex.args[1]
    jargs[2] = JJblock(ex.args[2])
    jexpr.args = jargs
    return jexpr
end

# Prerequisite: Expr head MUST be a :block
function JJblock(ex::Expr)
    l = length(ex.args)
    jexpr = Expr(:block)
    jargs = Array{Any}(l)

    for i = 1:l
        println(i)
        arg = ex.args[i]
        println(arg)
        if     isa(arg, Expr)
            jargs[i] = JJblockhelper(arg)
            println("jargs at $(i) is $(jargs[i])")
        else
            error("Inside JULIOGtoJulia hit a non-expression block argument")
        end
    end

    # remove :block exprs
    i = 1
    while i <= length(jargs)
        arg = jargs[i]
        if     isa(arg, Expr) && arg.head == :block
            splice!(jargs, i, arg.args)
        elseif isa(arg, Expr) && arg.head == :call
            i = i + 1
        else
            error("A non Expr, or an Expr not of type call is found in the Julia:\n$(jargs)")
        end
    end

    jexpr.args = jargs
    return jexpr
end

function JJblockhelper(ex::Expr)
    h = ex.head
    elseif h == :(=)
        return JJequals(ex)
    elseif h == :(:=)
        return JJassignment(ex)
    elseif h == :macrocall
        return JJmacro(ex)
    else
        error("Hit an unexplored expression head $(h) in expr:\n$(ex)")
    end
end

function JJequals(ex::Expr)
    global arrowDict

    # Get the name of the new wire or reg
    if isa(ex.args[1], Expr)
        # Is a ref
        name = ex.args[1].args[1]
        # MOON
    else
        name = ex.args[1]
    end

    # Get the arrow info
    a = arrowDict[name]

    # Get the wire assignment information
    asgn = a.Assigned
    # Get the wire type
    wt = a.WireType
    # Get the bit indices and string
    bc = a.BitCount

    if     asgn == 0
        jargs = Array{Any}(2)
        jargs[1] = Expr(:call, :Wire, name, bc)
        jargs[2] = Expr(:call, :CONST, name, fill!(Array{UIn8}(bc), UInt8('z')))
        jexpr = Expr(:block, jargs[1], jargs[2])
        return jexpr
    elseif asgn == 1
        if     wt == 0
            return Expr(:call, :Wire, name, bc)
        elseif wt == 1
            return Expr(:call, :Input, name, bc)
        elseif wt == 2
            return Expr(:call, :Output, name, bc)
        end
    elseif asgn == 2
        global nameDict
        jargs = Array{Any}(0)
        nameArray = Array{Symbol}(0)

        tmpasgn = a.AsgnArray[1]
        tmpi = 1
        bmin = min(a.BitLeft, a.BitRight)

        l = length(a.AsgnArray)
        for i = 2:l
            asgn = a.AsgnArray[i]
            if asgn != tmpasgn
                newname = Symbol(
                    String(String(name) * "!!b" *
                    dec(tmpi + bmin - 1) * "_" * dec(i + bmin - 2)))
                bc = i - tmpi
                append!(nameArray, fill(newname , bc))

                if   asgn == 0
                    append!(jargs, [Expr(:call, :Wire, newname, bc) ;
                        Expr(:call, :CONST, newname, fill(UInt8('z'), bc))])
                else
                    push!(jargs, Expr(:call, :Wire, newname, bc))
                end

                tmpasgn = asgn
                tmpi = i
            end
        end
        # Solve for last wire bits
        newname = Symbol(
            String(String(name) * "!!b" *
            dec(tmpi + bmin - 1) * "_" * dec(l + bmin - 1)))
        bc = l - tmpi + 1
        append!(nameArray, fill(newname , bc))

        if   asgn == 0
            append!(jargs, [Expr(:call, :Wire, newname, bc) ;
                Expr(:call, :CONST, newname, fill(UInt8('z'), bc))])
        else
            push!(jargs, Expr(:call, :Wire, newname, bc))
        end

        nameDict[name] = nameArray

        jexpr = Expr(:block)
        jexpr.args = jargs
        return jexpr
    end
end

function JJassignment(ex::Expr)

    # Deal with lhs
    if isa(ex.args[1], Expr)
        # Is a ref
        name = ex.args[1].args[1]
        # TODO get name for this assignment from new data structure
        # TODO account for reverse indexing

        # MOON
    else
        # Must just be a symbol and therefore the wire is not partitioned
        name = ex.args[1]
    end

    # Deal with rhs
    rhs = ex.args[2]
    # If declaring two wires as equal, make a symbolic link and save bit offset
    if     isa(rhs, Symbol)

    # Also a symbolic link
    elseif isa(rhs, Expr) && rhs.head == :ref && isa(rhs.args[1], Symbol)

    # Is a simple constant declaration
    elseif isa(rhs, Int)

    # Is a simple constant declaration
    elseif isa(rhs, Char)

    # Actually some logic involved
    else
        rhsArgs, rhsJulia = JJrhs(rhs)
    end

    # Assign to lhs
    if isa(ex.args[1], Expr)
        # We are indexing bit(s)
        sym = ex.args[1].args[1]

        # Get bits of interest
        if isa(ex.args[1].args[2], Expr)
            # Multiple Bits refererence
            bl = ex.args[1].args[2].args[1]
            br = ex.args[1].args[2].args[2]
            endi = bl > br
        else
            # Bit reference
            bl = ex.args[1].args[2]
            br = bl
            endi = true
        end

        # throw error if writing bits of undefined wire
        if !haskey(arrowdict, sym)
            error("Attempting to write to the bits of undefined wire: $(sym)")
        else
            # throw error if we are assigning to these wires an additional time
            arrow = arrowdict[sym]
            l = length(arrow.SubArrows)
            if l != 0
                if endi
                    # Big endian
                    for i = 1:l
                        sa = arrow.SubArrows[i]
                        if ((bl >= sa.bitleft) && (sa.bitleft >= br))
                            || ((bl >= sa.bitright) && (sa.bitright >= br))
                            error("Attempting to assign bits of wire $(sym) an additional time")
                        end
                    end
                else
                    # little endian
                    for i = 1:l
                        sa = arrow.SubArrows[i]
                        if ((br >= sa.bitleft) && (sa.bitleft >= bl))
                            || ((br >= sa.bitright) && (sa.bitright >= bl))
                            error("Attempting to assign bits of wire $(sym) an additional time")
                        end
                    end
                end
            end
        end

        # Done error checking
        # TODO create wire(s) as needed for correct bifurcation
        # TODO evaluation
        # something like
        # JJevaluation(sym , ex, arrowdict)

    else
        # Are not bit indexing
        sym = ex.args[1]

        if haskey(arrowdict, sym)
            # throw error if writing to a defined wire we have already written over
            # throw error if we are assigning to these wires an additional time
            arrow = arrowdict[sym]
            l = length(arrow.SubArrows)
            if l != 0
                error("Attempting to assign wire $(sym) an additional time")
            end

            # TODO evaluation
            # something like
            # JJevaluation(sym, ex, arrowdict)

        else
            # determine bit count of rhs
            # TODO

        end
    end

end

function JJrhs(rhs)
    if isa(rhs, Expr)
        bcrhs, endirhs = JSrhsexpr(rhs)
    else
        # is a Symbol or 0,1,'x','z'
        if isa(rhs, Symbol)
            if !haskey(arrowDict, rhs)
                error("Discovered uninstantiated symbol: $(rhs) in expr:\n$(ex)")
            end
            bcrhs =  arrowDict[rhs].BitCount
            endirhs = arrowDict[rhs].Endian
        else
            if     rhs == 0
                bcrhs = 1
                endirhs = 0
            elseif rhs == 1
                bcrhs = 1
                endirhs = 0
            elseif rhs == 'x'
                bcrhs = 1
                endirhs = 0
            elseif rhs == 'z'
                bcrhs = 1
                endirhs = 0
            else
                error("Attempting to assign wire to incorrectly syntaxed Right Hand Side:\n$(ex)")
            end
        end
    end
end

# The righthand side of an assignment
function JJrhsexpr(ex::Expr)
    h = ex.head
    if     h == :ref
        return JJrhsref(ex)
    elseif h == :vcat
        return JJrhsvcat(ex)
    elseif h == :call
        return JJrhscall(ex)
    elseif h == :if
        return JJrhsif(ex)
    else
        error("Hit unexpected expression symbol $(h) in JSreadrhs on expression:\n$(ex)")
    end
end

function JJevaluation(sym::Symbol, bl::Int, br::Int, endi::Bool, ex::Expr)
    jexpr = Expr(:block)
    newsym = Symbol(String(sym)  * "!" * dec(bl) * "_" * dec(br))
    b = abs(bl - br) + 1
    wirearg = [Symbol("@wire") ; newsym ; b]

    if !isa(ex.args[2], Expr)
        # is either a constant wire, or a complete assignment of one wire to another
        if isa(ex.args[2], Symbol)
            # Assignment of one wire to another
            jargs = Array{Any}(2)
            jargs[1] = wirearg
            jargs[2] = [Symbol("EQUALS") ; newsym ; ex.args[2]]
            jexpr.args = jargs
            return jexpr
        else
            # constant (unsigned) value
            val = ex.args[2]
            vt = typeof(val)
            if vt != Int && vt != UInt && vt != UInt8 && vt != Array{UInt8}
                error("Attempting to write a constant value which is of wrong type")
            end
            if vt == Int
                val = UInt8(val)
            end
            jargs = Array{Any}(2)
            newsym = Symbol(String(sym)  * "!" * dec(bl) * "_" * dec(br))
            jargs[1] = wirearg
            jargs[2] = [Symbol("CONST") ; newsym ; val]
            jexpr.args = jargs
            return jexpr
        end
    end

    # Must be an expr on rhs
    #  descend into expressions
    finalsuffix, bitcount, evalargs = JJevaluationhelper(0, newsym, ex.args[2])

    # raise block that was created
    # Make certain that lhs and rhs have same bitcount

    # Return expr and wire sym
end

# Always passed an expr, and always the rhs of an assignment
function JJevaluationhelper(suffix::Int, newsym::Symbol, endi::Bool, ex::Expr)
    h = ex.head
    if     h == :ref
        # Is a bitwise wire index
        return JJref(suffix, newsym, ex, arrowdict)
    elseif h == :vcat
        # Is a wire concatenation
        return JJvcat(suffix, newsym, endi, ex, arrowdict)
    elseif h == :call
        l = length(ex.args)
        s = Array{Int}(l+1)
        s[1] = suffix
        bt = Array{Int}(l)
        jargs = Array{Any}(0)
        for i = 2:l
            if     isa(ex.args[i], Expr)
                s[i+1], bt[i], evalargs = JJevaluationhelper(s[i], newsym, endi, ex)
                append!(jargs, evalargs)
            end
        end

        op = ex.args[1]
        if isa(op, Symbol)
            if     op == :(&)
                # Must have 2 operands
                if bt[1] != bt[2]
                    error("Attempting to AND wires of differing bit sizes: $(bt[1]) and $(bt[2])")
                end
                andargs = Array{Any}(2)

            elseif op == :(|)

            elseif op == :(~)

            elseif op == :(^)

            elseif op == :(+)

            elseif op == :(-)

            elseif op == :(*)

            elseif op == :(/)

            elseif op == :(%)

            else
                # must be the bit repeat of a wire

            end

        else
            # Must be a bit repeat of a bit reference

        end

    else
        error("In JJevaluationhelper, hit an unexpected subexpression head: $(h)")
    end
end

# TODO add support for matrix notation
function JJref(suffix::Int, newsym::Symbol, ex::Expr)
    # Since this is an expr, it must be a bit(s) reference
    refsym = ex.args[1]
    if !haskey(arrowdict, refsym)
        error("Attempting to reference an undefined wire: $(refsym)")
    end
    refarrow = arrowdict[refsym]

    # Create a symbol for this intermediate wire
    interwiresym = Symbol(String(newsym) * "!!" * dec(suffix))

    if isa(ex.args[2], Expr)
        # bits reference
        bl = ex.args[2].args[1]
        br = ex.args[2].args[2]

        if refarrow.Endian
            # Big-Endian

            # Check endian of reference
            if br > bl
                error("Attempting to reference $(refsym) in incorrect endian, reverse endian not implemented yet")
            end

            # Check boundary conditions
            if     bl > refarrow.BitLeft
                error("Attempting to reference bit $(bl) which is outside $(refsym)")
            elseif br < refarrow.BitRight
                error("Attempting to reference bit $(br) which is outside $(refsym)")
            end

            # Create new wire
            bt = (bl-br)+1
            wireexpr = Expr(:macrocall)
            wireexpr.args = [Symbol("@wire") ; interwiresym ; bt]

            # BITS new wire
            refexpr = Expr(:call)
            refexpr.args = [Symbol("BITS") ; interwiresym ; refsym ; bl-refarrow.BitRight ; br-refarrow.BitRight ]

            jargs = Array{Any}(2)
            jargs[1] = wireexpr
            jargs[2] = refexpr
            return (suffix+1, bt, jargs)

        else
            # Little-Endian

            # Check endian of reference
            if bl > br
                error("Attempting to reference $(refsym) in incorrect endian")
            end

            # Check boundary conditions
            if     bl < refarrow.BitLeft
                error("Attempting to reference bit $(bl) which is outside $(refsym)")
            elseif br > refarrow.BitRight
                error("Attempting to reference bit $(br) which is outside $(refsym)")
            end

            # Create new wire
            # Create expr
            bt = (br-b1)+1
            wireexpr = Expr(:macrocall)
            wireexpr.args = [Symbol("@wire") ; interwiresym ; bt]

            # BITS new wire
            refexpr = Expr(:call)
            refexpr.args = [Symbol("BITS") ; interwiresym ; refsym ; br-refarrow.BitLeft ; bl-refarrow.BitLeft ]

            jargs = Array{Any}(2)
            jargs[1] = wireexpr
            jargs[2] = refexpr
            return (suffix+1, bt, jargs)
        end

    else
        # single bit reference
        b = ex.args[2]

        # Check boundary conditions
        if refarrow.Endian
            if b > refarrow.BitLeft
                error("Attempting to reference bit $(b) which is outside $(refsym)")
            end

            # BIT new wire
            refexpr = Expr(:call)
            refexpr.args = [Symbol("BIT") ; interwiresym ; refsym ; b-refarrow.BitRight]

        else
            if b > refarrow.BitRight
                error("Attempting to reference bit $(b) which is outside $(refsym)")
            end
            # BIT new wire
            refexpr = Expr(:call)
            refexpr.args = [Symbol("BIT") ; interwiresym ; refsym ; b-refarrow.BitLeft]
        end

        # Create new wire
        wireexpr = Expr(:macrocall)
        wireexpr.args = [Symbol("@wire") ; interwiresym ; 1]

        jargs = Array{Any}(2)
        jargs[1] = wireexpr
        jargs[2] = refexpr
        return (suffix+1, 1, jargs)
    end
end

# TODO add support for bit repeat
# only accepts references, symbols, chars, ints
# need to know endian of wire being written
function JJvcat(suffix::Int, newsym::Symbol, endi::Bool, ex::Expr)
    l = length(ex.args)

    if l == 1
        error("Attempting to concatenate a single wire: $(ex.args[1])\n\tDon't do that.")
    end

    s = suffix

    totalbitcount = 0
    bitcount
    bitval = 0
    i = endi ? 1 : l
    jargs = Array{Any}(0)
    while (endi && (i <= l)) || (!endi && (i > 0))
        arg = ex.args[i]

        if     isa(arg, Int) && (arg == 0 || arg == 1)
            # accumulate bit string until something other than int
            bitcount = bitcount + 1
            bitval = (bitval << 1) | arg
            totalbitcount = totalbitcount + 1
        else
            if s == suffix && bitcount == 0
                catargs = Array{Any}(2)
                vcatsym = Symbol(String(newsym)  * "!!" * dec(s))
                if     isa(arg, Expr)
                    # Must be a bit reference
                    tmp, bt, refargs = JJref(s, newsym, ex)
                    catargs[1:2] = refargs
                    totalbitcount = totalbitcount + bt
                elseif isa(arg, Symbol)
                    # Must be a wire symbol
                    # TODO flesh this out
                    error("Until Reverse endian indexing is implemented, You must bit reference any wire in a vcat")
                else
                    # Must be an int not 0 or 1
                    argtype = typeof(arg)
                    if argtype != Int && argtype != Char && argtype != UInt && argtype != UInt8 && argtype != Array{UInt8}
                        error("Attempting to write a constant value which is of wrong type")
                    end
                    if argtype == Int || argtype == Char
                        val = UInt8(arg)
                    end

                    wireexpr = Expr(:macrocall)
                    wireexpr.args = [Symbol("@wire") ; vcatsym ; 1]
                    constexpr = Expr(:call)
                    constexpr.arg = [Symbol("CONST") ; vcatsym ; val]
                    catargs[1] = wireexpr
                    catargs[2] = constexpr
                    bitcount = 0
                    bitval = 0
                    totalbitcount = totalbitcount + 1
                end
                append!(jargs, catargs)
                s = s + 1
            else
                if bitcount != 0
                    # Means we have accumulated integers
                    vcatsym = Symbol(String(newsym)  * "!!" * dec(s))
                    wireexpr = Expr(:macrocall)
                    wireexpr.args = [Symbol("@wire") ; vcatsym ; bitcount]
                    constexpr = Expr(:call)
                    constexpr.arg = [Symbol("CONST") ; vcatsym ; bitval]

                    constargs = Array{Any}(2)
                    constargs[1] = wireexpr
                    constargs[2] = constexpr
                    append!(jargs, constargs)
                    s = s + 1

                    bitcount = 0
                    bitval = 0
                    totalbitcount = totalbitcount + 1
                end

                catargs = Array{Any}(4)
                vcatsym = Symbol(String(newsym)  * "!!" * dec(s))
                if     isa(arg, Expr)
                    # Must be a bit reference
                    tmp, bt, refargs = JJref(s, newsym, ex)
                    constargs[1:2] = refargs
                    totalbitcount = totalbitcount + bt
                elseif isa(arg, Symbol)
                    # Must be a wire symbol
                    error("Until Reverse endian indexing is implemented, You must bit reference any wire in a vcat")
                else
                    # Must be an int not 0 or 1
                    argtype = typeof(arg)
                    if argtype != Int && argtype != Char && argtype != UInt && argtype != UInt8 && argtype != Array{UInt8}
                        error("Attempting to write a constant value which is of wrong type")
                    end
                    if argtype == Int || argtype == Char
                        val = UInt8(arg)
                    end

                    wireexpr = Expr(:macrocall)
                    wireexpr.args = [Symbol("@wire") ; vcatsym ; 1]
                    constexpr = Expr(:call)
                    constexpr.arg = [Symbol("CONST") ; vcatsym ; val]
                    catargs[1] = wireexpr
                    catargs[2] = constexpr
                    totalbitcount = totalbitcount + 1
                end

                vcatsymtwo = Symbol(String(newsym)  * "!!" * dec(s+1))
                catwireexpr = Expr(:macrocall)
                catwireexpr.args = [Symbol("@wire") ; vcatsymtwo ; totalbitcount]
                catcallexpr = Expr(:call)
                catcallexpr.args = [Symbol("CAT") ; vcatsymtwo; vcatsym ; Symbol(String(newsym)  * "!!" * dec(s-1))]

                catargs[3] = catwireexpr
                catargs[4] = catcallexpr
                append!(jargs, catargs)
                s = s + 2
            end
        end

        if endi
            i = i + 1
        else
            i = i - 1
        end
    end

    # In case the concatenation ends with bit padding
    if bitcount != 0
        catargs = Array{Any}(4)
        vcatsym = Symbol(String(newsym)  * "!!" * dec(s))
        constwireexpr = Expr(:macrocall)
        constwireexpr.args = [Symbol("@wire") ; vcatsym ; bitcount]
        constcallexpr = Expr(:call)
        constcallexpr.arg = [Symbol("CONST") ; vcatsym ; bitval]

        vcatsymtwo = Symbol(String(newsym)  * "!!" * dec(s+1))
        catwireexpr = Expr(:macrocall)
        catwireexpr.args = [Symbol("@wire") ; vcatsymtwo ; totalbitcount]
        catcallexpr = Expr(:call)
        catcallexpr.args = [Symbol("CAT") ; vcatsymtwo; vcatsym ; Symbol(String(newsym)  * "!!" * dec(s-1))]

        catargs[1] = constwireexpr
        catargs[2] = constcallexpr
        catargs[3] = catwireexpr
        catargs[4] = catcallexpr
        append!(jargs, catargs)
    end

    # Return the expressions
    return (s, totalbitcount, jargs)
end

function JJmacro(ex::Expr)
    mc = ex.args[1]
    if     mc == Symbol("@async")
    error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@reg")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@posedge")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@negedge")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@delay")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@block")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@verilog")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@Julia")
        error("Hit an unfinished macrocall: $(mc)")
    else
        error("Hit an unexplored macrocall $(mc) in expr:\n$(ex)")
    end
end

#=

    # Can't be reduced to a single wire, so create a new wire
    wireName = "if!" * toString(ex.args[2]) * "_" * toString(ex.args[3]) * "!!"
    i = 0
    wireNamei = Symbol(wireName * dec(i))
    while haskey(arrowDict, wireNamei)
        i = i + 1
        wireNamei = Symbol(wireName * dec(i))
    end

    jargs = Array{Any}(2)
    jargs[1] = Expr(:(=), wireNamei, Expr(:(ref), Expr(:call, :Wire), 0))
    jargs[2] = Expr(:(:=), wireNamei, ex)

function toString(name)
    if isa(name ,Symbol)
        return String(name)
    elseif isa(name, Int)
        return dec(name)
    else
        return string(name)
    end
end

=#

max_stars_repo_name: trevortheblack/Juliog

max_issues_repo_name: trevortheblack/Juliog

max_forks_repo_name: trevortheblack/Juliog


text:
module QFT

using QXZoo.GateOps
using QXZoo.Circuit
using QXZoo.DefaultGates
using QXZoo.CompositeGates

export apply_qft!, apply_iqft!

"""
    apply_qft!(cct::Circuit.Circ, qubit_indices::Vector)

Apply the quantum Fourier transform over the given qubit index range.
"""
function apply_qft!(cct::Circuit.Circ, qubit_indices::Vector)
    for idx in 1:length(qubit_indices)
        cct << DefaultGates.h(qubit_indices[idx])
        for jdx in idx+1:length(qubit_indices)
            cct << DefaultGates.c_r_phase(qubit_indices[idx], qubit_indices[jdx], 2π/(2^(qubit_indices[jdx]-qubit_indices[idx]+1)))
        end
    end
    lhs, rhs = swap_idx(qubit_indices)
    for i in 1:length(lhs)
        cct << CompositeGates.swap(lhs[i], rhs[i])
    end
    ;
end
"""
    apply_qft!(cct::Circuit.Circ)

Apply the quantum Fourier transform over the entire qubit register
"""
apply_qft!(cct::Circuit.Circ) = apply_qft!(cct, collect(1:cct.num_qubits))

"""
    apply_qft!(cct::Circuit.Circ, qubit_indices::Vector)

Apply the inverse quantum Fourier transform over the given qubit index range.
"""
function apply_iqft!(cct::Circuit.Circ, qubit_indices::Vector)
    lhs, rhs = swap_idx(qubit_indices)

    for i in 1:length(lhs)
        cct << CompositeGates.swap(lhs[i], rhs[i])
    end

    for idx in length(qubit_indices):-1:1
        for jdx in length(qubit_indices):-1:idx+1
            cct << DefaultGates.c_r_phase(qubit_indices[idx], qubit_indices[jdx], -2π/(2^(qubit_indices[jdx]-qubit_indices[idx]+1)))
        end
        cct << DefaultGates.h(qubit_indices[idx])
    end

    ;
end
"""
    apply_iqft!(cct::Circuit.Circ)

Apply the inverse quantum Fourier transform over the entire qubit register
"""
apply_iqft!(cct::Circuit.Circ) = apply_iqft!(cct, collect(1:cct.num_qubits))

"""
    swap_idx(qubit_indices::Vector)

Return the indices of the qubits required to apply the swap operations as part of the transform.
"""
function swap_idx(qubit_indices::Vector)
    is_even = (length(qubit_indices) & 1 === 0)
    mp = convert(Int, floor(length(qubit_indices)//2))
    lhs = @view qubit_indices[1:mp]
    if is_even
        rhs = @view qubit_indices[end:-1:mp+1]
    else
        rhs = @view qubit_indices[end:-1:mp+2]
    end
    return lhs, rhs
end

end


text:
using Test
using QuantumOptics

# Implement custom type with AbstractArray interface
mutable struct TestData{T,N,X} <: AbstractArray{T,N}
    x::X
    function TestData(x::X) where X
        x_ = complex.(x)
        new{ComplexF64,length(axes(x)),typeof(x_)}(x_)
    end
end
Base.size(A::TestData) = size(A.x)
Base.getindex(A::TestData, inds...) = getindex(A.x, inds...)
Base.setindex!(A::TestData, val, inds...) = setindex!(A.x, val, inds...)
Base.similar(A::TestData, dims::Int...) = TestData(similar(A.x,dims...))

@testset "abstract-data" begin

###############
# Schrödinger #
N = 3
Ncutoff = 2
T = [0.:0.1:1.;]

Ω = [0 2 3;
     2 0 1;
     3 1 0]

ω = [1., 1.2, 1.5]

basis_fock = FockBasis(Ncutoff)
basis = tensor([basis_fock for i=1:N]...)

a = destroy(basis_fock)
at = create(basis_fock)
I = identityoperator(basis_fock)

psi0 = tensor([coherentstate(basis_fock, i%Ncutoff) for i=1:N]...)

# Interaction picture
Hrot = SparseOpType[]
for i=1:N, j=1:N
    if i==j
        continue
    end
    h = embed(basis, [i,j], [a, Ω[i,j]*at])
    push!(Hrot, h)
end
Hrot = sum(Hrot)

# Schroedinger picture
function f(t, psi)
    H = SparseOpType[embed(basis, i, ω[i]*at*a) for i=1:N]
    for i=1:N, j=1:N
        if i==j
            continue
        end
        h = embed(basis, [i,j], [a, exp(1im*(ω[i]-ω[j])*t)*Ω[i,j]*at])
        push!(H, h)
    end
    sum(H)
end

Hrot_ = Operator(Hrot.basis_l,Hrot.basis_r,TestData(Hrot.data))
tout, psi_rot_t = timeevolution.schroedinger(T, psi0, Hrot)
tout, psi_rot_t_ = timeevolution.schroedinger(T, psi0, Hrot_)
for (psi1,psi2)=zip(psi_rot_t,psi_rot_t_)
    @test psi1==psi2
end

psi0_ = Ket(basis, TestData(psi0.data))
tout, psi_rot_t_ = timeevolution.schroedinger(T, psi0_, Hrot_)
for (psi1,psi2)=zip(psi_rot_t,psi_rot_t_)
    @test psi1.data≈psi2.data
end

##########
# Master #

ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1

T = Float64[0.,1.]

fockbasis = FockBasis(10)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H_ = Ha + Hc + Hint
H = Operator(basis, TestData(H_.data))

Ja_unscaled = embed(basis, 1, sm)
Jc_unscaled = embed(basis, 2, destroy(fockbasis))
Junscaled = [Ja_unscaled, Jc_unscaled]

Ja_ = embed(basis, 1, sqrt(γ)*sm)
Jc_ = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
Ja = Operator(basis, TestData(Ja_.data))
Jc = Operator(basis, TestData(Jc_.data))
J = [Ja, Jc]
Jlazy = [LazyTensor(basis, 1, sqrt(γ)*Operator(spinbasis, TestData(sm.data))), Jc]

Hnh = H - 0.5im*sum([dagger(J[i])*J[i] for i=1:length(J)])

Hdense = dense(H)
Hlazy = LazySum(Ha, Hc, Hint)
Hnh_dense = dense(Hnh)
Junscaled_dense = map(dense, Junscaled)
Jdense = map(dense, J)

Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 5)
ρ₀ = Operator(basis, TestData(dm(Ψ₀).data))

# Test Liouvillian
L = liouvillian(H, J)
ρ = -1im*(H*ρ₀ - ρ₀*H)
for j=J
    ρ .+= j*ρ₀*dagger(j) - 0.5*dagger(j)*j*ρ₀ - 0.5*ρ₀*dagger(j)*j
end
@test tracedistance(L*ρ₀, ρ) < 1e-10

# Test master
tout, ρt = timeevolution.master(T, ρ₀, Hdense, Jdense; reltol=1e-7)
ρ = ρt[end]
@test isa(ρ.data, TestData)
@test tracedistance(dense(exp(dense(L)*T[end])*ρ₀), dense(ρ)) < 1e-6

@test isa(ρ₀.data, TestData) && isa(H.data,TestData) && all(isa(j.data,TestData) for j=J)
tout, ρt = timeevolution.master(T, ρ₀, H, J; reltol=1e-6)
@test isa(ρt[end].data, TestData)
@test tracedistance(dense(ρt[end]), dense(ρ)) < 1e-5

########
# MCWF #

# Define parameters for spin coupled to electric field mode.
ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1

Ntrajectories = 500
T = Float64[0.:0.1:10.;]

# Define operators
fockbasis = FockBasis(8)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

# Hamiltonian
Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H_ = Ha + Hc + Hint
H = Operator(basis, TestData(H_.data))
Hdense = dense(H)

# Jump operators
Ja = embed(basis, 1, sqrt(γ)*sm)
Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
J_ = [Ja, Jc]
J = [Operator(basis, TestData(j.data)) for j=J_]
Jdense = map(dense, J)

# Initial conditions
Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 5)

# Test mcwf
tout, Ψt = timeevolution.mcwf(T, Ψ₀, Hdense, Jdense; seed=UInt(1), reltol=1e-7)
tout2, Ψt2 = timeevolution.mcwf(T, Ψ₀, Hdense, Jdense; seed=UInt(1), reltol=1e-7)
@test Ψt == Ψt2
Ψ = Ψt[end]

tout, Ψt = timeevolution.mcwf(T, Ψ₀, H, J; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf(T, Ψ₀, H, J; seed=UInt(2), reltol=1e-6)
@test norm(Ψt[end]-Ψ) > 0.1

#####################
# Two-level example #

basis = SpinBasis(1//2)

# Random 2 level Hamiltonian
a1 = 0.5
a2 = 1.9
c = 1.3
d = -4.7

data = [a1 c-1im*d; c+1im*d a2]
H = DenseOperator(basis, data)

a = (a1 + a2)/2
b = (a1 - a2)/2
r = [c d b]

sigma_r = c*sigmax(basis) + d*sigmay(basis) + b*sigmaz(basis)

U(t) = exp(-1im*a*t)*(cos(norm(r)*t)*one(basis) - 1im*sin(norm(r)*t)*sigma_r/norm(r))

# Random initial state
psi0 = randstate(basis)
T = [0:0.5:1;]

f(t, psi::Ket) = @test 1e-5 > norm(psi - U(t)*psi0)
H_ = Operator(basis,TestData(H.data))
timeevolution.schroedinger(T, psi0, H_; fout=f)
timeevolution.mcwf(T, psi0, H_, []; fout=f)

f(t, rho) = @test 1e-5 > tracedistance(dense(rho), dm(U(t)*psi0))
rho0 = Operator(basis, TestData(dm(psi0).data))
timeevolution.master(T, rho0, H_, []; fout=f)

end # testset


text:
function mk_netsta(S::GphysData)
  netsta = Array{String,1}(undef, S.n)
  cha = Array{String,1}(undef, S.n)
  for i in 1:S.n
    id = split_id(S.id[i])
    netsta[i] = id[1]*"."*id[2]
    cha[i] = lowercase(id[4])
  end
  nsid = unique(netsta)
  return netsta, cha, nsid
end

function asdf_wsxml(xbuf::IOBuffer, S::GphysData, chans::Array{Int64,1}, sta::HDF5Group)
  seekstart(xbuf)
  mk_xml!(xbuf, S, chans)
  sta["StationXML"] = take!(xbuf)
  return nothing
end

function asdf_mktrace(S::GphysData, xml_buf::IO, chan_numbers::Array{Int64,1}, wav::HDF5Group, ts::Array{Int64,1}, te::Array{Int64,1}, len::Int64, v::Integer, tag::String)
  nc = length(chan_numbers)
  netsta, cha, nsid = mk_netsta(S)
  trace_names = Array{Array{String,1},1}(undef, nc)
  (v>2) && println("traces to create = channels ", chan_numbers)

  # Build trace_names
  d0 = zeros(Int64, nc)
  for (i,j) in enumerate(chan_numbers)
    trace_names[i] = String[]
    id = S.id[j]
    cc = isempty(tag) ? cha[j] : tag

    # divide t0 and t1 by the time length using div to get d0, d1
    d0[i] = len*(div(ts[i], len))
    d1 = len*(div(te[i], len))
    Δ = round(Int64, 1.0e9/S.fs[j])

    # range of times is then d0:d1
    for d in d0[i]:len:d1
      # ...but we must end Δ ns *before* d1 to prevent a one-sample overlap
      s0 = string(u2d(d/1000000000))
      s1 = string(u2d((d+len-Δ)/1000000000))

      # create string like CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_
      chan_str = join([id, s0, s1, cc], "__")
      push!(trace_names[i], chan_str)
    end
    (v > 2) && println("trace_names[", i, "] = ", trace_names[i])
  end

  # Check for trace names in each ns
  chans = falses(nc)
  for j in 1:length(nsid)
    id = nsid[j]
    fill!(chans, false)
    for (j,i) in enumerate(chan_numbers)
      if (netsta[i] == id) && (S.fs[i] != 0.0)
        chans[j] = true
      end
    end
    (maximum(chans) == false) && continue

    # does net.sta exist?
    if has(wav, id)
      sta = wav[id]

    else
      # Create Waveforms/net.sta
      (v > 0) && println("creating Waveforms/", id)
      sta = g_create(wav, id)

      # Write StationXML to sta
      asdf_wsxml(xml_buf, S, chan_numbers[chans], sta)
    end

    # Create Waveforms/net.sta/chan_str
    for i = 1:length(chans)
      if chans[i]
        j = chan_numbers[i]
        (v > 2) && println("S.id[", j, "] = ", S.id[j], ", trace_names[", i, "] = ", trace_names[i])
        T = eltype(S.x[j])
        fs = S.fs[j]
        Δ = round(Int64, 1.0e9/fs)
        nx = div(len, Δ)
        for k = 1:length(trace_names[i])
          chan_str = trace_names[i][k]

          # is there a trace with the corresponding string?
          if has(sta, chan_str) == false
            (v > 1) && println("creating trace:  Waveforms/", id, "/", chan_str)
            sta[chan_str] = ones(T, nx).*T(NaN)

            # create a trace of all NaNs
            attrs(sta[chan_str])["sampling_rate"] = fs
            attrs(sta[chan_str])["starttime"] = d0[i] + (k-1)*len
          end
        end
      end
    end
  end

  return nothing
end

function asdf_write_chan(S::GphysData, sta::HDF5Group, i::Int64, tag::String, eid::String, v::Integer)
  fs = S.fs[i]
  tx = S.t[i]
  t = t_win(tx, fs)
  xi = x_inds(tx)
  n_seg = size(t, 1)
  (n_seg == 0) && return

    for k = 1:n_seg
    t0 = t[k,1]
    t1 = t[k,2]
    s0 = string(u2d(div(t0, 1000000)))
    s1 = string(u2d(div(t1, 1000000)))

    # create string like CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_
    chan_str = join([S.id[i], s0, s1, tag], "__")
    if has(sta, chan_str)
      (v > 0) && println("incrementing tag of ", chan_str)
      j = 0x2f
      while has(sta, chan_str)
        j += 0x01
        tag1 = String(vcat(UInt8.(codeunits(tag)), j))
        chan_str = join([S.id[i], s0, s1, tag1], "__")
      end
    end
    sta[chan_str] = S.x[i][xi[k,1]:xi[k,2]]

    # set dictionary attributes
    D = attrs(sta[chan_str])
    D["sampling_rate"] = fs
    D["starttime"] = t0*1000

    if !isempty(eid)
      D["event_id"] = eid
    end
  end

  return nothing
end

function write_asdf( hdf_out::String,
                     S::GphysData,
                     chan_numbers::Array{Int64,1},
                     add::Bool,
                     evid::String,
                     ovr::Bool,
                     len::Period,
                     tag::String,
                     v::Integer)

  # "add" implies "ovr"
  if add == true
    ovr = true
  end

  # Precheck for degenerate time structs, time shift
  if ovr
    for i in chan_numbers
      if (size(S.t[i],1) < 2) && (S.fs[i] != 0.0)
        error(string(S.id[i] *
          ": malformed :t; can't write with ovr=true."))
      end

      # ensure each channel starts on an exact sample
      t0 = S.t[i][1,2]
      Δ = round(Int64, sμ/S.fs[i])
      t1 = div(t0,Δ)*Δ
      δ = t0-t1
      S.t[i][1,2] = t1
      S.misc[i]["tc"] = δ
    end
  end

  xml_buf = IOBuffer()
  netsta, cha, nsid = mk_netsta(S)

  if isfile(hdf_out)
    io = h5open(hdf_out, "r+")
    fmt = read(attrs(io)["file_format"])
    (fmt == "ASDF") || (close(io); error("invalid ASDF file!"))
    if has(io, "Waveforms")
      wav = io["Waveforms"]
    else
      wav = g_create(io, "Waveforms")
    end
  else
    io = h5open(hdf_out, "cw")
    attrs(io)["file_format"] = "ASDF"
    attrs(io)["file_format_version"] = "1.0.2"
    wav = g_create(io, "Waveforms")
  end

  # =======================================================================
  # Add "empty" traces of all NaNs
  if add
    nc = length(chan_numbers)
    ts = zeros(Int64, nc)
    te = zeros(Int64, nc)
    for (i,j) in enumerate(chan_numbers)
      ts[i] = S.t[j][1,2]*1000
      te[i] = endtime(S.t[j], S.fs[j])*1000
    end
    p = convert(Nanosecond, len).value
    asdf_mktrace(S, xml_buf, chan_numbers, wav, ts, te, p, v, tag)
  end

  # write channels to net.sta waveform groups
  for j in 1:length(nsid)
    (v > 0) && println("writing ", nsid)
    id = nsid[j]
    chans = Int64[]
    for i in chan_numbers
      if (netsta[i] == id) && (S.fs[i] != 0.0)
        push!(chans, i)
      end
    end
    nc = length(chans)

    # does net.sta exist?
    if has(wav, id)
      sta = wav[id]

      if ovr
        # ====================================================================
        # overwrite StationXML
        (v > 1) && println("merging XML")

        # read old XML
        SX = SeisData()
        sxml = String(UInt8.(read(sta["StationXML"])))
        read_station_xml!(SX, sxml, "0001-01-01T00:00:00", "9999-12-31T23:59:59", true, v)

        # merge S headers into SX, overwriting SX
        SM = SeisData(length(chan_numbers))
        for f in (:id, :name, :loc, :fs, :gain, :resp, :units)
          setfield!(SM, f, deepcopy(getindex(getfield(S, f), chan_numbers)))
        end
        sxml_mergehdr!(SX, SM, false, true, v)

        # remake channel list; SX ordering differs from S
        cc = Int64[]
        for i in 1:SX.n
          idx = split_id(SX.id[i])
          ns = idx[1]*"."*idx[2]
          if ns == id
            push!(cc, i)
          end
        end
        o_delete(sta, "StationXML")
        asdf_wsxml(xml_buf, SX, cc, sta)

        # ==================================================================
        # overwrite trace data
        trace_ids = S.id[chans]
        t = Array{Array{Int64,2},1}(undef,nc)
        for (i,j) in enumerate(chans)
          t[i] = t_win(S.t[j], S.fs[j])
          broadcast!(*, t[i], t[i], 1000)
        end
        (v > 2) && println("t = ", t)

        # loop over trace waveforms using id, start time, end time
        for n in names(sta)
          (n == "StationXML") && continue
          cha_id = String(split(n, "_", limit=2, keepempty=true)[1])
          x = sta[n]
          (v > 1) && println("checking ", n)
          nx = length(x)
          t0 = read(x["starttime"])
          fs = read(x["sampling_rate"])
          (v > 2) && println("t0 = ", t0, "; fs = ", fs, "; nx = ", nx)

          # convert fs to sampling interval in ns
          Δ = round(Int64, 1.0e9/fs)
          t1 = t0 + (nx-1)*Δ

          for i in 1:length(chans)
            (trace_ids[i] == cha_id) || continue
            nk = size(t[i], 1)

            for k in 1:nk
              (v > 2) && println("segment k = ", k, "/", nk)
              # overlap
              ts = t[i][k,1]
              te = t[i][k,2]
              if (ts ≤ t1) && (te ≥ t0)
                # set channel index j in S
                j = chans[i]
                (v > 2) && println("j = ", j)

                # check for fs mismatch
                trace_fs = S.fs[j]
                (trace_fs == fs) || (@warn(string("Can't write ", S.id[j], "; fs mismatch!")); continue)
                lx = div(te-ts, Δ)+1

                # determine start, end indices in x that are overwritten
                (v > 2) && println("ts = ", ts, "; te = ", te, "; t0 = ", t0, "; t1 = ", t1, "; nx = ", nx, "; lx = ", lx)
                i0, i1, t2 = get_trace_bounds(ts, te, t0, t1, Δ, nx)

                # determine start, end indices in X to copy
                si = get_trace_bound(t0, ts, Δ, lx)
                ei = si + min(i1-i0, lx-1)

                # overwrite
                if v > 2
                  println("writing ", si, ":", ei, " to ", i0, ":", i1)
                end
                save_data!(S.x[j], x, i0:i1, si:ei)
              end
            end
          end

        end
      else
        for i in chans
          asdf_write_chan(S, sta, i, cha[i], evid, v)
        end
      end

    elseif ovr == false
      # Create Waveforms/net.sta
      sta = g_create(wav, id)

      # Create Waveforms/net.sta/chan_str
      for i in chans
        asdf_write_chan(S, sta, i, isempty(tag) ? cha[i] : tag, evid, v)
      end

      # Write StationXML to sta
      asdf_wsxml(xml_buf, S, chans, sta)
    end
  end

  close(xml_buf)
  close(io)

  # Correct :t
  if ovr
    for i in chan_numbers
      δ = get(S.misc[i], "tc", 0)
      t0 = S.t[i][1,2]
      S.t[i][1,2] = t0+δ
      delete!(S.misc[i], "tc")
    end
  end
  return nothing
end


text:
using .Gadfly
using Colors
using Printf

const NeXLSpectrumStyle = style(
    background_color = nothing,
    panel_fill = RGB(253 / 255, 253 / 255, 241 / 255),
    grid_color = RGB(255 / 255, 223 / 255, 223 / 255),
    grid_color_focused = RGB(255 / 255, 200 / 255, 200 / 255),
    grid_line_style = :solid,
    major_label_color = RGB(32 / 255, 32 / 255, 32 / 255),
    major_label_font_size = 9pt,
    panel_stroke = RGB(32 / 255, 32 / 255, 32 / 255),
    plot_padding = [2pt],
    key_title_font_size = 9pt,
    key_position = :right, # :bottom
)

"""
    plot(
	    specs::Union{Spectrum...,AbstractVector{Spectrum{<:Real}}};
	    klms=[],
	    edges=[],
		escapes=[],
		coincidences=[],
	    autoklms = false,
	    xmin=0.0,
	    xmax=missing,
	    norm=:None,
	    yscale=1.05,
	    ytransform = identity,
		style=NeXLSpectrumStyle,
		palette=NeXLPalette
    )::Plot

Required:

	specs::AbstractVector{Spectrum};

Named:

    klms = [ Element &| CharXRay ]
	edges = [ Element &| AtomicSubShell ]
	escapes = [ CharXRay ],
	coincidences = [ CharXRay ]
	autoklms = false # Add KLMs based on elements in spectra
	xmin = 0.0 # Min energy (eV)
	xmax = missing # Max energy (eV) (defaults to max(:BeamEnergy))
	norm = NoScaling() | ScaleDoseWidth() | ScaleDose() | ScaleSum() | ScaleROISum() | ScalePeak() | (<: SpectrumScaling)()
	yscale = 1.05 # Fraction of max intensity for ymax over [max(lld,xmin):xmax]
	ytransform = identity | log10 | sqrt | ??? # How to transform the counts data before plotting
	style=NeXLSpectrumStyle (or another Gadfly.style)
	palette = NeXLPalette | Colorant[ ... ] # Colors for spectra...
    customlayers = Gadfly.Layer[] # Allows additional plot layers to be added

Plot a multiple spectra on a single plot using Gadfly.

    Gadfly.plot(
      ffr::FilterFitResult,
      roi::Union{Missing,UnitRange{Int}} = missing;
      palette = NeXLPalette,
      style = NeXLSpectrumStyle, xmax=missing)

Plot the spectrum, residual and k-ratios using Gadfly.

    Gadfly.plot(fr::FilteredReference)

Plot a filtered reference spectrum.
"""
Gadfly.plot( #
    specs::AbstractVector{Spectrum{<:Real}};
    klms = Union{Element,CharXRay}[],
    edges = AtomicSubShell[],
    escapes = CharXRay[],
    coincidences = CharXRay[],
    autoklms = false,
    xmin = 0.0,
    xmax = missing,
    norm = :None,
    yscale = 1.05,
    ytransform = identity,
    style = NeXLSpectrumStyle,
    palette = NeXLPalette,
)::Plot = plot( #
    specs...,
    klms = klms,
    edges = edges,
    escapes = escapes,
    coincidences = coincidences,
    autoklms = autoklms,
    xmin = xmin,
    xmax = xmax,
    norm = norm,
    yscale = yscale,
    ytransform = ytransform,
    style = style,
    palette = palette,
)

function Gadfly.plot(
    specs::Spectrum{<:Real}...;
    klms = Union{Element,CharXRay}[],
    edges = AtomicSubShell[],
    escapes = CharXRay[],
    coincidences = CharXRay[],
    autoklms = false,
    xmin = 0.0,
    xmax = missing,
    legend = true,
    norm = NoScaling(),
    yscale = 1.05,
    ytransform = identity,
    style = NeXLSpectrumStyle,
    palette = NeXLPalette,
    customlayers = Gadfly.Layer[],
    duanehunt = false,
)::Plot
    function klmLayer(specdata, cxrs::AbstractArray{CharXRay})
        d = Dict{Any,Vector{CharXRay}}()
        for cxr in cxrs
            d[(element(cxr), shell(cxr))] =
                push!(get(d, (element(cxr), shell(cxr)), []), cxr)
        end
        x, y, label = [], [], []
        for cs in values(d)
            br = brightest(cs)
            ich = maximum(
                get(specdata[i], channel(energy(br), specs[i]), -1.0) for i in eachindex(specs)
            )
            if ich > 0
                for c in cs
                    push!(x, energy(c))
                    push!(y, ytransform(ich * weight(c)))
                    push!(label, weight(c) > 0.1 ? "$(element(c).symbol)" : "")
                end
            end
        end
        return layer(
            x = x,
            y = y,
            label = label,
            Geom.hair,
            Geom.point,
            Geom.label(position = :above),
            Theme(default_color = colorant"gray"),
        )
    end
    function edgeLayer(maxI, ashs::AbstractArray{AtomicSubShell})
        d = Dict{Any,Vector{AtomicSubShell}}()
        for ash in ashs
            d[(element(ash), shell(ash))] =
                push!(get(d, (element(ash), shell(ash)), []), ash)
        end
        x, y, label = [], [], []
        for ass in values(d)
            br = ass[findmax(capacity.(ass))[2]]
            for ash in ass
                push!(x, energy(ash))
                push!(y, ytransform(maxI * capacity(ash) / capacity(br)))
                push!(label, "$(ash)")
            end
        end
        return layer(
            x = x,
            y = y,
            label = label,
            Geom.hair,
            Geom.label(position = :right),
            Theme(default_color = colorant"lightgray"),
        )
    end
    function siEscapeLayer(crxs)
        x, y, label = [], [], []
        for xrs in crxs
            eesc = energy(xrs) - energy(n"Si K-L3")
            if eesc > 0.0
                ich = maximum(
                    get(specdata[i], channel(eesc, specs[i]), 0.0) for i in eachindex(specs)
                )
                push!(x, eesc)
                push!(y, ytransform(ich))
                push!(label, "$(element(xrs).symbol)\nesc")
            end
        end
        return layer(
            x = x,
            y = y,
            label = label,
            Geom.hair,
            Geom.label(position = :above),
            Theme(default_color = colorant"black"),
        )
    end
    function sumPeaks(cxrs)
        x, y, label = [], [], []
        for (i, xrs1) in enumerate(cxrs)
            for xrs2 in cxrs[i:end] # place each pair once...
                eesc = energy(xrs1) + energy(xrs2)
                ich = maximum(
                    get(specdata[i], channel(eesc, specs[1]), 0.0) for i in eachindex(specs)
                )
                if ich > 0.0
                    push!(x, eesc)
                    push!(y, ytransform(ich))
                    push!(label, "$(element(xrs1).symbol)\n+\n$(element(xrs2).symbol)")
                end
            end
        end
        return layer(
            x = x,
            y = y,
            label = label,
            Geom.hair,
            Geom.label(position = :above),
            Theme(default_color = colorant"gray"),
        )
    end
    @assert length(specs) <= length(palette) "The palette must specify at least as many colors as spectra."
    specdata = [scaledcounts(norm, s) for s in specs]
    ylbl = repr(norm)
    maxI, maxE, maxE0 = 16, 1.0e3, 1.0e3
    names, layers = String[], Layer[]
    append!(layers, customlayers)
    for (i, spec) in enumerate(specs)
        mE =
            ismissing(xmax) ? get(spec, :BeamEnergy, energy(length(spec), spec)) :
            convert(Float64, xmax)
        mE0 = get(spec, :BeamEnergy, missing)
        chs =
            max(
                1,
                channel(convert(Float64, xmin), spec),
            ):min(length(spec), channel(mE, spec))
        mchs =
            max(channel(200.0, spec), chs[1], lld(spec)):min(length(specdata[i]), chs[end])  # Ignore zero strobe...
        maxI = max(maxI, maximum(specdata[i][mchs]))
        maxE = max(maxE, mE)
        maxE0 = ismissing(mE0) ? maxE : max(maxE, mE0)
        push!(names, spec[:Name])
        ly = Gadfly.layer(
            x = energyscale(spec)[chs],
            y = ytransform.(specdata[i][chs]), #
            Geom.step,
            Theme(default_color = palette[i]),
        )
        append!(layers, ly)
    end
    if duanehunt
        dhx, dhy = Float64[], Float64[]
        for i in eachindex(specs)
            append!(dhx, duane_hunt(specs[i]))
            append!(dhy, 0.1*maxI*yscale)
        end
        append!(layers, layer(x=dhx, y=dhy, color=palette[eachindex(specs)], Geom.hair(orientation=:vertical), Geom.point))
    end
    autoklms && append!(klms, mapreduce(s -> elms(s, true), union!, specs))
    if length(klms) > 0
        tr(elm::Element) = characteristic(elm, alltransitions, 1.0e-3, maxE0)
        tr(cxr::CharXRay) = [cxr]
        pklms = mapreduce(klm -> tr(klm), append!, klms)
        if length(pklms) > 0
            append!(layers, klmLayer(specdata, pklms))
        end
    end
    if length(edges) > 0
        shs(elm::Element) = atomicsubshells(elm, maxE0)
        shs(ash::AtomicSubShell) = [ash]
        pedges = mapreduce(ash -> shs(ash), append!, edges)
        if length(pedges) > 0
            append!(layers, edgeLayer(0.5 * maxI, pedges))
        end
    end
    if length(escapes) > 0
        append!(layers, siEscapeLayer(escapes))
    end
    if length(coincidences) > 0
        append!(layers, sumPeaks(coincidences))
    end
    Gadfly.with_theme(style) do
        leg =
            legend ?
            tuple(
                Guide.manual_color_key(
                    length(specs) > 1 ? "Spectra" : "Spectrum",
                    names,
                    palette[1:length(specs)],
                ),
            ) : tuple()
        try
            plot(
                layers...,
                Guide.XLabel("Energy (eV)"),
                Guide.YLabel(ylbl),
                Scale.x_continuous(format = :plain),
                Scale.y_continuous(format = :plain),
                Coord.Cartesian(
                    ymin = 0,
                    ymax = ytransform(yscale * maxI),
                    xmin = convert(Float64, xmin),
                    xmax = maxE,
                ),
                leg...,
            )
        catch
            plot(
                layers...,
                Guide.XLabel("Energy (eV)"),
                Guide.YLabel(ylbl),
                Scale.x_continuous(format = :plain),
                Scale.y_continuous(format = :plain),
                Coord.Cartesian(
                    ymin = 0,
                    ymax = ytransform(yscale * maxI),
                    xmin = convert(Float64, xmin),
                    xmax = maxE,
                ),
                leg...,
            )
        end
    end
end

function Gadfly.plot(
    ffr::FilterFitResult,
    roi::Union{Nothing,AbstractUnitRange{<:Integer}} = nothing;
    palette = NeXLPalette,
    style = NeXLSpectrumStyle,
    xmax::Union{AbstractFloat, Nothing} = nothing,
    comp::Union{Material, Nothing} = nothing,
    det::Union{EDSDetector, Nothing} = nothing,
    resp::Union{AbstractArray{<:AbstractFloat,2},Nothing} = nothing,
    yscale = 1.0
)
    function defroi(ffrr) # Compute a reasonable default display ROI
        tmp =
            minimum(
                lbl.roi[1] for lbl in keys(ffrr.kratios)
            ):maximum(lbl.roi[end] for lbl in keys(ffrr.kratios))
        return max(
            lld(ffr.label.spectrum),
            tmp[1] - length(ffrr.roi) ÷ 40,
        ):min(tmp[end] + length(ffrr.roi) ÷ 10, ffrr.roi[end])
    end
    roilt(l1, l2) = isless(l1.roi[1], l2.roi[1])
    roi = something(roi, defroi(ffr))
    layers = [
        layer(x = roi, y = ffr.residual[roi], Geom.step, Theme(default_color = palette[2])),
        layer(x = roi, y = ffr.raw[roi], Geom.step, Theme(default_color = palette[1])),
    ]
    # If the information is available,also model the continuum
    comp = isnothing(comp) ? get(spectrum(ffr), :Composition, nothing) : comp
    det = isnothing(det) ? get(spectrum(ffr), :Detector, nothing) : det
    if !any(isnothing.((comp, resp, det)))
        cc = fitcontinuum(spectrum(ffr), det, resp)
        push!(layers, layer(x=roi, y=cc[roi], Geom.line, Theme(default_color=palette[2])))
    end
    miny, maxy, prev, i =
        minimum(ffr.residual[roi]), 3.0 * yscale * maximum(ffr.residual[roi]), -1000, -1
    for lbl in sort(collect(keys(ffr.kratios)), lt = roilt)
        if NeXLUncertainties.value(ffr, lbl) > 0.0
            # This logic keeps the labels on different lines (mostly...)
            i, prev =
                (lbl.roi[1] > prev + length(roi) ÷ 10) || (i == 6) ? (0, lbl.roi[end]) :
                (i + 1, prev)
            labels = ["", name(lbl.xrays)]
            # Plot the ROI
            push!(
                layers,
                layer(
                    x = [lbl.roi[1], lbl.roi[end]],
                    y = maxy * [0.4 + 0.1 * i, 0.4 + 0.1 * i],
                    label = labels,
                    Geom.line,
                    Geom.point,
                    Geom.label(position = :right),
                    Theme(default_color = "gray"),
                ),
            )
            # Plot the k-ratio as a label above ROI
            push!(
                layers,
                layer(
                    x = [0.5 * (lbl.roi[1] + lbl.roi[end])],
                    y = maxy * [0.4 + 0.1 * i],
                    label = [@sprintf("%1.4f", NeXLUncertainties.value(ffr, lbl))],
                    Geom.label(position = :above),
                    Theme(default_color = "gray"),
                ),
            )
        end
    end
    Gadfly.with_theme(style) do
        plot(
            layers...,
            Coord.cartesian(
                xmin = roi[1],
                xmax = something(xmax, roi[end]),
                ymin = min(1.1 * miny, 0.0),
                ymax = maxy,
            ),
            Guide.XLabel("Channels"),
            Guide.YLabel("Counts"),
            Guide.title("$(ffr.label)"),
        )
    end
end

function Gadfly.plot(fr::FilteredReference; palette = NeXLPalette)
    roicolors = Colorant[RGB(0.9, 1.0, 0.9), RGB(0.95, 0.95, 1.0)]
    layers = [
        layer(x = fr.ffroi, y = fr.data, Theme(default_color = palette[1]), Geom.step),
        layer(x = fr.ffroi, y = fr.filtered, Theme(default_color = palette[2]), Geom.step),
        layer(x = fr.roi, y = fr.charonly, Theme(default_color = palette[3]), Geom.step),
        layer(
            xmin = [fr.ffroi[1], fr.roi[1]],
            xmax = [fr.ffroi[end], fr.roi[end]],
            Geom.vband,
            color = roicolors,
        ),
    ]
    try
        plot(
            layers...,
            Coord.cartesian(
                xmin = fr.ffroi[1] - length(fr.ffroi) ÷ 10,
                xmax = fr.ffroi[end] + length(fr.ffroi) ÷ 10,
            ),
            Guide.xlabel("Channel"),
            Guide.ylabel("Counts"),
            Guide.title(repr(fr.label)),
            Guide.manual_color_key(
                "Legend",
                ["Spectrum", "Filtered", "Char. Only", "Filter ROC", "Base ROC"],
                [palette[1:3]..., roicolors...],
            ),
        )
    catch
        plot(
            layers...,
            Coord.cartesian(
                xmin = fr.ffroi[1] - length(fr.ffroi) ÷ 10,
                xmax = fr.ffroi[end] + length(fr.ffroi) ÷ 10,
            ),
            Guide.xlabel("Channel"),
            Guide.ylabel("Counts"),
            Guide.title(repr(fr.label)),
            #    Guide.manual_color_key("Legend",["Spectrum", "Filtered", "Char. Only", "Filter ROC", "Base ROC"], [ palette[1:3]..., roicolors...] )
        )
    end
end

Gadfly.plot(ff::TopHatFilter, fr::FilteredReference) = spy(
    filterdata(ff, fr.ffroi),
    Guide.title(repr(fr.label)),
    Guide.xlabel("Channels"),
    Guide.ylabel("Channels"),
)

function Gadfly.plot(vq::VectorQuant, chs::UnitRange)
    colors = distinguishable_colors(
        size(vq.vectors, 1) + 2,
        [
            RGB(253 / 255, 253 / 255, 241 / 255),
            RGB(0, 0, 0),
            RGB(0 / 255, 168 / 255, 45 / 255),
        ],
        transform = deuteranopic,
    )[3:end]
    lyrs = mapreduce(
        i -> layer(
            x = chs,
            y = vq.vectors[i, chs],
            Theme(default_color = colors[i]),
            Geom.line,
        ),
        append!,
        eachindex(vq.references),
    )
    try
        plot(
            lyrs...,
            Guide.xlabel("Channel"),
            Guide.ylabel("Filtered"),
            Guide.manual_color_key(
                "Vector",
                [repr(r[1]) for r in vq.references],
                color = Colorant[colors...],
            ),
        )
    catch
        plot(
            lyrs...,
            Guide.xlabel("Channel"),
            Guide.ylabel("Filtered"),
            #    Guide.manual_color_key("Vector", [ repr(r[1]) for r in vq.references ], color=Colorant[colors...])
        )
    end
end

function Gadfly.plot(deteff::DetectorEfficiency, emax = 20.0e3)
    eff(ee) = efficiency(deteff, ee, π / 2)
    plot(eff, 100.0, emax)
end

function plotandimage(plot::Gadfly.Plot, image::Array)
    io = IOBuffer(maxsize = 10 * 1024 * 1024)
    save(Stream(format"PNG", io), image)
    pix = max(size(image, 1), size(image, 2))
    scaleX, scaleY = size(image, 1) / pix, size(image, 2) / pix
    return compose(
        context(),
        (context(0.0, 0.0, 0.8, 1.0), Gadfly.render(plot)),
        (
            context(0.8, 0.0, 0.2, 1.0),
            bitmap(
                "image/png",
                take!(io),
                0.5 * (1.0 - scaleX),
                0.5 * (1.0 - scaleY),
                scaleX,
                scaleY,
            ),
        ),
    )
end

Gadfly.plot(ffp::FilterFitPacket; kwargs...) =
    plot(unique(spectra(ffp))...; klms = collect(elms(ffp)), kwargs...)

"""
    plot_compare(specs::AbstractArray{<:Spectrum}, mode=:Plot; xmin=100.0, xmax=1.0, palette = NeXLPalette)

Plots a comparison of the channel-by-channel data from each individual spectrum relative to the dose-corrected
mean of the other spectra.  Count statistics are taken into account so if the spectra agree to within count
statistics we expect a mean of 0.0 and a standard deviation of 1.0 over all channels. Note: xmax is relative
to the :BeamEnergy.
"""
function plot_compare(specs::AbstractArray{<:Spectrum}, mode=:Plot; xmin=100.0, xmax=1.0, palette = NeXLPalette)
    channels(spec) = channel(100.0, spec):channel(xmax*get(spec,:BeamEnergy,20.0e3),spec)
    if mode==:Plot
        layers = [
            layer(x=energyscale(specs[i],channels(specs[i])), y = sigma(specs[i], specs, channels(specs[i])),
                Theme(default_color = palette[i], alphas=[0.4]))
                    for i in eachindex(specs)
        ]
        plot(layers..., Guide.xlabel("Energy (eV)"), Guide.ylabel("σ"),
        Guide.manual_color_key("Spectra", [String(spec[:Name]) for spec in specs], palette[eachindex(specs)]),
        Coord.cartesian(xmin=100.0, xmax=xmax*maximum(get(spec, :BeamEnergy, 20.0e3) for spec in specs))
        )
    elseif mode==:Histogram
        layers = [
            layer(x = sigma(specs[i], specs, channels(specs[i])), Geom.histogram(),
                Theme(default_color = palette[i], alphas=[0.2]))
                    for i in eachindex(specs)
        ]
        plot(layers..., Guide.xlabel("σ"), Guide.manual_color_key("Spectra", [String(spec[:Name]) for spec in specs], palette[eachindex(specs)]))
    else
        error("Unknown mode in plot_compare(...):  Not :Plot or :Histogram.")
    end
end

"""
    plot_multicompare(specs::AbstractArray{Spectrum{T}}; minE=200.0, maxE=0.5*specs[1][:BeamEnergy]) where { T<: Real}

Compare spectra collected simultaneously on multiple detectors in a single acquisition.
"""
function plot_multicompare(specs::AbstractArray{Spectrum{T}}; minE=200.0, maxE=0.5*specs[1][:BeamEnergy]) where { T<: Real}
    s, mcs = specs[1], multicompare(specs)
    chs = max(1,channel(minE, s)): min(channel(maxE, s), length(s))
    xx = map(i->energy(i,s), chs)
    plot(
        (layer(x=xx, y=view(mc,chs), Geom.line, Theme(default_color=c)) for (c,mc) in zip(NeXLPalette[1:length(mcs)],mcs))...,
        Guide.xlabel("Energy (eV)"), Guide.ylabel("Ratio")
    )
end

@info "Loading Gadfly support into NeXLSpectrum."


text:
abstract type Algorithm end
# Base.copy(o::Algorithm) = deepcopy(o)
init!(o::Algorithm, p) = o
update!(o::Algorithm, gx) = nothing
Base.merge!(o::T, o2::T, γ) where {T<:Algorithm} = o
Base.show(io::IO, o::Algorithm) = print(io, name(o, false, false))

abstract type SGAlgorithm <: Algorithm end

#-----------------------------------------------------------------------# SGD
"""
    SGD()

Stochastic Gradient Descent.
"""
struct SGD <: SGAlgorithm end

# #-----------------------------------------------------------------------# NSGD
# """
#     NSGD(α)

# Nesterov accelerated Proximal Stochastic Gradient Descent.
# """
# struct NSGD <: SGUpdater
#     α::Float64
#     v::Vector{Float64}
#     θ::Vector{Float64}
#     NSGD(α = 0.0, p = 0) = new(α, zeros(p), zeros(p))
# end
# init(u::NSGD, p::Int) = NSGD(u.α, p)
# function Base.merge!(o::NSGD, o2::NSGD, γ::Float64)
#     o.α == o2.α || error("Merge Failed.  NSGD objects use different α.")
#     smooth!(o.v, o2.v, γ)
#     smooth!(o.θ, o2.θ, γ)
# end

#-----------------------------------------------------------------------# ADAGRAD
"""
    ADAGRAD()

A variation of [`SGD`](@ref) with element-wise weights generated by the average of the squared gradients.
"""
mutable struct ADAGRAD <: SGAlgorithm
    h::Vector{Float64}
    n::Int
    ADAGRAD() = new(zeros(0), 0)
end
init!(o::ADAGRAD, p) = (o.h = zeros(p))
Base.merge!(o::ADAGRAD, o2::ADAGRAD, γ) = (smooth!(o.h, o2.h, γ); o)
function update!(o::ADAGRAD, gx)
    o.n += 1
    for j in eachindex(o.h)
        o.h[j] = smooth(o.h[j], gx[j]^2, 1 / o.n)
    end
end

#-----------------------------------------------------------------------# RMSPROP
"""
    RMSPROP(α = .9)

A Variation of [`ADAGRAD`](@ref) that uses element-wise weights generated by an exponentially
weighted mean of the squared gradients.
"""
mutable struct RMSPROP <: SGAlgorithm
    α::Float64
    h::Vector{Float64}
    RMSPROP(α = .9) = new(α, zeros(0))
end
init!(a::RMSPROP, p) = (a.h = zeros(p))
function Base.merge!(o::RMSPROP, o2::RMSPROP, γ)
    o.α = smooth(o.α, o2.α, γ)
    smooth!(o.h, o2.h, γ)
end
function update!(o::RMSPROP, gx)
    for j in eachindex(o.h)
        o.h[j] = smooth(gx[j]^2, o.h[j], o.α)
    end
end

#-----------------------------------------------------------------------# ADADELTA
"""
    ADADELTA(ρ = .95)

An extension of [`ADAGRAD`](@ref).
"""
mutable struct ADADELTA <: SGAlgorithm
    v::Vector{Float64}
    Δ::Vector{Float64}
    δ::Vector{Float64}
    ρ::Float64
    ADADELTA(ρ = .95) = new(Float64[], Float64[], Float64[], ρ)
end
init!(o::ADADELTA, p) = (o.v = zeros(p); o.δ = zeros(p); o.Δ = zeros(p))
function Base.merge!(o::ADADELTA, o2::ADADELTA, γ)
    smooth!(o.v, o2.v, γ)
    smooth!(o.Δ, o2.Δ, γ)
    o.ρ = smooth(o.ρ, o2.ρ, γ)
    o
end
function update!(o::ADADELTA, gx)
    for j in eachindex(o.v)
        g2 = gx[j] ^ 2
        o.v[j] = smooth(g2 + ϵ, o.v[j], o.ρ)
        o.δ[j] = sqrt(o.Δ[j] / o.v[j] + ϵ)  # not multiplied by gx[j] as in paper
        o.Δ[j] = smooth(o.δ[j]^2 * g2, o.Δ[j], o.ρ)
    end
end

#-----------------------------------------------------------------------# ADAM
"""
    ADAM(β1 = .99, β2 = .999)

A variant of [`SGD`](@ref) with element-wise learning rates generated by exponentially weighted
first and second moments of the gradient.
"""
mutable struct ADAM <: SGAlgorithm
    m::Vector{Float64}
    v::Vector{Float64}
    β1::Float64
    β2::Float64
    ADAM(β1 = .99, β2 = .999) = new(Float64[], Float64[], β1, β2)
end
init!(o::ADAM, p) = (o.m = zeros(p); o.v = zeros(p))
function Base.merge!(o::ADAM, o2::ADAM, γ)
    smooth!(o.m, o2.m, γ)
    smooth!(o.v, o2.v, γ)
    o.β1 = smooth(o.β1, o2.β1, γ)
    o.β2 = smooth(o.β2, o2.β2, γ)
    o
end
function update!(o::ADAM, gx)
    for j in eachindex(o.m)
        g = gx[j]
        o.m[j] = smooth(g, o.m[j], o.β1)
        o.v[j] = smooth(g * g, o.v[j], o.β2)
    end
end

#-----------------------------------------------------------------------# ADAMAX
"""
    ADAMAX(η, β1 = .9, β2 = .999)

ADAMAX with momentum parameters `β1`, `β2`.  ADAMAX is an extension of [`ADAM`](@ref).
"""
mutable struct ADAMAX <: SGAlgorithm
    β1::Float64
    β2::Float64
    m::Vector{Float64}
    v::Vector{Float64}
    function ADAMAX(β1::Float64 = 0.9, β2::Float64 = .999)
        0 < β1 < 1 || error("β1 must be in (0, 1)")
        0 < β2 < 1 || error("β2 must be in (0, 1)")
        new(β1, β2, zeros(0), zeros(0))
    end
end
init!(o::ADAMAX, p) = (o.m = zeros(p); o.v = zeros(p))
function Base.merge!(o::ADAMAX, o2::ADAMAX, γ)
    smooth!(o.m, o2.m, γ)
    smooth!(o.v, o2.v, γ)
    o.β1 = smooth(o.β1, o2.β1, γ)
    o.β2 = smooth(o.β2, o2.β2, γ)
    o
end
function update!(o::ADAMAX, gx)
    for j in eachindex(o.m)
        o.m[j] = smooth(gx[j], o.m[j], o.β1)
        o.v[j] = max(abs(gx[j]) + ϵ, o.β2 * o.v[j])
    end
end

# #-----------------------------------------------------------------------# NADAM
# """
#     NADAM(α1 = .99, α2 = .999)

# Adaptive Moment Estimation with momentum parameters `α1` and `α2`.
# """
# mutable struct NADAM <: SGUpdater
#     β1::Float64
#     β2::Float64
#     M::Vector{Float64}
#     V::Vector{Float64}
#     nups::Int
#     function NADAM(β1::Float64 = 0.99, β2::Float64 = .999, p::Integer = 0)
#         0 < β1 < 1 || error("β1 must be in (0, 1)")
#         0 < β2 < 1 || error("β2 must be in (0, 1)")
#         new(β1, β2, zeros(p), zeros(p), 0)
#     end
# end
# init(u::NADAM, p) = NADAM(u.β1, u.β2, p)
# function Base.merge!(o::NADAM, o2::NADAM, γ::Float64)
#     (o.β1 == o2.β1) && (o.β2 == o2.β2) ||
#         error("Merge failed.  NADAM objects use different momentum parameters.")
#     o.nups += o2.nups
#     smooth!(o.M, o2.M, γ)
#     smooth!(o.V, o2.V, γ)
# end

#-----------------------------------------------------------------------# MSPI, OMAS, OMAP
"""
    MSPI()

Majorized Stochastic Proximal Iteration.
"""
struct MSPI <: Algorithm end

"""
    OMAP()

Online MM via Averaged Parameter.
"""
struct OMAP <: Algorithm end

"""
    OMAS()

Online MM via Averaged Surrogate.
"""
mutable struct OMAS <: Algorithm
    a::Vector{Float64}
    b::Vector{Float64}
    A::Matrix{Float64}
    B::Matrix{Float64}
    OMAS() = new(zeros(0), zeros(0), zeros(0, 0), zeros(0, 0))
end
init!(o::OMAS, p) = (o.a = zeros(p); o.b = zeros(p))


text:
# basic pin hole camera model

#TODO: create a package that defines a camera model, there are lots of different ones all over the place
"""
	$(TYPEDEF)
Data structure for a Camera model with parameters.
Use `CameraModel(width,height,fc,cc,skew,kc)` for easy construction.
"""
struct CameraModel
    width::Int		# image width
    height::Int		# image height
    fc::Vector{Float64}	# focal length in x and y
	f::Float64			# if one focal lenght is used
    cc::Vector{Float64}	# camera center
    skew::Float64	    # skew value
    kc::Vector{Float64} # distortion coefficients up to fifth order or NTuple{5,Float64}
    K::Matrix{Float64} # 3x3 camera calibration matrix (Camera intrinsics)
    Ki::Matrix{Float64} # inverse of a 3x3 camera calibratio matrix
end

"""
	$(SIGNATURES)
Constructor helper for creating a camera model.
"""
function CameraModel(width::Int,height::Int,fc::Vector{Float64},cc::Vector{Float64},skew::Float64,kc::Vector{Float64})::CameraModel
    KK = [fc[1]      skew  cc[1];
             0       fc[2] cc[2];
             0		    0     1]
    Ki = inv(KK)
    CameraModel(width,height,fc, mean(fc),cc,skew,kc,KK,Ki)
end

"""
	$(SIGNATURES)
Constructor helper for creating a camera model.
"""
function CameraModel(width::Int,height::Int,f::Float64,cc::Vector{Float64})::CameraModel
    KK = [f 0 cc[1];
          0 f cc[2];
          0	0     1.]
    Ki = inv(KK)
    return CameraModel(width,height,[f,f], f, cc, 0.0, [0.], KK,Ki)
end

mutable struct DepthCamera
  K::Array{Float64,2}
  shape::Tuple{Int, Int}
  skip::Int
  D::Vector{Float64}
  xs::Array{Float64}
  ys::Array{Float64}
  DepthCamera(K::Array{Float64,2};
      shape::Tuple{Int,Int}=(480,640),
      skip::Int=1,
      D::Vector{Float64}=zeros(5) ) = new( K, shape, skip, D, Array{Float64,2}(), Array{Float64,2}() )
end


text:
"""
    get_svg(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String

Get a SVG depiction of the mol.

```julia
svg = get_svg(mol)
```
"""
function get_svg(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String
    details_json::String = jsonify_details(details)
    val::Cstring = ccall((:get_svg, librdkitcffi), Cstring, (Cstring, Csize_t, Cstring), mol.mol[], mol.mol_size[], details_json)
    svg = unsafe_string_and_free(val)
    return svg
end


text:
function readdata(path::String, worddict::Dict, chardict::Dict, tagdict::Dict)
    unkwordid = worddict["UNKNOWN"]
    data_w, data_c, data_t = Var[], Var[], Var[]

    lines = open(readlines, path)
    i = 1
    while i <= length(lines)
        j = i
        while j < length(lines) && !isempty(lines[j+1])
            j += 1
        end
        w, cs, t = Int[], Vector{Int}[], Int[]
        for k = i:j
            items = split(lines[k], '\t')
            word, tag = items[2], items[5]
            word0 = replace(word, r"[0-9]", '0')
            wordid = get(worddict, lowercase(word0), unkwordid)
            chars = Vector{Char}(word0)
            charids = map(c -> get!(chardict,c,length(chardict)+1), chars)
            tagid = get!(tagdict, tag, length(tagdict)+1)
            push!(w, wordid)
            push!(cs, charids)
            push!(t, tagid)
        end
        push!(data_w, Var(w))
        batchdims = map(length, cs)
        c = cat(1, cs...)
        push!(data_c, Var(c,batchdims))
        push!(data_t, Var(t))
        i = j + 2
    end
    data_w, data_c, data_t
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
# ====================================================================
# Component structure
mutable struct OffsetSlope_1D <: AbstractComponent
    offset::Parameter
    x0::Parameter
    slope::Parameter

    function OffsetSlope_1D(offset::Number, x0::Number, slope::Number)
        out = new(Parameter(offset), Parameter(x0), Parameter(slope))
        out.x0.free = false
        return out
    end
end

mutable struct OffsetSlope_2D <: AbstractComponent
    offset::Parameter
    x0::Parameter
    y0::Parameter
    slopeX::Parameter
    slopeY::Parameter

    function OffsetSlope_2D(offset::Number, x0::Number, y0::Number, slopeX::Number, slopeY::Number)
        out = new(Parameter(offset), Parameter(x0), Parameter(y0), Parameter(slopeX), Parameter(slopeY))
        out.x0.free = false
        out.y0.free = false
        return out
    end
end

OffsetSlope(offset, x0, slope) = OffsetSlope_1D(offset, x0, slope)
OffsetSlope(offset, x0, y0, slopeX, slopeY) = OffsetSlope_2D(offset, x0, y0, slopeX, slopeY)

ceval_data(domain::Domain_1D, comp::OffsetSlope_1D) = (nothing, length(domain))
ceval_data(domain::Domain_2D, comp::OffsetSlope_2D) = (nothing, length(domain))

# ====================================================================
# Evaluate component
function evaluate(c::CompEval{Domain_1D, OffsetSlope_1D},
                  offset, x0, slope)
    @. (c.eval = slope * (c.domain[1] - x0) + offset)
end

function evaluate(c::CompEval{Domain_2D, OffsetSlope_2D},
                   offset, x0, y0, slopeX, slopeY)
    x = c.domain[1]
    y = c.domain[2]
    @. (c.eval =
        slopeX * (x - x0) +
        slopeY * (y - y0) +
        offset)
end


text:
immutable LevenshteinDistance <: Distance
end

"""
The Levenshtein distance between two sequences is the minimum number of
single element edits, which may be insertions, deletions, or substitutions,
that are required to transform one sequence into the other.
"""
levenshtein() = LevenshteinDistance()

dist(::LevenshteinDistance, x::DirectIndexString, y::DirectIndexString) =
  error("Not implemented!")


text: using Pkg; pkg"registry add General https://github.com/legend-exp/LegendJuliaRegistry.git"


text:
using BioinformaticsBISC195
using Test

@testset "BioinformaticsBISC195" begin

@testset "Using Strings" begin

    @testset "normalizeDNA" begin
        @test normalizeDNA("aatgn") == "AATGN"
        @test_throws Exception normalizeDNA("ZCA")
        @test_throws Exception normalizeDNA(42)
        c = normalizeDNA('C')
        @test c == "C"
        @test typeof(c) == String
    end # normalizeDNA

    @testset "composition" begin
        seq = rand(['A','T','G','C','N'], 20) |> join
        bc = composition(seq)
        @test bc isa Tuple

        @test bc[1] == count(x-> x == 'A', seq)
        @test bc[3] == count(x-> x == 'C', seq)
        @test bc[2] == count(x-> x == 'G', seq)
        @test bc[4] == count(x-> x == 'T', seq)
        @test bc[5] == count(x-> x == 'N', seq)

        bc = composition(lowercase(seq))

        @test bc[1] == count(x-> x == 'A', seq)
        @test bc[3] == count(x-> x == 'C', seq)
        @test bc[2] == count(x-> x == 'G', seq)
        @test bc[4] == count(x-> x == 'T', seq)
        @test bc[5] == count(x-> x == 'N', seq)
    end # composition

    @testset "gc_content" begin
        @test gc_content("ANTG") == 0.25
        @test gc_content("cccggg") * 100 == 100.0
        @test gc_content("ATta") == 0.0
        @test_throws Exception gc_content("ATtz")
    end # gc_content

    @testset "complement" begin
        @test complement("ATTAN") == "TAATN"
        @test complement("gcta") == "CGAT"
        @test complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception complement("AZC")
    end # complement

    @testset "reverse_complement" begin
        @test reverse_complement("ATTAN") == "NTAAT"
        @test reverse_complement("gcta") == "TAGC"
        @test reverse_complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception reverse_complement("AZC")
    end # reverse_complement

    @testset "parse_fasta" begin
        testpath = normpath(joinpath(@__DIR__, "..", "data"))
        genomes = joinpath(testpath, "cov2_genomes.fasta")
        ex1_path = joinpath(testpath, "ex1.fasta")
        ex2_path = joinpath(testpath, "ex2.fasta")

        ex1 = parse_fasta(ex1_path)
        @test ex1 isa Tuple
        @test all(x-> x isa String, ex1[1])
        @test all(x-> x isa String, ex1[2])

        @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
        @test ex1[2] == ["AATTATAGC", "CGCCCCCCAGTCGGATT"]

        @test_throws Exception parse_fasta(ex2_path)

        cov2 = parse_fasta(genomes)
        @test length(cov2[1]) == 8
        @test length(cov2[2]) == 8
    end #parse_fasta

    @testset "removeUnder25k" begin
        parse1 = parse_fasta("../data/test1.fasta")
        parse2 = parse_fasta("../data/test2.fasta")
        cleaned1 = removeUnder25k(parse1)
        cleaned2 = removeUnder25k(parse2)
        @test length(findall(seq -> length(seq)<25000, cleaned2[2])) == 0
        @test length(cleaned1) ==length(parse1[2])
        @test length(cleaned2[1]) == length(cleaned2[2])
    end

    @testset "uniqueKmers" begin
        testex1 = uniqueKmers("AGTCTCTCTCGGGGTA", 4)
        @test length(testex1) == 10
        @test testex1 isa Vector
        @test testex1[end] isa String
        @test testex1[end] == "GGTA"
        @test length(testex1[3]) == 4

        testex2 = uniqueKmers("ATTTATATAGCTNNA", 3)
        @test length(testex2) == 8
        @test length(testex2[1]) == 3
        @test testex2[begin] == "ATT"
        @test testex2[3] == "TTA"
        @test testex2 isa Vector
    end

    @testset "compareKmerSets" begin
        testex1 = uniqueKmers("AGTCTCTCTCGGGGTA", 4)
        testex2 = uniqueKmers("ATTTATATAGCTNNA", 3)
        @test compareKmerSets(testex1, testex2) == 1.0
        @test compareKmerSets(testex1, testex2) isa Float64

        testex3 = ["ATT", "AGT", "CTC", "AGC"]
        testex4 = ["ATT", "CTC", "AGG", "AGA"]
        @test compareKmerSets(testex3, testex4) == 0.6666666666666667
    end

    @testset "nwscore" begin
        @test_throws ArgumentError nwscore(nothing, nothing)
        @test nwscore('A', 'G') == -1
        @test nwscore('A', 'A'; match = 10) == 10
        @test nwscore('A', 'T'; mismatch = -4) !=-1
        @test_throws MethodError nwscore("A", "C")
    end

    @testset "nwsetupmatrix" begin
        test5 = nwsetupmatrix("AGAGAGACT", "AGACAGCCT")
        @test test5 isa Matrix
        nrow, ncol = size(test5)
        @test nrow == ncol
        @test nrow == 10
        @test test5[1:10] isa Vector
        @test test5[1:10][10] == -9.0
        test6 = nwsetupmatrix("AGCT", "AGTTGAGAGTAT"; gap = -3)
        nrow2, ncol2 = size(test6)
        @test test6[1, 13] == -3 * (length("AGTTGAGAGTAT"))
        @test test6[5] == -12.0
        @test nrow2 == length("AGCT") + 1
        @test ncol2 == length("AGTTGAGAGTAT") + 1
    end

    @testset "nwscorematrix" begin
        ex7 = nwscorematrix("AGCT", "AGTTGAGAGTAT")
        @test ex7 isa Dict
        @test ex7[(4, 6)]==[-1.0, "left"]
        @test ex7[(2, 8)]==[-5.0, "left"]
        @test ex7[(3, 3)][1] isa Float64
        @test ex7[(3, 3)][2] isa String
        @test ex7[(5, 12)][2] == "left"
        @test collect(keys(ex7))[1] isa Tuple
        @test collect(keys(ex7))[2] == (5, 5)
        @test length(ex7) == 48
    end

    @testset "nwalign" begin
        ex8 = nwalign("AGCT", "AGCCA"; match = 2, mismatch = -3)
        @test ex8 isa Tuple
        @test ex8 == ("AGC--T", "AGCCA-")
        @test nwalign("AGCT", "AGCCA")[1] == "AGCT-"
        @test nwalign("AGCT", "AGCCA")[2] == "AGCCA"

        ex9 = nwalign("GTGCCGA", "GAACCGTA"; gap=-2)
        @test ex9 == ("GTGCCG-A", "GAACCGTA")
    end

    @testset "alignmentHeatmap" begin
        ex10 = alignmentHeatmap("GTA", "ATA")
        nrow, ncol = size(ex10)
        @test nrow == ncol == 3
        @test ex10[3, 1] == 1.0
        @test ex10[2, 2] == 1.0
        @test ex10[3, 3] == 1.0

        ex11 = alignmentHeatmap("ATAAA", "ATAAA")
        @test ex11[1, 1] == ex11[2, 2] == ex11[3, 3] == ex11[4, 4] == ex11[5, 5]== 1.0
    end

    @testset "mis_matchSeq" begin
        aligned1 = nwalign("ATATAGCAAAAAA", "ATAGGCTATAAAA")
        ex12 = mis_matchSeq(aligned1)
        @test ex12 isa Dict
        @test ex12["Match#"] isa Int64
        @test ex12["Mismatch#"] isa Int64
        @test ex12["Mis/Match%"] isa Tuple
        @test ex12["Match#"] == 10
        @test ex12["Mismatch#"] == 4
        @test round(ex12["Mis/Match%"][1], digits=6) == 0.714286
    end

    @testset "mis_matchSubsequences" begin
        aligned1 = nwalign("ATATAGCAAAAAA", "ATAGGCTATAAAA")
        ex14 = mis_matchSubsequences(aligned1)
        @test ex14 isa Vector
        @test ex14[1].type == "match"
        @test length(ex14[1].sequences) == 1
        @test length(ex14[7].sequences[1]) == 4
        @test ex14[4].positions isa UnitRange
        @test ex14[4].positions == 8:8
        @test length(ex14[4].sequences) == 2
    end
end # strings

# @testset "Using BioSequences" begin

#     @testset "normalizeDNA" begin
#         @test normalizeDNA("aatgn") == dna"AATGN"
#         @test_throws Exception normalizeDNA("ZCA")
#         @test_throws Exception normalizeDNA(42)
#         c = normalizeDNA('C')
#         @test c == dna"c"
#         @test c isa LongSequence
#     end #  normalizeDNA

#     @testset "gc_content" begin
#         @test gc_content(dna"ANTG") == 0.25
#         @test gc_content(dna"cccggg") * 100 == 100.0
#         @test gc_content(dna"ATta") == 0.0
#     end #  composition

#     @testset "composition" begin
#         seq = rand(['A','T','G','C','N'], 20) |> join |> LongDNASeq
#         bc = composition(seq)

#         @test bc[DNA_A] == count(==(DNA_A), collect(seq))
#         @test bc[DNA_C] == count(==(DNA_C), collect(seq))
#         @test bc[DNA_G] == count(==(DNA_G), collect(seq))
#         @test bc[DNA_T] == count(==(DNA_T), collect(seq))
#         @test bc[DNA_N] == count(==(DNA_N), collect(seq))
#     end #  gc_content

#     @testset "complement" begin
#         @test complement(dna"ATTAN") == dna"TAATN"
#         @test complement(dna"gcta") == dna"CGAT"
#         @test complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  complement

#     @testset "reverse_complement" begin
#         @test reverse_complement(dna"ATTAN") == dna"NTAAT"
#         @test reverse_complement(dna"gcta") == dna"TAGC"
#         @test reverse_complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  reverse_complement

#     @testset "parse_fasta" begin
#         testpath = normpath(joinpath(@__DIR__, "..", "data"))
#         genomes = joinpath(testpath, "cov2_genomes.fasta")
#         ex1_path = joinpath(testpath, "ex1.fasta")
#         ex2_path = joinpath(testpath, "ex2.fasta")

#         ex1 = parse_fasta(ex1_path)
#         @test ex1 isa Tuple
#         @test all(x-> x isa String, ex1[1])
#         @test all(x-> x isa LongSequence, ex1[2])

#         @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
#         @test ex1[2] == [dna"AATTATAGC", dna"CGCCCCCCAGTCGGATT"]

#         @test_throws Exception parse_fasta(ex2_path)

#         cov2 = parse_fasta(genomes)
#         @test length(cov2[1]) == 8
#         @test length(cov2[2]) == 8
#     end # parse_fasta

# end # BioSequences

end # BioinformaticsBISC195


text:
# Types

#= ALWB has a lot of data sources - with varying
degrees of nesting. Here wee use types to navigate
and structure the nesting, instead of just using
symbols.

Symbols may be better, idk. They don't need to be imported which is good,
but there is more room for mistakes/confusion around grounpings
and categories with symbols - which this dataset has a lot of.
=#

const ALWB_URL = "http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS"

abstract type DataMode end

"""
Get as the regular measured values.
"""
struct Values <: DataMode end
"""
Get the dataset in relative deciles.
"""
struct Deciles <: DataMode end

const ALWBperiod = Union{Day,Month,Year}

struct ALWB{M<:DataMode,D<:ALWBperiod} <: SDMDataSource end

# Precipitation
struct Precipiation end

# SoilMoisture
struct SoilMoisture{X} end
struct Upper end
struct Lower end
struct Deep end
struct RootZone end

# Runoff
"""
    Runoff

Runoff in mm. Only one dataset.
Available in `Values` and `Deciles` for `Day`, `Month`, and `Year`.
"""
struct Runoff end

# Evapotranspiration

"""
    Evapotrans{X}

[`Actual`](@ref), [`Potential`](@ref), [`Refcrop`](@ref) datasets.
"""
struct Evapotrans{X} end

struct Actual end

struct Potential{X} end
struct Landscape end
struct Areal end
struct SyntheticPan end

struct RefCrop{X} end
struct Short end
struct Tall end

# Evaporation
struct Evaporation{X} end
struct OpenWater end

# Deep Drainage
struct DeepDrainage end

# Interface methods

rasterlayers(::Type{<:ALWB}) =
    (SoilMoisture{Upper}, SoilMoisture{Lower}, SoilMoisture{Deep}, SoilMoisture{RootZone},
     Evapotrans{Actual}, Evapotrans{Potential{Landscape}}, Evapotrans{Potential{Areal}},
     Evapotrans{Potential{SyntheticPan}}, Evapotrans{RefCrop{Short}}, Evapotrans{RefCrop{Tall}},
     Precipiation, Runoff, Evaporation{OpenWater}, DeepDrainage)

rastername(T::Type{<:ALWB{M,P}}, layer, date) where {M,P} =
    _pathsegment(layer) * _pathsegment(P, date) * ".nc"

rasterpath(::Type{ALWB}) = joinpath(rasterpath(), "ALWB")
rasterpath(::Type{ALWB{M,P}}) where {M,P} =
    joinpath(joinpath(rasterpath(), "ALWB"), map(_pathsegment, (M, P))...)
rasterpath(T::Type{<:ALWB}, layer) = joinpath(rasterpath(T), _pathsegment(layer))
rasterpath(T::Type{<:ALWB}, layer, date) =
    joinpath(rasterpath(T), rastername(T, layer, date))

rasterurl(T::Type{<:ALWB{M,P}}, layer, date) where {M,P} =
    joinpath(ALWB_URL, _pathsegment(T), rastername(T, layer, date))

download_raster(T::Type{<:ALWB}, layers::Tuple=rasterlayers(T); kwargs...) =
    map(l -> download_raster(T, l; kwargs...), layers)
function download_raster(T::Type{<:ALWB{M,P}}, layer::Type; dates) where {M,P}
    dates = _date_sequence(dates, P(1))
    mkpath(rasterpath(T, layer))
    raster_paths = String[]
    for d in dates
        url = rasterurl(T, layer, d)
        path = rasterpath(T, layer, d)
        _maybe_download(url, path)
        push!(raster_paths, path)
    end
    raster_paths
end

# Utilitiy methods

_date_sequence(dates::AbstractArray, step) = dates
_date_sequence(dates::Tuple, step) = first(dates):step:last(dates)

_date2string(t, date) = Dates.format(date, _dateformat(t))
_string2date(t, d::AbstractString) = Date(d, _dateformat(t))

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/rain_day_2017.nc
_pathsegment(::Type{Precipiation}) = "rain_day"

_pathsegment(::Type{SoilMoisture{Upper}}) = "s0_pct"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/ss_pct_2017.nc
_pathsegment(::Type{SoilMoisture{Lower}}) = "ss_pct"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/sd_pct_2017.nc
_pathsegment(::SoilMoisture{Deep}) = "sd_pct"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/sm_pct_2017.nc
_pathsegment(::Type{SoilMoisture{RootZone}}) = "sm_pct"

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/qtot_2017.nc
_pathsegment(::Type{Runoff}) = "qtot"

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/etot_2017.nc
_pathsegment(::Type{Evapotrans{Actual}}) = "etot"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/e0_2017.nc
_pathsegment(::Type{Evapotrans{Potential{Landscape}}}) = "e0"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/ma_wet_2017.nc
_pathsegment(::Type{Evapotrans{Potential{Areal}}}) = "ma_wet"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/pen_pet_2017.nc
_pathsegment(::Type{Evapotrans{Potential{SyntheticPan}}}) = "pen_pet"

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/fao_pet_2017.nc
_pathsegment(::Type{Evapotrans{RefCrop{Short}}}) = "fao_pet"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/asce_pet_2017.nc
_pathsegment(::Type{Evapotrans{RefCrop{Tall}}}) = "asce_pet"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/etot_2017.nc

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/msl_wet_2017.nc
_pathsegment(::Type{Evaporation{OpenWater}}) = "msl_wet"

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/dd_2017.nc
_pathsegment(::Type{DeepDrainage}) = "dd"

_pathsegment(::Type{ALWB{M,P}}) where {M,P} = joinpath(_pathsegment(M), _pathsegment(P))
_pathsegment(::Type{Values}) = "values"
_pathsegment(::Type{Deciles}) = "deciles"
_pathsegment(::Type{Day}) = "day"
_pathsegment(::Type{Month}) = "month"
_pathsegment(::Type{Year}) = "year"
# Days are in whole-year files
_pathsegment(::Type{Day}, date) = "_" * string(year(date))
# Months and years are all in one file
_pathsegment(::Type{<:Union{Year,Month}}, date) = ""


text:
module FlowDiagnostics

export RichardsonNumber, RossbyNumber
export ErtelPotentialVorticity, ThermalWindPotentialVorticity
export IsotropicBuoyancyMixingRate, AnisotropicBuoyancyMixingRate
export IsotropicTracerVarianceDissipationRate, AnisotropicTracerVarianceDissipationRate

using Oceananigans
using Oceananigans.Operators
using Oceananigans.AbstractOperations
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.Grids: Center, Face

# Some useful operators
@inline fψ²(i, j, k, grid, f, ψ) = @inbounds f(i, j, k, grid, ψ)^2

function RichardsonNumber(model; N²_bg=0, dUdz_bg=0, dVdz_bg=0)
    u, v, w = model.velocities
    b = model.tracers.b

    dBdz_tot = ∂z(b) + N²_bg
    dUdz_tot = ∂z(u) + dUdz_bg
    dVdz_tot = ∂z(v) + dVdz_bg

    return dBdz_tot / (dUdz_tot^2 + dVdz_tot^2)
end

function RossbyNumber(model; dUdy_bg=0, dVdx_bg=0, f=nothing)
    u, v, w = model.velocities
    if f==nothing
        f = model.coriolis.f
    end

    dUdy_tot = ∂y(u) + dUdy_bg
    dVdx_tot = ∂x(v) + dVdx_bg

    return (dVdx_tot - dUdy_tot) / f
end

#++++ Potential vorticity
function potential_vorticity_in_thermal_wind_fff(i, j, k, grid, u, v, b, f)

    dVdx =  ℑzᵃᵃᶠ(i, j, k, grid, ∂xᶠᶠᵃ, v) # F, F, C → F, F, F
    dUdy =  ℑzᵃᵃᶠ(i, j, k, grid, ∂yᶠᶠᵃ, u) # F, F, C → F, F, F
    dbdz = ℑxyᶠᶠᵃ(i, j, k, grid, ∂zᵃᵃᶠ, b) # C, C, F → F, F, F

    pv_barot = (f + dVdx - dUdy) * dbdz

    dUdz = ℑyᵃᶠᵃ(i, j, k, grid, ∂zᵃᵃᶠ, u) # F, C, F → F, F, F
    dVdz = ℑxᶠᵃᵃ(i, j, k, grid, ∂zᵃᵃᶠ, v) # C, F, F → F, F, F

    pv_baroc = -f*(dUdz^2 + dVdz^2)

    return pv_barot + pv_baroc
end

function ThermalWindPotentialVorticity(model; f=nothing)
    u, v, w = model.velocities
    b = model.tracers.b
    if f==nothing
        f = model.coriolis.f
    end
    return KernelFunctionOperation{Face, Face, Face}(potential_vorticity_in_thermal_wind_fff, model.grid;
                                                     computed_dependencies=(u, v, w, b), parameters=f)
end

function ertel_potential_vorticity_fff(i, j, k, grid, u, v, w, b, params)

    dWdy =  ℑxᶠᵃᵃ(i, j, k, grid, ∂yᶜᶠᵃ, w) # C, C, F  → C, F, F  → F, F, F
    dVdz =  ℑxᶠᵃᵃ(i, j, k, grid, ∂zᵃᵃᶠ, v) # C, F, C  → C, F, F  → F, F, F
    dbdx = ℑyzᵃᶠᶠ(i, j, k, grid, ∂xᶠᶜᵃ, b) # C, C, C  → F, C, C  → F, F, F
    pv_x = (params.fx + dWdy - dVdz) * dbdx # F, F, F

    dUdz =  ℑyᵃᶠᵃ(i, j, k, grid, ∂zᵃᵃᶠ, u) # F, C, C  → F, C, F → F, F, F
    dWdx =  ℑyᵃᶠᵃ(i, j, k, grid, ∂xᶠᶜᵃ, w) # C, C, F  → F, C, F → F, F, F
    dbdy = ℑxzᶠᵃᶠ(i, j, k, grid, ∂yᶜᶠᵃ, b) # C, C, C  → C, F, C → F, F, F
    pv_y = (params.fy + dUdz - dWdx) * dbdy # F, F, F

    dVdx =  ℑzᵃᵃᶠ(i, j, k, grid, ∂xᶠᶠᵃ, v) # C, F, C  → F, F, C → F, F, F
    dUdy =  ℑzᵃᵃᶠ(i, j, k, grid, ∂yᶠᶠᵃ, u) # F, C, C  → F, F, C → F, F, F
    dbdz = ℑxyᶠᶠᵃ(i, j, k, grid, ∂zᵃᵃᶠ, b) # C, C, C  → C, C, F → F, F, F
    pv_z = (params.fz + dVdx - dUdy) * dbdz

    return pv_x + pv_y + pv_z
end

function ErtelPotentialVorticity(model; location = (Face, Face, Face))
    if location != (Face, Face, Face)
        throw(ArgumentError("ErtelPotentialVorticity only implemented at location (Face, Face, Face) for now."))
    end

    u, v, w = model.velocities
    if ~(model.background_fields.velocities.u isa Oceananigans.Fields.ZeroField)
        u += model.background_fields.velocities.u
    end
    if ~(model.background_fields.velocities.v isa Oceananigans.Fields.ZeroField)
        v += model.background_fields.velocities.v
    end
    if ~(model.background_fields.velocities.w isa Oceananigans.Fields.ZeroField)
        w += model.background_fields.velocities.w
    end

    b = model.tracers.b
    if ~(model.background_fields.tracers.b isa Oceananigans.Fields.ZeroField)
        b += model.background_fields.tracers.b
    end

    coriolis = model.coriolis
    if coriolis isa FPlane
        fx = fy = fz = model.coriolis.f
    elseif coriolis isa ConstantCartesianCoriolis
        fx = coriolis.fx
        fy = coriolis.fy
        fz = coriolis.fz
    else
        throw(ArgumentError("ErtelPotentialVorticity only implemented for FPlane and ConstantCartesianCoriolis"))
    end

    return KernelFunctionOperation{Face, Face, Face}(ertel_potential_vorticity_fff, model.grid;
                                                     computed_dependencies=(u, v, w, b), parameters=(; fx, fy, fz))
end
#----

#+++++ Tracer variance dissipation
function isotropic_tracer_variance_dissipation_rate_ccc(i, j, k, grid, b, κᵇ)
    dbdx² = ℑxᶜᵃᵃ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, b) # C, C, C  → F, C, C  → C, C, C
    dbdy² = ℑyᵃᶜᵃ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, b) # C, C, C  → C, F, C  → C, C, C
    dbdz² = ℑzᵃᵃᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, b) # C, C, C  → C, C, F  → C, C, C

    return 2 * κᵇ[i,j,k] * (dbdx² + dbdy² + dbdz²)
end
function IsotropicTracerVarianceDissipationRate(model, b, κᵇ; location = (Center, Center, Center))
    if location == (Center, Center, Center)
        return KernelFunctionOperation{Center, Center, Center}(isotropic_tracer_variance_dissipation_rate_ccc, model.grid;
                                       computed_dependencies=(b, κᵇ))
    else
        throw(Exception)
    end
end

function anisotropic_tracer_variance_dissipation_rate_ccc(i, j, k, grid, b, params)
    dbdx² = ℑxᶜᵃᵃ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, b) # C, C, C  → F, C, C  → C, C, C
    dbdy² = ℑyᵃᶜᵃ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, b) # C, C, C  → C, F, C  → C, C, C
    dbdz² = ℑzᵃᵃᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, b) # C, C, C  → C, C, F  → C, C, C

    return 2 * (params.κx*dbdx² + params.κy*dbdy² + params.κz*dbdz²)
end
function AnisotropicTracerVarianceDissipationRate(model, b, κx, κy, κz; location = (Center, Center, Center))
    if location == (Center, Center, Center)
        return KernelFunctionOperation{Center, Center, Center}(anisotropic_tracer_variance_dissipation_rate_ccc, model.grid;
                                       computed_dependencies=(b,),
                                       parameters=(κx=κx, κy=κy, κz=κz),)
    else
        throw(Exception)
    end
end
#-----

end # module


text:
module LidJul

    export BoundaryCondition,dirichlet,neumann

    @enum BoundaryCondition dirichlet=0 neumann=1

    include("laplacian1D.jl")
    include("laplacian2D.jl")
    # include("poisson2D_ADI.jl")
    # include("poisson2D_CG.jl")
    # include("poisson2D_PCG.jl")
    include("poisson2D_TT.jl")
    # include("poisson2DSparse.jl")
    include("poisson2DSparseLU.jl")
    include("poisson2DSparseAMG.jl")
    include("poisson2DSparseCGILU.jl")
    include("poisson2DSparseIterative.jl")
    # include("ADISmoother.jl")
    include("GSSmoother.jl")
    # include("GSRBSmoother.jl")
    # include("CGSmoother.jl")
    # include("BlockJacobiSmoother.jl")

    include("poisson2DGMG.jl")
    include("poisson2DGMG_new.jl")
    # include("RedBlackArray.jl")
    # include("poisson2DGMGRB.jl")
    include("TensorialOperator.jl")

    # include("poisson2DSparse.jl")
end # module