text:
# MOON rewrite AsgnArray for matrix, verilog bus notation
mutable struct Arrow
    WireType::Int
        # 0 -> wire
        # 1 -> input
        # 2 -> output
    Reg::Bool
    Endian::Int
        # 0 -> Undeclared
        # 1 -> Big Endian
        # 2 -> Little Endian
    BitLeft::Int
    BitRight::Int
    BitCount::Int
    Assigned::Int
        # 0 -> Unassigned
        # 1 -> Assigned
        # 2 -> Bitwise Assigned
    AsgnArray::Array{Int}
end

# MOON Add the ex so that the print statements can include it
function pushArrowJSequals(name::Symbol, wt::Int, endi::Int, bl::Int, br::Int, bc::Int)
    global arrowDict
    if haskey(arrowDict, name)
        error("Attempting to create wire $(name) after creation")
    end

    # Creating a whole new wire
    arrowDict[name] = Arrow(wt, false, endi, bl, br, bc, 0, fill!(Array{Bool}(bc), 1))
end

function pushArrowJSassignment(name::Symbol, endi::Int, bl::Int, br::Int, bc::Int)
    global arrowDict
    if haskey(arrowDict, name)
        a = arrowDict[name]
        if a.Assigned != 0
            error("Attempting to assign a whole wire after assignment: $(name)")
        end
        if a.WireType == 1
            error("Attempting to assign an Input: $(name)")
        end

        # Must be assigning a previously created wire
        arrowDict[name].Assigned = 1
        arrowDict[name].AsgnArray = fill!(Array{Int}(bc), 1)
    else
        # Creating a whole new wire
        arrowDict[name] = Arrow(0, false, endi, bl, br, bc, 1, fill!(Array{Int}(bc), 1))
    end
end

function pushArrowBitsJSassignment(name::Symbol, endi::Int, bl::Int, br::Int, bc::Int)
    global arrowDict
    if haskey(arrowDict, name)
        a = arrowDict[name]
        if     a.Assigned == 1
            error("Attempting to assign individual bits after a whole wire assignment: $(name)")
        elseif a.Assigned == 0
            # this is likely an erroneous error statement, and likely needs removing
            if     a.Endian == 0 && bc > 1
                error("Attempting to assign multiple bits of a wire which lacks Endianness: $(name)")
            end
            arrowDict[name].Assigned = 2
        end

        minb = min(bl, br)
        maxb = max(bl, br)
        arrowminb = min(a.BitRight, a.BitLeft)
        arrowmaxb = max(a.BitRight, a.BitLeft)
        if     minb < arrowminb
            error("Attempting to assign bit $(minb) which is outside $(name)")
        elseif maxb > arrowmaxb
            error("Attempting to assign bit $(maxb) which is outside $(name)")
        end

        if a.Assigned == 2
            overwrites = find(arrowDict[name].AsgnArray[(minb-arrowminb+1):(maxb-arrowminb+1)] != 0)
            if overwrites != []
                error("Assigning to bits $(overwrites + minb + arrowminb) after they were already assigned: $(name)")
            end
        end

        tmp = maximum(arrowDict[name].AsgnArray) + 1
        arrowDict[name].AsgnArray[(minb-arrowminb+1):(maxb-arrowminb+1)] = tmp
    else
        # Creating a whole new wire
        arrowDict[name] = Arrow(0, false, endi, bl, br, bc, 2, fill!(Array{Int}(bc), 1))
    end
end

function updateArrowReg(name::Symbol)
    global arrowDict
    arrowDict[name].Reg = true
end

function JSarrowDictCleaning()
    global arrowDict

    for (k, v) in arrowDict
        if v.WireType != 1
            # Check to see if any arrows in arrowDict are created but not assigned
            if v.Assigned == 0
                warn("Arrow $(k) was created but never assigned")
            end

            # Check to see if any bits in an arrow are created but not assigned
            bmin = min(v.BitLeft, v.BitRight)
            bmax = max(v.BitLeft, v.BitRight)
            for i = 1:length(v.AsgnArray)
                bit = v.AsgnArray[i]
                if bit == 0
                    warn("Bit $(bit + bmin) in Arrow $(k) was created but never assigned")
                    arrowDict[k].Assigned = 2
                end
            end

            # Change assigned to 1 from 2 if all bits are assigned at once
            if v.Assigned == 2
                multipleAsgns = find(v.AsgnArray != 1)
                if multipleAsgns == []
                    arrowDict[k].Assigned = 1
                end
            end
        end

        # Check to see if any arrows have an Undeclared Endianness
        if v.Endian == 0
            warn("Endianness of Arrow $(k) is never properly defined")
        end
    end
end

# Create Arrows for each wire
# Initialize wires if they lack it
# Make certain the user is not attempting to
#   write or read outside of a wires access
#   Read a wire in the incorrect endian
#   write a wire which has already been written to
#   write to an input
# Throws a warning if a wire is created but never assigned (except for inputs)
# Removes lines where BitRepeat is zero
function sanitize(ex::Expr)
    # Removes lines and exprs which are a BitRepeat of zero
    # MOON

    # Create a global arrowDict
    global arrowDict = Dict{Symbol, Arrow}()

    #  Sanatize expr block and fill the new arrowDict
    jexpr = Expr(:function)
    jargs = Array{Any}(2)
    jargs[1] = Expr(:call)
    jargs[1].args = ex.args[1]
    jargs[2] = JSblock(ex.args[2])
    jexpr.args = jargs

    JSarrowDictCleaning()

    return jexpr, arrowDict
end

function JSblock(ex::Expr)
    l = length(ex.args)
    i = 1
    while i <= l
        arg = ex.args[i]
        if     isa(arg, Expr)
            jexpr = JSblockhelper(arg)
            if     isa(jexpr, Expr)
                if jexpr.head == :block
                    ex.args[i] = jexpr.args[1]
                    m = length(jexpr.args)
                    for j = 2:m
                        insert!(ex.args, i + 1, jexpr.args[j])
                        i = i + 1
                        l = l + 1
                    end
                else
                    ex.args[i] = jexpr
                end
            else
                error("Something Weird:\n$(ex)\n$(jexpr)")
            end
        else
            error("Found a non-expr in Juliog after parameterization:\n$(arg)")
        end
        i = i + 1
    end
    return ex
end

# MOON check for matrix notation
function JSblockhelper(ex::Expr)
    h = ex.head
    if     h == :(=)
        return JSequals(ex)
    elseif h == :(:=)
        return JSassignment(ex)
    elseif h == :if
        return JSif(ex)
    elseif h == :macrocall
        error("Macro statements are not completed in JuliogSanitize")
        return JSmacro(ex)
    else
        error("hit unexpected call $(j) in JSblockhelper:\n$(ex)")
    end
end

function JSequals(ex::Expr)
    # Checking to see if this expression creates a new wire, input, or output
    if isa(ex.args[2], Expr)
        if isa(ex.args[2].args[1], Expr)
            sym = ex.args[2].args[1].args[1]
            wireType = -1
            if     sym == :Wire
                wireType = 0
            elseif sym == :Input
                wireType = 1
            elseif sym == :Output
                wireType = 2
            end

            if wireType != -1
                if isa(ex.args[2].args[2], Expr)
                    # Multiple bit creation
                    bl = ex.args[2].args[2].args[1]
                    br = ex.args[2].args[2].args[2]
                    if bl >= br
                        pushArrowJSequals(ex.args[1], wireType, 1, bl, br, bl-br+1)
                    else
                        pushArrowJSequals(ex.args[1], wireType, 2, bl, br, br-bl+1)
                    end
                else
                    # must be a single bit creation
                    b = ex.args[2].args[2]
                    pushArrowJSequals(ex.args[1], wireType, 0, b, b, 1)
                end
            end
            return ex
        end
    end

    return JSassignment(ex)
end

function JSassignment(ex::Expr)
    global arrowDict

    # Check righthand side for read errors and bit counts
    rhs = ex.args[2]

    bcrhs, endirhs = JSrhs(rhs)

    # Use endirhs and bcrhs to assign the new wire    jexpr = Expr(:block)
    if isa(ex.args[1], Expr)
        # lhs is a reference
        return JSassignmentref(ex, endirhs, bcrhs)
    else
        # lhs is just a symbol for a wire
        return JSassignmentsymbol(ex, endirhs, bcrhs)
    end
end

function JSrhs(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

function JSassignmentref(ex::Expr, endirhs::Int, bcrhs::Int)
    global arrowDict

    # LHS is a reference
    #   Maybe a vcat or tuple, can't do anything about that yet
    name = ex.args[1].args[1]

    # Check to see if the rhs and lhs bitcounts agree
    #  Throw a warning if that is not the case
    arg = ex.args[1].args[2]
    if isa(arg, Expr)
        # Multiple bit or single bit
        bl = arg.args[1]
        l = length(arg.args)
        if l == 3
            br = arg.args[3]
        else
            br = arg.args[2]
        end

        # Determine the endianness and bitcount of the lhs
        if bl == br
            bclhs = 1
            endilhs = 0
        elseif bl > br
            bclhs = bl - br + 1
            endilhs = 1
        else
            bclhs = br - bl + 1
            endilhs = 2
        end

        if     bclhs > bcrhs
            warn("Zero Bit Padding necessitated: $(bclhs) bits on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        elseif bclhs < bcrhs
            warn("Truncation of Bits necessitated: $(bclhs) bits on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        end

        if     endilhs == 1 && endirhs == 2
            error("Endian Mismatch: LHS is Big Endian, RHS is Little Endian:\n$(ex)")
        elseif endilhs == 2 && endilhs == 1
            error("Endian Mismatch: LHS is Little Endian, RHS is Big Endian:\n$(ex)")
        end

        retExpr = Expr(:(=), name, Expr(:(ref), Expr(:call, :Wire), Expr(:(:), bl,br)))

    else
        # Single bit
        bclhs = 1
        bl = arg
        br = arg
        endilhs = 0
        if     1 > bcrhs
            warn("Zero Bit Padding necessitated: $(1) bit on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        elseif 1 < bcrhs
            warn("Truncation of Bits necessitated: $(1) bit on lh, and $(bcrhs) bits on rhs:\n$(ex)")
        end

        retExpr = Expr(:(=), name, Expr(:(ref), Expr(:call, :Wire), 0))
    end

    if haskey(arrowDict, name)
        jexpr = ex
        jexpr.head = :(:=)
    else
        jexpr = Expr(:block)
        jargs = Array{Any}(2)
        jargs[1] = retExpr
        jargs[2] = ex
        jargs[2].head = :(:=)
        jexpr.args = jargs
    end

    pushArrowBitsJSassignment(name, endilhs, bl,br,bclhs)
    return jexpr
end

function JSassignmentsymbol(ex::Expr, endirhs::Int, bcrhs::Int)
    global arrowDict

    # Must be a wire name
    name = ex.args[1]

    # If the lhs and rhs bits don't agree, throw out warnings
    if haskey(arrowDict, name)
        bclhs = arrowDict[name].BitCount
        if     bclhs > bcrhs
            warn("Zero Bit Padding necessitated: $(bclhs) bits on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        elseif bclhs < bcrhs
            warn("Truncation of Bits necessitated: $(bclhs) bit on lh, and $(bcrhs) bits on rhs:\n$(ex)")
        end
    end

    # Add the new linewire to the arrowDict
    # and return the new wire instantiation expr so it can be added to the block

    #  Must be a wireType of wire
    #  Is not a reg since it is outside of an if or macro
    #  The endian is inherited from rhs with associated bitleft and bitright
    #  LineWire by definition
    #  It is being assigned this line
    if bcrhs == 1
        bl = 0
        br = 0
        retExpr = Expr(:(=), name, Expr(:(ref), Expr(:call, :Wire), 0))
    else
        if endirhs == 2
            bl = 0
            br = bcrhs - 1
        else
            # if enddianness if undeclared, then big notation will be adopted
            #  in such a case big endian operation may not be accurate
            bl = bcrhs - 1
            br = 0
        end
        retExpr = Expr(:(=), name, Expr(:(ref), Expr(:call, :Wire), Expr(:(:), bl,br)))
    end

    if haskey(arrowDict, name)
        jexpr = ex
        jexpr.head = :(:=)
    else
        jexpr = Expr(:block)
        jargs = Array{Any}(2)
        jargs[1] = retExpr
        jargs[2] = ex
        jargs[2].head = :(:=)
        jexpr.args = jargs
    end

    pushArrowJSassignment(name, endirhs, bl,br,bcrhs)
    return jexpr
end

# The righthand side of an assignment
# Makes sure symbols are instantiated wires
#  and that all bit reads satisfy wire bitcounts, bitends, and bitstarts
function JSrhsexpr(ex::Expr)
    h = ex.head
    if     h == :ref
        return JSrhsref(ex)
    elseif h == :vcat
        return JSrhsvcat(ex)
    elseif h == :call
        return JSrhscall(ex)
    elseif h == :if
        return JSrhsif(ex)
    else
        error("Hit unexpected expression symbol $(h) in JSreadrhs on expression:\n$(ex)")
    end
end

# MOON add checking for (and removing?) verilog notation
function JSrhsref(ex)
    global arrowDict

    # Check to see if the rhs is calling a symbol which is not in arrowDict
    name = ex.args[1]
    if !haskey(arrowDict, name)
        error("Discovered uninstantiated symbol: $(name) in expr:\n$(ex)")
    end
    arrow = arrowDict[name]

    # Throw a warning if the rhs is calling a wire which has not been assigned
    #  but only if that wire is not also an input
    if arrow.Assigned == 0 && arrow.WireType != 1
        warn("Wire $(name) is called before being assigned:\n$(ex)")
    end
    arrowbl = arrow.BitLeft
    arrowbr = arrow.BitRight
    arrowendi = arrow.Endian

    # Check to see if the reference is single-bit or multiple-bit

    arg = ex.args[2]
    if isa(arg, Expr)
        # Must be multiple-bit reference or single-bit double-qouted

        l = length(arg.args)
        bl = arg.args[1]

        # Check to see if reverse indexing
        # And then check to see indexing in the correct Endian
        if l == 3
            # Reverse indexing
            br = arg.args[3]
            if arg.args[2] != -1
                error("Second Number in $(ex.args[2]) is not -1, Failed attempt to reverse indexing in:\n$(ex)")
            end
            br = arg.args[3]
            if (arrowendi == 1 && bl > br) || (arrowendi == 2 && br > bl)
                error("Attempting to reference $(name) in incorrect Endian:\n$(ex)")
            end

        elseif l == 2
            br = arg.args[2]
            br = ex.args[2].args[2]
            if (arrowendi == 1 && br > bl) || (arrowendi == 2 && bl > br)
                error("Attempting to reference $(name) in incorrect Endian:\n$(ex)")
            end
        else
            error("In JSrhsref, have a ref with irregular length:\n$(ex)")
        end

        # Check to see if referening bits outside wire limits
        if arrowendi == 2
            if br > arrowbr
                error("Attempting to reference bit $(bl) which is outside $(name):\n$(ex)")
            end
            if bl < arrowbl
                error("Attempting to reference bit $(br) which is outside $(name):\n$(ex)")
            end
        else
            if br < arrowbr
                error("Attempting to reference bit $(br) which is outside $(name):\n$(ex)")
            end
            if bl > arrowbl
                error("Attempting to reference bit $(bl) which is outside $(name):\n$(ex)")
            end
        end

        # If nothing has failed, then get final bc and endi
        if     bl > br
            bc =  bl - br + 1
            endi = 1
        elseif br > bl
            bc = br - bl + 1
            endi = 2
        else
            bc = 1
            endi = 0
        end
    else
        # Must be single-bit reference
        b = ex.args[2]

        if arrowendi == 2
            if b > arrowbr
                error("Attempting to reference bit $(b) which is outside $(name):\n$(ex)")
            end
            if b < arrowbl
                error("Attempting to reference bit $(b) which is outside $(name):\n$(ex)")
            end
        else
            if b < arrowbr
                error("Attempting to reference bit $(b) which is outside $(name):\n$(ex)")
            end
            if b > arrowbl
                error("Attempting to reference bit $(b) which is outside $(name):\n$(ex)")
            end
        end

        # if nothing has failed, take default bc and endi
        bc = 1
        endi = 0
    end

    return bc, endi
end

# Helps to get rid of Verilog bus notation if present
function JSrhsrefhelper(ex::Expr)
    if isa(ex.args[1], Expr) && ex.args[1].head == :ref
        jargs = JPrefhelper(ex.args[1])
        append!(jargs,ex.args[2:end])
        return jargs
    end

    return ex.args
end

function JSrhsvcat(ex)
    global arrowDict
    l = length(ex.args)
    bctotal = 0
    enditotal = 0
        # 0 is undefined
        # 1 is big
        # 2 is little
    for i = 1:l
        # Establish Endianness and bitcount for this particular index
        # Throw warnings if endianness conflicts with previous endianness
        arg = ex.args[i]
        if     isa(arg, Expr)
            # is a ref or bit repeat
            if arg.head == :call
                # is a bit repeat
                bc, endi = JSrhscall(arg)
            else
                # is a ref
                bc, endi = JSrhsref(arg)
            end
            if enditotal != 0 && endi != enditotal
                warn("Attempting to concatenate BigEndian and LittleEndian wires:\n$(ex)")
            end
            enditotal = endi
            bctotal = bctotal + bc
        elseif isa(arg, Symbol)
            if !haskey(arrowDict, arg)
                error("Discovered uninstantiated symbol: $(arg) in expr:\n$(ex)")
            end
            endi = arrowDict[arg].Endian
            bc = arrowDict[arg].BitCount
            if enditotal != 0 && endi != enditotal
                warn("Attempting to concatenate BigEndian and LittleEndian wires:\n$(ex)")
            end
            enditotal = endi
            bctotal = bctotal + bc
        elseif arg==0 || arg==1 || arg == 'x' || arg == 'z'
            endi = 0
            bctotal = bctotal + 1
        else
            error("Attempting to vcat with an incorrect wire entry:$(arg) in:\n$(ex)")
        end
    end

    return bctotal, 0
end

function JSrhscall(ex)

    # Is a bit repeat
    if isa(ex.args[1], Expr)
        bc, endi = JSrhsref(ex.args[1])
        rep = ex.args[2]
        if rep < 0
            error("Attempting to BitRepeat $(ex) with a negative value")
        end
        bc = bc * rep
        return bc, endi
    end

    # is a logic
    opd1 = ex.args[2]
    if     isa(opd1, Expr)
        bc1, endi1 = JSrhsexpr(opd1)
    elseif isa(opd1, Symbol)
        global arrowDict
        endi1 = arrowDict[opd1].Endian
        bc1 = arrowDict[opd1].BitCount
    elseif isa(opd1, Int)
        bc1 = 1
        endi1 = 0
    elseif isa(opd1, Char)
        bc1 = 1
        endi1 = 0
    end

    # Check to see if it is a 2-op logic operation
    if length(ex.args) == 3
        opd2 = ex.args[3]
        if     isa(opd2, Expr)
            bc2, endi2 = JSrhsexpr(opd2)
        elseif isa(opd2, Symbol)
            global arrowDict
            endi2 = arrowDict[opd2].Endian
            bc2 = arrowDict[opd2].BitCount
        elseif isa(opd2, Int)
            bc2 = 1
            endi2 = 0
        elseif isa(opd2, Char)
            bc2 = 1
            endi2 = 0
        end

        if     endi1 == 0 && endi2 == 0
            endipredicted = 0
        elseif endi1 == 0
            endipredicted = endi2
        elseif endi2 == 0
            endipredicted = endi1
        elseif endi1 != endi2
            endipredicted = 0
        else
            endipredicted = endi1
        end
    end

    o = Symbol(ex.args[1])
    if     o == :&
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        bctotal = max(bc1, bc2)
        return bctotal, endipredicted
    elseif o == :|
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        bctotal = max(bc1, bc2)
        return bctotal, endipredicted
    elseif o == :~
        return bc1, endi1
    elseif o == :+
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        bctotal = max(bc1, bc2) + 1
        return bctotal, endipredicted
    elseif o == :-
        bctotal = max(bc1, bc2) + 1
        return bctotal, endipredicted
    elseif o == :*
        JSrctwoloosest(o, ex, endi1, endi2)
        bctotal = bc1 + bc2
        return bctotal, endipredicted
    elseif o == :/
        JSrctwoloosest(o, ex, endi1, endi2)
        bctotal = bc1 + bc2
        return bctotal, endipredicted
    elseif o == :%
        JSrctwoloosest(o, ex, endi1, endi2)
        bctotal = bc1 + bc2
        return bctotal, endipredicted
    elseif o == :(==)
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        return 1, endipredicted
    elseif o == :(!=)
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        return 1, endipredicted
    elseif o == :>
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        return 1, endipredicted
    elseif o == :<
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        return 1, endipredicted
    elseif o == :(>=)
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        return 1, endipredicted
    elseif o == :(<=)
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        return 1, endipredicted
    elseif o == :(<<)
        JSrcshift(o, ex, bc1, opd2)
        return bc1, endi1
    elseif o == :(>>)
        JSrcshift(o, ex, bc1, opd2)
        return bc1, endi1
    elseif o == :(>>>)
        JSrcshift(o, ex, bc1, opd2)
        return bc1, endi1
    elseif o == :^
        JSrctwoloose(o, ex, bc1, bc2, endi1, endi2)
        bctotal = max(bc1, bc2)
        return bctotal, endipredicted
    else
        error("Hit unexpected head $(o) in $(ex)")
    end
end

function JSrctwostrict(op::Symbol, ex::Expr, bc1::Int, bc2::Int, endi1::Int, endi2::Int)
    if bc1 != bc2
        error("Attempting to $(op) on wires of different bit count:\n$(ex)")
    end

    if endi1 != endi2
        warn("Attempting to $(op) on wires of different Endiannes:\n$(ex)")
    end
end

function JSrctwoloose(op::Symbol, ex::Expr, bc1::Int, bc2::Int, endi1::Int, endi2::Int)
    if bc1 != bc2
        warn("Attempting to $(op) on wires of different bit count:\n$(ex)")
    end

    if endi1 != endi2
        warn("Attempting to $(op) on wires of different Endiannes:\n$(ex)")
    end
end

function JSrctwoloosest(op::Symbol, ex::Expr, endi1::Int, endi2::Int)
    if endi1 != endi2
        warn("Attempting to $(ex.args[1]) on wires of different Endiannes:\n$(ex)")
    end
end

function JSrcshift(op::Symbol, ex::Expr, bc1::Int, shift::Int)
    if shift >= bc1
        warn("Attempting to shift a wire more than it has bits:\n$(ex)")
    end
end

# Even though this looks overly simplistic, this is likely the correct solution
# Must be a tertiary command on the rhs
function JSrhsif(ex::Expr)
    l = length(ex.args)

    condition = ex.args[1]
    if isa(condition, Expr)
        bc, endi = JSrhsexpr(condition)
        if bc != 1
            error("Condition of if statement does not equal 1 bit:\n$(ex)")
        end
    else
        # is a symbol
        global arrowDict
        if isa(condition, Symbol)
            if !haskey(arrowDict, condition)
                error("Discovered uninstantiated symbol: $(arg) in expr:\n$(ex)")
            end
        end
        if arrowDict[condition].BitCount != 1
            error("Condition of if statement does not equal 1 bit:\n$(ex)")
        end
    end

    arg = ex.args[2]
    if isa(arg, Expr)
        bc1, endi1 = JSrhsexpr(arg)
    else
        # is a 0,1,'x','z' or a symbol
        if isa(arg, Symbol)
            global arrowDict
            if !haskey(arrowDict, arg)
                error("Discovered uninstantiated symbol: $(arg) in expr:\n$(ex)")
            end
            a = arrowDict[arg]
            bc1 = a.BitCount
            endi1 = a.Endian
        elseif arg != 0 && arg != 1 && arg != 'x' && arg != 'z'
            error("Attempting to write a bit to an unauthorized value: $(arg)\n$(ex)")
        else
            bc1 = 1
            endi = 0
        end
    end

    arg = ex.args[3]
    if isa(arg, Expr)
        bc2, endi2 = JSrhsexpr(arg)
    else
        # is a 0,1,'x','z' or a symbol
        if isa(arg, Symbol)
            global arrowDict
            if !haskey(arrowDict, arg)
                error("Discovered uninstantiated symbol: $(arg) in expr:\n$(ex)")
            end
            a = arrowDict[arg]
            bc2 = a.BitCount
            endi2 = a.Endian
        elseif arg != 0 && arg != 1 && arg != 'x' && arg != 'z'
            error("Attempting to write a bit to an unauthorized value: $(arg)\n$(ex)")
        else
            bc2 = 1
            endi2 = 0
        end
    end

    bctotal = bc1
    enditotal = endi1
    if bc1 != bc2
        warn("Bit counts differ between if and else statement in:\n$(ex)")
        bctotal = max(bc1, bc2)
    end
    if endi1 != endi2
        warn("Endianness differ between if and else statement in:\n$(ex)")
        enditotal = 0
    end

    return bctotal, enditotal
end

max_stars_repo_name: trevortheblack/Juliog

max_issues_repo_name: trevortheblack/Juliog

max_forks_repo_name: trevortheblack/Juliog


text:
module RoundFast

import Base:(+),(-),(*),(/),sqrt

include("adjacent.jl")
include("eftArith.jl")
include("eftRound.jl")
include("roundedOps.jl")

end # module

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.DAE (ModiaMath/DAE/_module.jl)
#

# Provide functions to handle time and state events in user models.
@eval using Printf

mutable struct EventHandler
   # Input values for the event functions
    time::Float64              # Current simulation time
    initial::Bool              # = true, if model is called at initialization
                              #         (if initial, event=true)
    terminal::Bool             # = true, if model is called for termination (close files, streams, visualization, ...)
    event::Bool                # = true, if model is called at an event
    afterSimulationStart::Bool # = true, if model is called after simulation start
    crossing::Bool             # = true, if model is called to compute crossing function

    # Computed by the event functions
    # For time events:
    maxTime::Float64        # Integrate at most up to maxTime
                            # (if integrateToEvent == true, maxTime = nextEventTime,
                            #  otherwise maxTime > nextEventTime)
    nextEventTime::Float64  # Next time event instant (= realmax(Time) if no event)
    integrateToEvent::Bool  # = true, if integrator shall integrate to nextEventTime
                            # = false, if integrator can integrate beyond nextEventTime
    restart::EventRestart   # Type of integrator restart at current event
    newEventIteration::Bool   # = true, if another event iteration; = false, if no event iteration anymore
    firstEventIteration::Bool # = true, if first iteration at an event.

   # For state events:
    nz::Int                 # Number of event indicators
    z::Vector{Float64}      # Vector of event indicators (zero crossings). If one of z[i] passes
                            # zero, that is beforeEvent(z[i])*z[i] < 0, an event is triggered
                            # (allocated during instanciation according to nz).
    zPositive::Vector{Bool} # = true if z > 0 at the last event instant, otherwise false
    zDir::Vector{Cint}      # zDir[i] =  0: Root is reported for both crossing directions
                            #         =  1: Root is reported when crossing from negative to positive direction
                            #         = -1: Root is reported when crossing from positive to negative direction

    # Logger
    logger::ModiaMath.Logger

    function EventHandler(;nz::Int=0)
        @assert(nz >= 0)
        new(0.0, false, false, false, false, false, floatmax(Float64), floatmax(Float64),
            true, NoRestart, false, false, nz, ones(nz), fill(false, nz), fill(0, nz))
    end
end

positiveCrossingAsString(positive::Bool) = positive ? " (became > 0)" : " (became <= 0)"
negativeCrossingAsString(negative::Bool) = negative ? " (became < 0)" : " (became >= 0)"

function initEventIteration!(h::EventHandler, t::Float64)
    h.time          = t
    h.restart       = NoRestart
    h.maxTime       = floatmax(Float64)
    h.nextEventTime = floatmax(Float64)
    h.newEventIteration   = false
    h.firstEventIteration = true
end

function terminateEventIteration!(h::EventHandler)
    result                = h.initial ? false : !h.newEventIteration
    h.firstEventIteration = false
    h.newEventIteration   = false
    return result
end

isInitial(h::EventHandler)              = h.initial
isTerminal(h::EventHandler)             = h.terminal
isEvent(h::EventHandler)                = h.event
isAfterSimulationStart(h::EventHandler) = h.afterSimulationStart
isEventIndicator(h::EventHandler) = h.crossing
const zEps = 1.e-14

function setNextEvent!(h::EventHandler, nextEventTime::Float64;
                       integrateToEvent::Bool=true,
                       restart::EventRestart=Restart)
    if (h.event && nextEventTime > h.time) || (h.initial && nextEventTime >= h.time)
        if integrateToEvent
            h.maxTime = min(h.maxTime, nextEventTime)
        end
        if nextEventTime < h.nextEventTime
            h.nextEventTime    = nextEventTime
            h.integrateToEvent = integrateToEvent
            if ModiaMath.isLogEvents(h.logger)
                @printf("        nextEventTime = %.6g s, integrateToEvent = %s\n",
                   nextEventTime, integrateToEvent ? "true" : "false")
            end
        end
        h.restart = max(h.restart, restart)
    end
    return nothing
end

function positive!(h::EventHandler, nr::Int, crossing::Float64, crossingAsString::String;
                   restart::EventRestart=Restart)::Bool
    if h.initial
        h.zPositive[nr] = crossing > 0.0
        if ModiaMath.isLogEvents(h.logger)
            println("        ", crossingAsString, " = ", crossing, positiveCrossingAsString(h.zPositive[nr]))
        end

    elseif h.event
        # println("... nr = ", nr, ", crossing = ", crossing)
        new_zPositive = crossing > 0.0
        change = (h.zPositive[nr] && !new_zPositive) || (!h.zPositive[nr] && new_zPositive)
        h.zPositive[nr] = new_zPositive

        if change
            h.restart = max(h.restart, restart)
            if ModiaMath.isLogEvents(h.logger)
                println("        ", crossingAsString, " = ", crossing, positiveCrossingAsString(h.zPositive[nr]))
            end
            h.newEventIteration = true
        end
    end
    h.z[nr] = crossing + (h.zPositive[nr] ? zEps : -zEps)

    return h.zPositive[nr]
end

function negative!(h::EventHandler, nr::Int, crossing::Float64, crossingAsString::String;
                   restart::EventRestart=Restart)::Bool
    if h.initial
        h.zPositive[nr] = crossing >= 0.0
        if ModiaMath.isLogEvents(h.logger)
            println("        ", crossingAsString, " = ", crossing, negativeCrossingAsString(!h.zPositive[nr]))
        end

    elseif h.event
        new_zPositive = crossing >= 0.0
        change = (h.zPositive[nr] && !new_zPositive) || (!h.zPositive[nr] && new_zPositive)
        h.zPositive[nr] = new_zPositive

        if change
            h.restart = max(h.restart, restart)
            if ModiaMath.isLogEvents(h.logger)
                println("        ", crossingAsString, " = ", crossing, negativeCrossingAsString(!h.zPositive[nr]))
            end
            h.newEventIteration = true
        end
    end
    h.z[nr] = crossing + (h.zPositive[nr] ? zEps : -zEps)

    return !h.zPositive[nr]
end

function change!(h::EventHandler, nr::Int, crossing::Float64, crossingAsString::String;
                 restart::EventRestart=Restart)::Bool
    h.z[nr] = crossing
    if h.initial
        h.zPositive[nr] = crossing > 0.0
        if ModiaMath.isLogEvents(h.logger)
            println("        ", crossingAsString, " = ", crossing, positiveCrossingAsString(h.zPositive[nr]))
        end
        return false

    elseif h.event
        new_zPositive = crossing > 0.0
        change = (h.zPositive[nr] && !new_zPositive) || (!h.zPositive[nr] && new_zPositive)
        h.zPositive[nr] = new_zPositive

        if change
            h.restart = max(h.restart, restart)
            if ModiaMath.isLogEvents(h.logger)
                println("        ", crossingAsString, " = ", crossing, positiveCrossingAsString(h.zPositive[nr]))
            end
            h.newEventIteration = true
            return true
        end
    end

    return false
end

function edge!(h::EventHandler, nr::Int, crossing::Float64, crossingAsString::String;
               restart::EventRestart=Restart)::Bool
    h.z[nr] = crossing

    if h.initial
        h.zPositive[nr] = crossing > 0.0
        h.zDir[nr] = 1

        if ModiaMath.isLogEvents(h.logger)
            println("        ", crossingAsString, " = ", crossing, positiveCrossingAsString(h.zPositive[nr]))
        end
        return false

    elseif h.event
        new_zPositive = crossing > 0.0
        edge = !h.zPositive[nr] && new_zPositive
        h.zPositive[nr] = new_zPositive

        if edge
            h.restart = max(h.restart, restart)
            if ModiaMath.isLogEvents(h.logger)
                println("        ", crossingAsString, " = ", crossing, " (became > 0)")
            end
            h.newEventIteration = true
            return true
        end

    elseif h.crossing && h.zPositive[nr] && crossing < 0.0
        # Eventually need to improve later
        # (need to check, whether h.get_z is good enough to switch h.zPositive)
        h.zPositive[nr] = false
    end
    return false
end


text:
export QSFitResults, QSFitMultiResults

function estimate_line_EWs(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: AbstractRecipe
    EW = OrderedDict{Symbol, Float64}()

    cont = deepcopy(model())
    for (lname, lc) in pspec.lcs
        haskey(model, lname) || continue
        cont .-= model(lname)
    end
    for (lname, lc) in pspec.lcs
        haskey(model, lname) || continue
        EW[lname] = int_tabulated(domain(model)[:],
                                  model(lname) ./ cont)[1]
    end
    return EW
end

struct QSFitResults{T}
    source::QSO{T}
    pspec::PreparedSpectrum
    model::Model
    fitres::GFit.FitResult
    EW::OrderedDict{Symbol, Float64}
end

QSFitResults(source::QSO{T}, pspec::PreparedSpectrum, model::Model, fitres::GFit.FitResult) where T <: AbstractRecipe =
    QSFitResults{T}(source, pspec, model, fitres, estimate_line_EWs(source, pspec, model))

struct QSFitMultiResults{T}
    source::QSO{T}
    pspecs::Vector{PreparedSpectrum}
    multi::MultiModel
    fitres::GFit.FitResult
    EW::Vector{OrderedDict{Symbol, Float64}}
end

QSFitMultiResults(source::QSO{T}, pspecs::Vector{PreparedSpectrum}, multi::MultiModel, fitres::GFit.FitResult) where T <: AbstractRecipe =
    QSFitMultiResults{T}(source, pspecs, multi, fitres, [estimate_line_EWs(source, pspecs[id], multi[id]) for id in 1:length(multi)])


text:
module Gallery

using Pkg
import GeoPhyInv.IO
import GeoPhyInv.Models
import GeoPhyInv.Acquisition
import GeoPhyInv.FWI
import GeoPhyInv.Utils
using Statistics

global marmousi_folder=joinpath(@__DIR__, "marmousi2")

"""
Gallery of `Seismic` models.

# Arguments
* `attrib::Symbol` : name of the model

# Optional Arguments
* `δ::Real` : spatial sampling to resample the models

# Outputs
* `attrib=:acou_homo1` : an homogeneous acoustic model with `vp0=2000` and `ρ0=2000`
* `attrib=:acou_homo2` : same as above, but with spatial sampling as 40 m (faster testing)
* `attrib=:seismic_marmousi2` : marmousi model with lower resolution; ideal for surface seismic experiments
* `attrib=:seismic_marmousi2_high_res` : marmousi model high resolution; slower to load
* `attrib=:seismic_marmousi2_xwell` : boxed marmousi model ideal for crosswell experiments
* `attrib=:seismic_marmousi2_surf` : boxed marmousi2 for surface seismic experiments
* `attrib=:seismic_marmousi2_downhole` : boxed marmousi2 for downhole seismic experiments

"""
function Seismic(attrib::Symbol, δ::Real=0.0; verbose=false)
	bfrac=0.1;
	δ=Float64(δ)
	if((attrib == :acou_homo1))
		vp0 = [1500., 3500.] # bounds for vp
		vs0 = [1.0, 1.0] # dummy
		ρ0 = [1500., 3500.] # density bounds
		mgrid = repeat([range(-1000.0,stop=1000.0,length=201)],2)
		nz,nx=length.(mgrid)
		model= Models.Seismic(vp0, vs0, ρ0,
		      fill(0.0, (nz, nx)),
		      fill(0.0, (nz, nx)),
		      fill(0.0, (nz, nx)),
		      mgrid)

	elseif((attrib == :acou_homo2))
		vp0 = [1700., 2300.] # bounds for vp
		vs0 = [1.0, 1.0] # dummy
		ρ0 = [1700., 2300.] # density bounds
		mgrid = repeat([range(-1000.0,stop=1000.0,length=51)],2)
		nz,nx=length.(mgrid)
		model= Models.Seismic(vp0, vs0, ρ0,
		      fill(0.0, (nz, nx)),
		      fill(0.0, (nz, nx)),
		      fill(0.0, (nz, nx)),
		      mgrid)

	elseif(attrib == :seismic_marmousi2)
		vp, h= IO.readsu(joinpath(marmousi_folder,"vp_marmousi-ii_0.1.su"))
		vs, h= IO.readsu(joinpath(marmousi_folder,"vs_marmousi-ii_0.1.su"))
		ρ,  h= IO.readsu(joinpath(marmousi_folder,"density_marmousi-ii_0.1.su"))
		vp .*= 1000.; vs .*= 1000.; #ρ .*=1000
		vp0=Models.bounds(vp,bfrac);
		vs0=Models.bounds(vs,bfrac);
		ρ0=Models.bounds(ρ, bfrac);
		mgrid=[range(0.,stop=3500.,length=size(vp,1)),range(0., stop=17000., length=size(vp,2))]
		model= Models.Seismic(vp0, vs0, ρ0, Models.χ(vp,mean(vp0),1), Models.χ(vs,mean(vs0),1), Models.χ(ρ,mean(ρ0),1), mgrid)
	elseif(attrib == :seismic_marmousi2_high_res)
		vp, h= IO.readsu(joinpath(marmousi_folder,"vp_marmousi-ii.su"))
		vs, h= IO.readsu(joinpath(marmousi_folder,"vs_marmousi-ii.su"))
		ρ,  h= IO.readsu(joinpath(marmousi_folder,"density_marmousi-ii.su"))
		vp .*= 1000.; vs .*= 1000.; #ρ .*=1000
		vp0=Models.bounds(vp,bfrac);
		vs0=Models.bounds(vs,bfrac);
		ρ0=Models.bounds(ρ, bfrac);
		mgrid=[range(0.,stop=3500.,length=size(vp,1)),range(0., stop=17000., length=size(vp,2))]
		model= Models.Seismic(vp0, vs0, ρ0, Models.χ(vp,mean(vp0),1), Models.χ(vs,mean(vs0),1), Models.χ(ρ,mean(ρ0),1), mgrid)

	elseif(attrib == :seismic_marmousi2_xwell)
		model=Models.Seismic_trun(Seismic(:seismic_marmousi2_high_res),
				     zmin=1000., zmax=2000., xmin=8500., xmax=9500.,)
		Models.adjust_bounds!(model, bfrac) # adjuts bounds just inside the bounds
	elseif(attrib == :seismic_marmousi2_surf)
		model=Models.Seismic_trun(Seismic(:seismic_marmousi2_high_res),
				     xmin=6000., xmax=12000.,)
		Models.adjust_bounds!(model, bfrac) # adjust bounds just inside the bounds
	elseif(attrib == :seismic_marmousi2_downhole)
		model=Models.Seismic_trun(Seismic(:seismic_marmousi2_high_res),
				     xmin=9025., xmax=9125., zmin=1400., zmax=1600.,)
		Models.adjust_bounds!(model, bfrac) # adjust bounds just inside the bounds
	elseif(attrib == :seismic_marmousi2_rvsp)
		model=Models.Seismic_trun(Seismic(:seismic_marmousi2_high_res),
				     xmin=8000., xmax=10000., zmax=1700.,zmin=500.)
		Models.adjust_bounds!(model, bfrac) # adjust bounds just inside the bounds

	else
		error("invalid attrib")
	end
	if(δ==0.0)
		verbose && Models.print(model,string(attrib))
		return model
	elseif(δ > 0.0)
		mgrid_out=broadcast(x->range(x[1],stop=x[end],step=δ),model.mgrid)
		model_out=Models.Seismic_zeros(mgrid_out)
		Models.adjust_bounds!(model_out,model)
		Models.interp_spray!(model, model_out, :interp, :B1)
		verbose && Models.print(model_out,string(attrib))
		return model_out
	else
		error("invalid δ")
	end

end

"""
Gallery of acquisition geometries `Geom` using an input mesh.
The sources and receivers are not placed anywhere on the edges of the mesh.

# Arguments
* `mgrid` : a 2-D mesh
* `attrib::Symbol` : attribute decides output
  * `=:xwell` cross-well acquisition
  * `=:surf` surface acquisition
  * `=:vsp` vertical seismic profiling
  * `=:rvsp`  reverse vertical seismic profiling
  * `=:downhole` downhole sources and receivers

# Keyword Arguments

* `nss=2` : number of supersources
* `nr=2` : number of receivers per supersource
* `rand_flags::Vector{Bool}=[false, false]` : randomly or equally spaced supersources and receivers.
"""
function Geom(mgrid, attrib::Symbol; nss=2, nr=2, rand_flags=[false, false])
	otx=(0.9*mgrid[2][1]+0.1*mgrid[2][end]); ntx=(0.1*mgrid[2][1]+0.9*mgrid[2][end]);
	otwx=(0.95*mgrid[2][1]+0.05*mgrid[2][end]); ntwx=(0.05*mgrid[2][1]+0.95*mgrid[2][end]);
	otz=(0.9*mgrid[1][1]+0.1*mgrid[1][end]); ntz=(0.1*mgrid[1][1]+0.9*mgrid[1][end]);
	otwz=(0.95*mgrid[1][1]+0.05*mgrid[1][end]); ntwz=(0.05*mgrid[1][1]+0.95*mgrid[1][end]);
	quatx = (0.75*mgrid[2][1]+0.25*mgrid[2][end]); quatz = (0.75*mgrid[1][1]+0.25*mgrid[1][end])
	tquatx = (0.25*mgrid[2][1]+0.75*mgrid[2][end]); tquatz = (0.25*mgrid[1][1]+0.75*mgrid[1][end])
	halfx = 0.5*(mgrid[2][1]+mgrid[2][end]);	halfz = 0.5*(mgrid[1][1]+mgrid[1][end]);
	if(attrib == :xwell)
		geom=Acquisition.Geom_fixed(otz, ntz, otx, otz, ntz, ntx, nss, nr, :vertical, :vertical, rand_flags)
	elseif(attrib == :surf)
		geom=Acquisition.Geom_fixed(otx, ntx, otwz, otx, ntx, otwz, nss, nr, :horizontal, :horizontal, rand_flags)
	elseif(attrib == :vsp)
		geom=Acquisition.Geom_fixed(otx, ntx, otwz, quatz, ntz, otx, nss, nr, :horizontal, :vertical, rand_flags)
	elseif(attrib == :rvsp)
		geom=Acquisition.Geom_fixed(quatz, ntz, otx, otx, ntx, otwz, nss, nr, :vertical, :horizontal, rand_flags)
	elseif(attrib == :downhole)
		geom=Acquisition.Geom_fixed(quatz, otz, quatx, quatz, otz, quatx, nss, nr, :vertical, :vertical, rand_flags)
	else
		error("invalid attrib")
	end
	Acquisition.print(geom, string(attrib))
	return geom
end

"""
Gallery of source signals `Src`.

# Arguments
* `attrib::Symbol` :
* `nss::Int64=1` : number of supersources

# Outputs
* `attrib=:acou_homo1` :
"""
function Src(attrib::Symbol, nss::Int64=1)
	if(attrib == :acou_homo1)
		tgrid=range(0.0,stop=2.0,length=1000)
		wav = Utils.Wavelets.ricker(10.0, tgrid, tpeak=0.25, )
		return Acquisition.Src_fixed(nss, 1, [:P], wav, tgrid)
	elseif(attrib == :acou_homo2)
		tgrid=range(0.0,stop=2.0,length=250)
		wav = Utils.Wavelets.ricker(3.0, tgrid, tpeak=0.3, )
		return Acquisition.Src_fixed(nss, 1, [:P], wav, tgrid)
	elseif(attrib == :vecacou_homo1)
		tgrid=range(0.0,stop=2.0,length=1000)
		wav = Utils.Wavelets.ricker(10.0, tgrid, tpeak=0.25, )
		return Acquisition.Src_fixed(nss, 1, [:P, :Vx, :Vz], wav, tgrid)
	end
end

for file in ["fwi.jl", "fdtd.jl"]
	fn=joinpath(@__DIR__, file)
	include(fn)
	#Revise.track(@__MODULE__,fn)
end

end # module


text:
module Folly

using KernelAbstractions: @kernel, @index, CPU, CUDA
export CPU, CUDA

export Outer, Inner, Block
export IndVar, Loop
export Native, KATarget

abstract type Kind end

# @index(Group, NTuple)[dim]

# Assumed to be parallel
struct Outer <: Kind end
struct Inner <: Kind end

# @index(Local, NTuple)[dim]
struct Block <: Kind end

# Recurrence relationship
# IndVar{Outer}(:gI)
# IndVar{Block}(:lI)
# IndVar{Dependent}(:I, ..., :(gI * N + lI))
# struct Dependent <: Kind end

# TODO: add metadata like unroll factor,
# The induction variable is most often an index
struct IndVar{Kind}
    name::Symbol
    kind::Kind
    recurrence::Expr
    lowerbound::Union{Number, Expr}
    upperbound::Union{Number, Expr}
end

# Represents an loop over a variable
# the body can contain statements and/or other loops
struct Loop
    indvar::Symbol
    body::Vector{Any}
end

# TODO:
# - From recurrence to affine-map
# - Dependent IndVars
# - Compilation to KernelAbstractions
#  - map IndVars to NDrange dimensions

abstract type Target end
struct KATarget <: Target
    device
    dims::Dict{IndVar, Int}
    KATarget(device) = new(device, Dict{IndVar, Int}())
end
struct Native <: Target end

construct(::Target, expr::Expr, indvars) = expr
function construct(target::Target, loop::Loop, indvars)
    idx = findfirst(indv->indv.name == loop.indvar, indvars)
    @assert idx !== nothing
    indvar = indvars[idx]
    construct(target, indvar, loop, indvars)
end

function construct(target::Native, indvar::IndVar, loop::Loop, indvars)
    body = map(expr->construct(target, expr, indvars), loop.body)
    name = indvar.name
    quote
        $name = $(indvar.lowerbound)
        while $name <= $(indvar.upperbound)
            $(body...)
            $name = $(indvar.recurrence)
        end
    end
end

function construct(target::KATarget, indvar::IndVar{Inner}, loop::Loop, indvars)
    body = map(expr->construct(target, expr, indvars), loop.body)
    name = indvar.name
    quote
        $name = $(indvar.lowerbound)
        while $name <= $(indvar.upperbound)
            $(body...)
        end
    end
end

index_expr(::IndVar{Outer}, dim) = :(@index(Group, NTuple)[$dim])
index_expr(::IndVar{Block}, dim) = :(@index(Local, NTuple)[$dim])

function construct(target::KATarget, indvar::IndVar, loop::Loop, indvars)
    body = map(expr->construct(target, expr, indvars), loop.body)
    name = indvar.name
    dim = target.dims[indvar]
    quote
        $name = $(index_expr(indvar, dim))
        if $(indvar.lowerbound) <= $name <= $(indvar.upperbound)
            $(body...)
        end
    end
end

function compile(target::Native, loop::Loop, indvars, freevars)
    body = construct(target, loop, indvars)
    expr = quote
        function $(gensym(:folly))($(freevars...))
            $(body)
        end
    end
    eval(expr)
end

# TODO: Sort indvars according to depth
function compile(target::KATarget, loop::Loop, indvars, freevars)
    outer_indvars = filter(iv -> iv.kind isa Outer, indvars)
    block_indvars = filter(iv -> iv.kind isa Block, indvars)

    @assert length(block_indvars) == 0
    dim = 1
    for iv in block_indvars
        target.dims[iv] = dim
        dim += 1
    end
    for iv in outer_indvars
        target.dims[iv] = dim
        dim += 1
    end
    body = construct(target, loop, indvars)

    workgroupsize = append!(
                        map(iv->iv.upperbound, block_indvars),
                        map(iv->1, outer_indvars))
    ndrange       = append!(
                        map(iv->iv.upperbound, block_indvars),
                        map(iv->iv.upperbound, outer_indvars))
    ka_func = quote
    end

    expr = quote
        function $(gensym(:folly))($(freevars...); dependencies=nothing)
            @kernel function kernel($(freevars...))
                $(body)
            end
            workgroupsize = ($(workgroupsize...),)
            ndrange       = ($(ndrange...),)

            k = kernel($(target.device), workgroupsize)
            k($(freevars...), ndrange=ndrange, dependencies=dependencies)
        end

    end
    eval(expr)
end

function freevars(loop::Loop, indvars::Set{Symbol} = Set{Symbol}(), free::Set{Symbol} = Set{Symbol}())
    @assert loop.indvar ∉ indvars
    push!(indvars, loop.indvar)
    for expr in loop.body
        freevars(expr, indvars, free)
    end
    delete!(indvars, loop.indvar)
    return free
end

function freevars(name::Symbol, indvars::Set{Symbol}, free::Set{Symbol})
    if name ∉ indvars
        push!(free, name)
    end
    return free
end

function freevars(expr::Expr, indvars::Set{Symbol}, free::Set{Symbol})
    if expr.head === :ref || expr.head === :(=) || expr.head === :(+=)
        for arg in expr.args
            freevars(arg, indvars, free)
        end
    elseif expr.head === :call
        # skip first arg
        for i in 2:length(expr.args)
            freevars(expr.args[i], indvars, free)
        end
    else
        @warn "Unable to handle expr" expr.head
    end
    return free
end

# Very simple fusion, we expect that the indvars match
# this does not check that the loops can be fused.
# i.e. reductions and cross loop dependencies can be broken
function fuse!(loop::Loop, other::Loop)
    if loop.indvar == other.indvar
        for expr in other.body
            if expr isa Loop
                n = expr.name
                idxs = findall(e -> e isa Loop && e.name == n, loop.body)
                if length(idxs) == 1
                    fuse(loop.body[idxs[1]], expr)
                    continue
                elseif length(idxs) > 1
                    error("IndVar not unique on loop-level")
                end
            end
            push!(loop.body, expr)
        end
    end
    return nothing
end

end # module


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

function main()
    # Define beam
    λ = LinRange(200,5000,1000) # wavelength range [nm]
    θ = LinRange(0,90,90)
    beam = PlaneWave(λ,θ)
    # Define layers
    l0 = LayerTMMO(RIdb.air(beam.λ))
    l1 = LayerTMMO(RIdb.dummy(beam.λ,1.45,0.0); type=:OT, d=1/4.)
    l2 = LayerTMMO(RIdb.dummy(beam.λ,3.45,0.0); type=:OT, d=1/4.)
    layers = [l0, l1, l2, l1, l2, l1, l2, l1, l2, l0]
    # Reference wavelenth
    λ0 = 700.0
    # call main script
    sol = tmm_optics(beam, layers; λ0=λ0, emfflag=true, h=10, pbgflag=true)
    return sol
end

function plotEMF(sol, ϕ)
    ϕ_ = findmin(abs.(sol.beam.θ .- ϕ))[2][1]
    plot(EMF2D(),
              sol.beam.λ, sol.Misc.ℓ, log10.(sol.Field.emfp[:,ϕ_,:]),
              title="Log EMF intesnsity, p-TM",
    )
    gui()
    plot(EMF2D(),
              sol.beam.λ, sol.Misc.ℓ, log10.(sol.Field.emfs[:,ϕ_,:]),
              title="Log EMF intesnsity, s-TE",
    )
    gui()
    return nothing
end

sol = main()

# plot the R, T and A spectra
plot(Spectrum2D(),
        sol.beam.λ, sol.beam.θ, sol.Spectra.Rp,
        title="Reflectance, p-TM", clims=(0.0, 1.0),
)
gui()
plot(Spectrum2D(),
        sol.beam.λ, sol.beam.θ, sol.Spectra.Rs,
        title="Reflectance, s-TE", clims=(0.0, 1.0),
)
gui()

# Plot the EMF at normal incidence
plotEMF(sol, 0.0)

# Plot the EMF at 15 degrees
plotEMF(sol, 15.0)

# plot the refractive index profile
plot(RIprofile(), sol; plotemf=true)
gui()
### Note that if you choose θ outside the range it will show the EMF for one extrema. Same for λ.

# plot the photonic dispersion with custom function
plot(PBGDispersion2D(), sol.Bloch; wave=:p)
gui()

plot(PBGDispersion2D(), sol.Bloch; wave=:s)
gui()

plot(PBGDispersion2Dalt(), sol.Bloch)
gui()


text:
"""
    IndexedGraph{T<:Integer} <: AbstractIndexedGraph{T}

A type representing a sparse undirected graph.

### FIELDS

- `A` -- square adjacency matrix. `A[i,j] == A[j,i]` contains the unique index associated to unidrected edge `(i,j)`
"""
struct IndexedGraph{T<:Integer} <: AbstractIndexedGraph{T}
    A :: SparseMatrixCSC{T, T}   # Adjacency matrix. Values are unique edge id's

    function IndexedGraph(A::SparseMatrixCSC{T}) where T
        _checksquare(A)
        _check_selfloops(A)
        M = SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, ones(Int, nnz(A)))
        idx_map = inverse_edge_idx(M)
        cnt = 1
        # give increasing indices to upper triangular elements
        for j in 1:size(M, 2)
            for k in nzrange(M, j)
                i = M.rowval[k]
                if i < j
                    M.nzval[k] = cnt; cnt += 1
                else
                    break
                end
            end
        end
        @assert cnt - 1 == (nnz(M) / 2)  # sanity check
        # give (i=>j) the same idx as (j=>i)
        for j in 1:size(M, 2)
            for k in nzrange(M, j)
                i = M.rowval[k]
                if i > j
                    M.nzval[k] = M.nzval[idx_map[k]]
                end
            end
        end
        new{Int}(M)
    end
end

"""
    IndexedGraph(A::AbstractMatrix)

Construct an `IndexedGraph` from symmetric adjacency matrix A.
"""
IndexedGraph(A::AbstractMatrix) = IndexedGraph(convert(SparseMatrixCSC, A))

"""
    IndexedGraph(A::SimpleGraph)

Construct an `IndexedGraph` from undirected `SimpleGraph` (Graphs.jl).
"""
IndexedGraph(sg::SimpleGraph) = IndexedGraph(adjacency_matrix(sg))

function edges(g::IndexedGraph)
    (IndexedEdge{Int}(extrema((i, g.A.rowval[k]))..., g.A.nzval[k])
        for i=1:size(g.A,2) for k=nzrange(g.A,i) if i > g.A.rowval[k])
end

ne(g::IndexedGraph) = Int( nnz(g.A) / 2 )

neighbors(g::IndexedGraph, i::Integer) = outneighbors(g, i)

inneighbors(g::IndexedGraph, i::Integer) = outneighbors(g, i)

is_directed(g::IndexedGraph) = false

is_directed(::Type{IndexedGraph{T}}) where T = false

Base.zero(g::IndexedGraph) = IndexedGraph(zero(g.A))

"""
    edges(g::IndexedGraph, i::Integer)

Return a lazy iterators to the edges incident to `i`.

By default unordered edges sort source and destination nodes in increasing order.
See [outedges](@ref) and [inedges](@ref) if you need otherwise.
"""
function edges(g::IndexedGraph, i::Integer)
    (IndexedEdge(extrema((i, g.A.rowval[k]))..., g.A.nzval[k])
        for k in nzrange(g.A, i))
end

"""
    outedges(g::IndexedGraph, i::Integer)

Return a lazy iterators to the edges incident to `i` with `i` as the source.
"""
function outedges(g::IndexedGraph, i::Integer)
    (IndexedEdge(i, g.A.rowval[k], g.A.nzval[k]) for k in nzrange(g.A, i))
end

"""
    inedges(g::IndexedGraph, i::Integer)

Return a lazy iterators to the edges incident to `i` with `i` as the destination.
"""
function inedges(g::IndexedGraph, i::Integer)
    (IndexedEdge(g.A.rowval[k], i, g.A.nzval[k]) for k in nzrange(g.A, i))
end

function edge_idx(g::IndexedGraph, src::Integer, dst::Integer)
    k = nzindex(g.A, src, dst)
    g.A.nzval[k]
end
function edge_src_dst(g::IndexedGraph, id::Integer)
    k = findfirst(isequal(id), g.A.nzval)
    i, j = nzindex(g.A, k)
    return extrema((i,j))    # return sorted
end

"""
    get_edge(g::IndexedGraph, src::Integer, dst::Integer)
    get_edge(g::IndexedGraph, id::Integer)

Get edge given source and destination or given edge index.
"""
function get_edge(g::IndexedGraph, src::Integer, dst::Integer)
    id = edge_idx(g, src, dst)
    IndexedEdge(src, dst, id)
end

function get_edge(g::IndexedGraph, id::Integer)
    i, j = edge_src_dst(g, id)
    IndexedEdge(i, j, id)
end


text:
# Autogenerated wrapper script for ITSOL_2_jll for aarch64-linux-gnu-libgfortran5
export libITSOL_2

using LAPACK_jll
using OpenBLAS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ITSOL_2")
JLLWrappers.@declare_library_product(libITSOL_2, "libITSOL_2.so")
function __init__()
    JLLWrappers.@generate_init_header(LAPACK_jll, OpenBLAS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libITSOL_2,
        "lib/libITSOL_2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"Given the previous scattering interface and current layer information, return what type of scattering interface is nexts"
function get_scattering_interface(scattering_interface, scatter, iz)

    # First layer (TOA)
    if (iz == 1)

        # If scattering, 4. If non-scattering, 1.
        scattering_interface = scatter ? ScatteringInterface_11() : ScatteringInterface_00()

    # Not the first layer (not TOA)
    elseif !(scattering_interface isa ScatteringInterface_00)

        # If kn was 1, then toggle between 0/0 and 0/1
        # Else, toggle between 1/0 and 1/1
        scattering_interface = (scattering_interface isa ScatteringInterface_00) ?
                                    (!scatter ? ScatteringInterface_00() : ScatteringInterface_01()) :
                                    (!scatter ? ScatteringInterface_10() : ScatteringInterface_11())
    end

    return scattering_interface
end

"Minimum number of doublings needed to reach an optical depth τ_end, starting with an optical depth dτ.
The starting optical depth dτ is also determined from its maximum possible value, τ"
function doubling_number(dτ_max, τ_end) # check if τ_end can be replaced by τ_end*ϖ for absorbing atmospheres
    FT = eltype(dτ_max)

    # minimum number of doublings needed to reach an optical depth τ_end, starting with an optical depth dτ.
    # The starting optical depth dτ is also determined from its maximum possible value, dτ_max
    if τ_end <= dτ_max
        dτ = τ_end
        ndoubl = 0
        return dτ, ndoubl
    else
        q1 = log10(2.0)
        q2 = log10(dτ_max)
        q3 = log10(τ_end)
        tlimit = (q3 - q2) / q1
        nlimit = floor(Int, tlimit)
        diff = tlimit - nlimit
        if diff < eps(FT)
            dτ = dτ_max
            ndoubl = nlimit
        else
            ndoubl = nlimit + 1
            x = q3 - q1 * ndoubl
            dτ = 10.0^x
        end
        return dτ, ndoubl
    end
end

"Finds index i of f_array (i) which is nearest point to f"
nearest_point(f_array, f) = argmin(abs.(f_array.-f))

"Get indices scaled according to pol_type"
function get_indices(iμ::Integer, pol_type::AbstractPolarizationType)

    st_iμ = (iμ - 1) * pol_type.n
    istart = st_iμ + 1
    iend   = st_iμ + pol_type.n

    return st_iμ, istart, iend
end

"Default matrix in RT calculation (zeros)"
default_matrix(FT, arr_type, dims, nSpec)   = arr_type(zeros(FT, tuple(dims[1], dims[2], nSpec)))

"Default J matrix in RT calculation (zeros)"
default_J_matrix(FT, arr_type, dims, nSpec) = arr_type(zeros(FT, tuple(dims[1], 1, nSpec)))

"Default matrix in RT calculation (random)"
default_matrix_rand(FT, arr_type, dims, nSpec)   = arr_type(randn(FT, tuple(dims[1], dims[2], nSpec)))

"Default J matrix in RT calculation (random)"
default_J_matrix_rand(FT, arr_type, dims, nSpec) = arr_type(randn(FT, tuple(dims[1], 1, nSpec)))

"Make an added layer, supplying all default matrices"
make_added_layer(FT, arr_type, dims, nSpec)     = AddedLayer(default_matrix(FT, arr_type, dims, nSpec),
                                                         default_matrix(FT, arr_type, dims, nSpec),
                                                         default_matrix(FT, arr_type, dims, nSpec),
                                                         default_matrix(FT, arr_type, dims, nSpec),
                                                         default_J_matrix(FT, arr_type, dims, nSpec),
                                                         default_J_matrix(FT, arr_type, dims, nSpec)
                                                         )

"Make a random added layer, supplying all random matrices"
make_added_layer_rand(FT, arr_type, dims, nSpec)     = AddedLayer(default_matrix_rand(FT, arr_type, dims, nSpec),
                                                         default_matrix_rand(FT, arr_type, dims, nSpec),
                                                         default_matrix_rand(FT, arr_type, dims, nSpec),
                                                         default_matrix_rand(FT, arr_type, dims, nSpec),
                                                         default_J_matrix_rand(FT, arr_type, dims, nSpec),
                                                         default_J_matrix_rand(FT, arr_type, dims, nSpec)
                                                         )

"Make a composite layer, supplying all default matrices"
make_composite_layer(FT, arr_type, dims, nSpec) = CompositeLayer(default_matrix(FT, arr_type, dims, nSpec),
                                                                 default_matrix(FT, arr_type, dims, nSpec),
                                                                 default_matrix(FT, arr_type, dims, nSpec),
                                                                 default_matrix(FT, arr_type, dims, nSpec),
                                                                 default_J_matrix(FT, arr_type, dims, nSpec),
                                                                 default_J_matrix(FT, arr_type, dims, nSpec)
                                                            )

"Given a ComputedAtmosphereProperties object, extract a ComputedLayerProperties object using data from the iz index of all arrays in the ComputedAtmosphereProperties"
function get_layer_properties(computed_atmospheric_properties::ComputedAtmosphereProperties, iz, arr_type)

    @unpack τ_λ_all, ϖ_λ_all, τ_all, ϖ_all, Z⁺⁺_all, Z⁻⁺_all, dτ_max_all, dτ_all, ndoubl_all, dτ_λ_all, expk_all, scatter_all, τ_sum_all, scattering_interfaces_all = computed_atmospheric_properties

    τ_λ = arr_type(τ_λ_all[:, iz])
    ϖ_λ = arr_type(ϖ_λ_all[:, iz])
    τ   = τ_all[iz]
    ϖ   = ϖ_all[iz]
    Z⁺⁺ = arr_type(Z⁺⁺_all[:,:,iz])
    Z⁻⁺ = arr_type(Z⁻⁺_all[:,:,iz])

    dτ_max = dτ_max_all[iz]
    dτ     = dτ_all[iz]
    ndoubl = ndoubl_all[iz]
    dτ_λ   = arr_type(dτ_λ_all[:, iz])
    expk   = arr_type(expk_all[:, iz])
    scatter = scatter_all[iz]
    τ_sum = arr_type(τ_sum_all[:,iz])
    scattering_interface = scattering_interfaces_all[iz]

    # τ * ϖ should remain constant even though they individually change over wavelength
    # @assert all(i -> (i ≈ τ * ϖ), τ_λ .* ϖ_λ)

    return ComputedLayerProperties(τ_λ, ϖ_λ, τ, ϖ, Z⁺⁺, Z⁻⁺, dτ_max, dτ, ndoubl, dτ_λ, expk, scatter, τ_sum, scattering_interface)
end


text:
#InspectDR: Support for pan/zoom
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const PAN_STEPRATIO = 0.25 #Percentage of current extents
const ZOOM_STEPRATIO = 2.0 #How much to zoom in/out with mousewheel + keybindings

#==Types
===============================================================================#
#User input states
mutable struct ISPanningData <: InputState #Typically with mouse
	p0::Point2D #Start
	pmouse::Point2D #move
	ext_start::PExtents2D #Extents @ start of operation
	#Store ext_start to avoid accumulation of numerical errors.
	istrip::Int
end
mutable struct ISSelectingArea <: InputState #Box-zoom
	bb::BoundingBox
	istrip::Int
end

#==Drawing functions
===============================================================================#
function selectionbox_draw(ctx::CairoContext, selbb::BoundingBox, graphbb::BoundingBox, hallowed::Bool, vallowed::Bool)
	xmin = selbb.xmin; xmax = selbb.xmax
	ymin = selbb.ymin; ymax = selbb.ymax
	if !hallowed
		xmin = graphbb.xmin
		xmax = graphbb.xmax
	end
	if !vallowed
		ymin = graphbb.ymin
		ymax = graphbb.ymax
	end

	Cairo.save(ctx) #-----
	setlinestyle(ctx, LineStyle(:dash, 1.0, COLOR_BLACK))
	Cairo.rectangle(ctx, BoundingBox(xmin, xmax, ymin, ymax))
	Cairo.stroke(ctx)
	Cairo.restore(ctx) #-----

	nothing
end

function selectionbox_draw(ctx::CairoContext, w::PlotWidget)
	if !isa(w.state, ISSelectingArea); return; end
	graphbb = w.plotinfo.strips[activestrip(w)].graphbb
	selectionbox_draw(ctx, w.state.bb, graphbb, w.hallowed, w.vallowed)
end

#==Switching to h/v lock (or releasing both)
===============================================================================#
function locdir_h(pwidget::PlotWidget)
	pwidget.hallowed = true
	pwidget.vallowed = false
	refresh(pwidget, refreshdata=false)
end
function locdir_v(pwidget::PlotWidget)
	pwidget.hallowed = false
	pwidget.vallowed = true
	refresh(pwidget, refreshdata=false)
end
function locdir_any(pwidget::PlotWidget)
	pwidget.hallowed = true
	pwidget.vallowed = true
	refresh(pwidget, refreshdata=false)
end

#==Scale/position widget control
===============================================================================#
function scalectrl_enabled(pwidget::PlotWidget)
	return Gtk.GAccessor.sensitive(pwidget.w_xscale)
end
function scalectrl_enabled(pwidget::PlotWidget, v::Bool)
	Gtk.GAccessor.sensitive(pwidget.w_xscale, v)
	Gtk.GAccessor.sensitive(pwidget.w_xpos, v)
end
#Apply current scale/position scrollbar values to plot extents:
function scalectrl_apply(pwidget::PlotWidget)
	xscale = get_gtk_property(pwidget.xscale, :value, Int)
	xpos = get_gtk_property(pwidget.xpos, :value, Float64)

	plot = pwidget.src #WANTCONST
	ixf = InputXfrm1D(plot.xscale) #WANTCONST

	#Use transformed coordinate system:
	xext_full = read2axis(getxextents_full(plot), ixf)
	span = xext_full.max - xext_full.min
	center = (xext_full.max + xext_full.min) / 2
	vspan = span/xscale #Visible span
	xmin = center + span*xpos - vspan/2
	xmax = xmin + vspan

	#Update extents & redraw
	xext_new = PExtents1D(xmin, xmax)
	xext = merge(getxextents_axis(plot), xext_new)
	setxextents_axis(plot, xext)

	refresh(pwidget, refreshdata=true)
end

#==Basic zoom control
===============================================================================#
#Zoom to bounding box (in device coordinates):
function zoom(pwidget::PlotWidget, bb::BoundingBox, istrip::Int)
	p1 = Point2D(bb.xmin, bb.ymin)
	p2 = Point2D(bb.xmax, bb.ymax)

	xf = Transform2D(pwidget.plotinfo, istrip)
	p1 = map2axis(xf, p1)
	p2 = map2axis(xf, p2)

	if pwidget.hallowed
		setxextents_axis(pwidget.src, PExtents1D(min(p1.x, p2.x), max(p1.x, p2.x)))
	end

	if pwidget.vallowed
		setyextents_axis(pwidget.src, PExtents1D(min(p1.y, p2.y), max(p1.y, p2.y)), istrip)
	end

	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	refresh(pwidget, refreshdata=true)
end

#Zoom in/out @ point (pt in plot coordinates)
function zoom(pwidget::PlotWidget, ext::PExtents2D, pt::Point2D, ratio::Float64, istrip::Int)
	xspan = ext.xmax - ext.xmin
	yspan = ext.ymax - ext.ymin
	Δx = pt.x - ext.xmin
	Δy = pt.y - ext.ymin
	xmin = pt.x-ratio*Δx
	ymin = pt.y-ratio*Δy

	setxextents_axis(pwidget.src, PExtents1D(xmin, xmin + ratio*xspan))
	setyextents_axis(pwidget.src, PExtents1D(ymin, ymin + ratio*yspan), istrip)

	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	refresh(pwidget, refreshdata=true)
end

#Zoom in/out, centered on current extents
function zoom(pwidget::PlotWidget, ratio::Float64, istrip::Int)
	ext = getextents_axis(pwidget.src, istrip)
	pt = Point2D((ext.xmin+ext.xmax)/2, (ext.ymin+ext.ymax)/2)
	zoom(pwidget, ext, pt, ratio, istrip)
end
zoom_out(pwidget::PlotWidget, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, stepratio, activestrip(pwidget))
zoom_in(pwidget::PlotWidget, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, 1/stepratio, activestrip(pwidget))

#Zoom in/out around specified device coordinates:
function zoom(pwidget::PlotWidget, x::Float64, y::Float64, ratio::Float64, istrip::Int)
	pt = Point2D(x, y)
	ext = getextents_axis(pwidget.src, istrip)
	xf = Transform2D(pwidget.plotinfo, istrip)
	pt = map2axis(xf, pt)
	zoom(pwidget, ext, pt, ratio, istrip)
end
zoom_out(pwidget::PlotWidget, x::Float64, y::Float64, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, x, y, stepratio, activestrip(pwidget))
zoom_in(pwidget::PlotWidget, x::Float64, y::Float64, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, x, y, 1/stepratio, activestrip(pwidget))

function zoom_xfull(pwidget::PlotWidget)
	#TODO
end
function zoom_full(pwidget::PlotWidget, hallowed::Bool, vallowed::Bool, istrip::Int)
	#Reset desired extents:
	if hallowed
		setxextents(pwidget.src, PExtents1D())
	end
	if vallowed
		setyextents(pwidget.src, PExtents1D(), istrip)
	end

	if hallowed
		scalectrl_enabled(pwidget, false) #Suppress updates from set_gtk_property!
		set_gtk_property!(pwidget.xscale, :value, Int(1))
		set_gtk_property!(pwidget.xpos, :value, Float64(0))
		scalectrl_enabled(pwidget, true)
		scalectrl_apply(pwidget)
		#Will refresh
	else
		refresh(pwidget, refreshdata=true)
	end
end
zoom_full(pwidget::PlotWidget, hallowed::Bool=true, vallowed::Bool=true) =
	zoom_full(pwidget, hallowed, vallowed, activestrip(pwidget))
zoom_hfull(pwidget::PlotWidget) = zoom_full(pwidget, true, false)
zoom_vfull(pwidget::PlotWidget) = zoom_full(pwidget, false, true)

#==Box-zoom control
===============================================================================#
function boxzoom_setstart(pwidget::PlotWidget, x::Float64, y::Float64)
	locdir_any(pwidget)
	bb = BoundingBox(x, x, y, y)
	istrip = activestrip(pwidget)
	gdk_window_set_cursor(pwidget.canvas, CURSOR_BOXSELECT)
	pwidget.state = ISSelectingArea(bb, istrip)
	return
end
function boxzoom_cancel(pwidget::PlotWidget)
	setstate_normal(pwidget)
	Gtk.draw(pwidget.canvas)
	return
end
function boxzoom_complete(s::ISSelectingArea, pwidget::PlotWidget, x::Float64, y::Float64)
	setstate_normal(pwidget)
	bb = BoundingBox(s.bb.xmin, x, s.bb.ymin, y)
	zoom(pwidget, bb, s.istrip)
	return
end
#Change end point of boxzoom area:
function boxzoom_setend(s::ISSelectingArea, pwidget::PlotWidget, x::Float64, y::Float64)
	s.bb = BoundingBox(s.bb.xmin, x, s.bb.ymin, y)
	Gtk.draw(pwidget.canvas)
end

#==Basic pan control
===============================================================================#
function pan_xratio(pwidget::PlotWidget, panstepratio::Float64)
	xext = getxextents_axis(pwidget.src)
	panstep = panstepratio*(xext.max-xext.min)
	setxextents_axis(pwidget.src,
		PExtents1D(xext.min+panstep, xext.max+panstep))
	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	refresh(pwidget, refreshdata=true)
end
function pan_yratio(pwidget::PlotWidget, panstepratio::Float64, istrip::Int)
	yext = getyextents_axis(pwidget.src, istrip)
	panstep = panstepratio*(yext.max-yext.min)
	setyextents_axis(pwidget.src,
		PExtents1D(yext.min+panstep, yext.max+panstep), istrip)
	refresh(pwidget, refreshdata=true)
end
pan_yratio(pwidget::PlotWidget, panstepratio::Float64) =
	pan_yratio(pwidget, panstepratio, activestrip(pwidget))

pan_left(pwidget::PlotWidget) = pan_xratio(pwidget, -PAN_STEPRATIO)
pan_right(pwidget::PlotWidget) = pan_xratio(pwidget, PAN_STEPRATIO)
pan_up(pwidget::PlotWidget) = pan_yratio(pwidget, PAN_STEPRATIO)
pan_down(pwidget::PlotWidget) = pan_yratio(pwidget, -PAN_STEPRATIO)

#==Mouse-pan control
===============================================================================#
#Δy/Δy: in device coordinates
function mousepan_delta(pwidget::PlotWidget, ext::PExtents2D, Δx::Float64, Δy::Float64, istrip::Int)
	#Convert to plot coordinates:
	xf = Transform2D(ext, pwidget.plotinfo.strips[istrip].graphbb)
	Δvec = map2axis(xf, Vector2D(-Δx, -Δy))

	setextents_axis(pwidget.src, ext, istrip) #Restore original extents before overwriting

	if pwidget.hallowed
		setxextents_axis(pwidget.src, PExtents1D(ext.xmin+Δvec.x, ext.xmax+Δvec.x))
	end
	if pwidget.vallowed
		setyextents_axis(pwidget.src, PExtents1D(ext.ymin+Δvec.y, ext.ymax+Δvec.y), istrip)
	end

	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	refresh(pwidget, refreshdata=true)
end
function mousepan_setstart(pwidget::PlotWidget, x::Float64, y::Float64)
	locdir_any(pwidget)
	istrip = activestrip(pwidget)
	ext_start = getextents_axis(pwidget.src, istrip)
	p0 = Point2D(x, y)
	gdk_window_set_cursor(pwidget.canvas, CURSOR_PAN)
	pwidget.state = ISPanningData(p0, p0, ext_start, istrip)
	return
end
function mousepan_cancel(s::ISPanningData, pwidget::PlotWidget)
	mousepan_delta(pwidget, s.ext_start, 0.0, 0.0, s.istrip)
	setstate_normal(pwidget)
	return
end
function mousepan_complete(pwidget::PlotWidget, x::Float64, y::Float64)
	#Already panned.  Just revert to normal state:
	setstate_normal(pwidget)
	return
end
#Set new point of mousepan operation:
function mousepan_move(s::ISPanningData, pwidget::PlotWidget, x::Float64, y::Float64)
	s.pmouse = Point2D(x, y)
	Δx = s.pmouse.x-s.p0.x; Δy = s.pmouse.y-s.p0.y
	mousepan_delta(pwidget, s.ext_start, Δx, Δy, s.istrip)
end

#==Event handlers for state: ISPanningData (Panning with mouse)
===============================================================================#
function handleevent_keypress(s::ISPanningData, pwidget::PlotWidget, event::Gtk.GdkEventKey)
	if GdkKeySyms.Escape == event.keyval
		mousepan_cancel(s, pwidget)
	elseif Int('h') == event.keyval || Int('H') == event.keyval
		locdir_h(pwidget)
	elseif Int('v') == event.keyval || Int('V') == event.keyval
		locdir_v(pwidget)
	elseif Int('b') == event.keyval || Int('B') == event.keyval
		locdir_any(pwidget)
	end
end
function handleevent_mouserelease(::ISPanningData, pwidget::PlotWidget, event::Gtk.GdkEventButton)
	if 1==event.button
		mousepan_complete(pwidget, event.x, event.y)
	end
end
function handleevent_mousemove(s::ISPanningData, pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	handleevent_plothover(pwidget, event)
	mousepan_move(s, pwidget, event.x, event.y)
end

#==Event handlers for state: ISSelectingArea (for box-zoom)
===============================================================================#
function handleevent_keypress(::ISSelectingArea, pwidget::PlotWidget, event::Gtk.GdkEventKey)
	if GdkKeySyms.Escape == event.keyval
		boxzoom_cancel(pwidget)
	elseif Int('h') == event.keyval
		locdir_h(pwidget)
	elseif Int('v') == event.keyval
		locdir_v(pwidget)
	elseif Int('b') == event.keyval
		locdir_any(pwidget)
	end
end
function handleevent_mouserelease(s::ISSelectingArea, pwidget::PlotWidget, event::Gtk.GdkEventButton)
	if 3==event.button
		boxzoom_complete(s, pwidget, event.x, event.y)
	end
end
function handleevent_mousemove(s::ISSelectingArea, pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	handleevent_plothover(pwidget, event)
	boxzoom_setend(s, pwidget, event.x, event.y)
end

#Last line


text:
include("../src/MsgPackRpcClient.jl")

using MsgPackRpcClient

session = MsgPackRpcClientSession.create()

the_number_of_call = 10000
results = {}

@sync begin
@async for i in 1:the_number_of_call
  future = call(session, "hello"; sync = false)
  push!(results, get(future))
end
end

@sync begin
@async for r in results
  println(r)
end
end

session.destroy()

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia


text:
module Domains

using  Glob, JSON, Proj4
import Base: show

export Domain,
       lonlat2lcc,
       lcc2lonlat,
       in,
       getgridpoints

const moduledir = @__DIR__
const jsondir = joinpath(moduledir, "json")
const jsonschemafile = joinpath(moduledir, "jsonschema/domain.schema.json")
const domains = getindex.(splitext.(basename.(Glob.glob("*.json", jsondir))), 1) # Array with domain names
const Rearth = 6.37122e6

"""
   Domain definitions
"""
struct Domain
    NAME::String
    TSTEP::Int
    NLON::Int
    NLAT::Int
    LONC::Float64
    LATC::Float64
    LON0::Float64
    LAT0::Float64
    GSIZE::Float64
    EZONE::Int
end

"""
    Domain(d::Dict)

Constructor for domains from Dictionaries
"""
Domain(d::Dict) = Domain(
       d["NAME"],d["TSTEP"],d["NLON"], d["NLAT"], d["LONC"],
       d["LATC"], d["LON0"], d["LAT0"],d["GSIZE"], d["EZONE"]
    )

"""
    Domain(name::String)

Constructor for `Domain`
"""
Domain(name::String) =  Domain(JSON.parsefile(joinpath(jsondir, "$name.json")))

function Base.show(io::IO,::MIME"text/plain",d::Domain)
    println(io, "NAME = $(d.NAME)")
    println(io, "NLON = $(d.NLON)")
    println(io, "NLAT = $(d.NLAT)")
    println(io, "LON0 = $(d.LON0)")
    println(io, "LAT0 = $(d.LAT0)")
    println(io, "LONC = $(d.LONC)")
    println(io, "LATC = $(d.LATC)")

    #println(io, ), d.LON0, d.LAT0])

end

"""
    Plcc(d)

Returns a Lambert Conformal Conic `Proj4.Projection` based on domain definitions in `d`
"""
Plcc(d::Domain) = Proj4.Projection("+proj=lcc +R=$Rearth +lat_0=$(d.LAT0) +lon_0=$(d.LON0) +lat_1=$(d.LAT0) +lat_2=$(d.LAT0) +a=$Rearth +b=$Rearth")

"""
    Plonlat()

Returns longlat `Proj4.Projection`
"""
Plonlat() = Proj4.Projection("+proj=longlat +R=$Rearth")

"""
    lonlat2lcc(d,lonlat)

Returns Lambert Conformal Conic projection coordinates for `lonlat` using domain definition from `d`
"""
lonlat2lcc(d::Domain, lonlat) = Proj4.transform(Plonlat(), Plcc(d),  lonlat)

"""
    lcc2lonlat(d,xy)

Returns lonlat coordinates for Lambert Conformal Conic projection coordinates `xy` using domain definitions from `d`
"""
lcc2lonlat(d::Domain, xy) = Proj4.transform(Plcc(d), Plonlat(),  xy)

"""
    getgridpoints(d; gsize=d.GSIZE)

Returns an array with the lonlat coordinates of the grid points
"""
function getgridpoints(d; gsize = d.GSIZE)
    v = get_lcc_val(d)
    xval = range(v.xl, stop = v.xr, step = gsize)
    yval = range(v.yb, stop = v.yt, step = gsize)
    lcc = Plcc(d)
    lonlat = Plonlat()
    return [Proj4.transform(lcc, lonlat, [x, y]) for x in xval for y in yval]
end

include("in.jl")
include("get_lcc_val.jl")

end # module Domains


text:
## Author: Pearson
## Date: April 2020
## Purpose: Functions for Simulation of EKNR one city case
## Julia version: 1.3.1

function Fun_1Cty_Levels!(
    res_lev::AbstractArray,
    guess_lev::AbstractArray,
    init_lev::NamedTuple,
    exos_lev::NamedTuple,
    params_lev::NamedTuple,
)
    @unpack K₁ = init_lev
    @unpack χ, Aᴰ, L, Kₛₛ = exos_lev
    @unpack βᴷ, βᴸ, δ, ρ, B, T = params_lev
    Y = similar(guess_lev)
    K = similar(guess_lev)
    Y[:] = guess_lev
    K[1] = K₁
    for tt = 2:T
        K[tt] = χ[tt-1]*Aᴰ[tt-1]*B*K[tt-1]*(K[tt-1]/L[tt-1])^(βᴷ-1)*(Y[tt-1]-1)/Y[tt-1]+
            (1-δ)*K[tt-1]
        res_lev[tt-1] = abs(Y[tt]/(Y[tt-1]/(χ[tt-1]*Aᴰ[tt-1]*B*K[tt-1])*(K[tt-1]/L[tt-1])^βᴸ/ρ*
            K[tt]/((1-δ)/(χ[tt]*Aᴰ[tt]*B)*(K[tt]/L[tt])^βᴸ+βᴷ))-1)
    end
    res_lev[T] = abs(K[T]/Kₛₛ-1)
    return res_lev, K, Y
end

function Fun_1Cty_Changes!(
    res_hat::AbstractArray,
    guess_hat::AbstractArray,
    init_hat::NamedTuple,
    exos_hat::NamedTuple,
    params_hat::NamedTuple,
)
    @unpack Y₁ = init_hat
    @unpack χ̂, Âᴰ, L̂ = exos_hat
    @unpack βᴷ, βᴸ, δ, ρ, B, T = params_hat
    K̂=similar(guess_hat)
    Ŷ=similar(guess_hat)
    Y=similar(guess_hat)
    K̂[1] = guess_hat[1]
    Ŷ[1:T-1] = guess_hat[2:T]
    Y[1] = Y₁
    for tt = 2:T
        Y[tt] = Y[tt-1]*Ŷ[tt-1]
    end
    for tt = 2:T
        K̂[tt] = χ̂[tt-1]*Âᴰ[tt-1]*(K̂[tt-1]/L̂[tt-1])^(βᴷ-1)*(Y[tt-1]*Ŷ[tt-1]-1)/
            ((Y[tt-1]-1)*Ŷ[tt-1])*(K̂[tt-1]-(1-δ))+(1-δ)
        res_hat[tt-1] = abs(Ŷ[tt-1]/(K̂[tt-1]/ρ/((1-δ)/(χ̂[tt-1]*Âᴰ[tt-1])*(K̂[tt-1]/L̂[tt-1])^βᴸ+
            βᴷ*(K̂[tt-1]-(1-δ))*(Y[tt-1]/(Y[tt-1]-1))))-1)
    end
    res_hat[T] = abs(K̂[T]/1-1)
    return res_hat, K̂, Ŷ
end


text:
using LinearAlgebra
using LogExpFunctions
using Statistics
using Tullio

"""
    gpd_fit(
        sample::AbstractVector{T<:Real},
        r_eff::T = 1;
        wip::Bool=true,
        min_grid_pts::Integer=30,
        sort_sample::Bool=false
    ) -> (ξ::T, σ::T)

Return a named list of estimates for the parameters ξ (shape) and σ (scale) of the
generalized Pareto distribution (GPD), assuming the location parameter is 0.

# Arguments

  - `sample::AbstractVector`: A numeric vector. The sample from which to estimate
    the parameters.
  - `wip::Bool = true`: Logical indicating whether to adjust ξ based on a weakly informative
    Gaussian prior centered on 0.5. Defaults to `true`.
  - `min_grid_pts::Integer = 30`: The minimum number of grid points used in the fitting
    algorithm. The actual number used is `min_grid_pts + ⌊sqrt(length(sample))⌋`.

# Note

Estimation method taken from Zhang, J. and Stephens, M.A. (2009). The parameter ξ is the
negative of k.
"""
function gpd_fit(
    sample::AbstractVector{T},
    r_eff::T=1;
    wip::Bool=true,
    min_grid_pts::Integer=30,
    sort_sample::Bool=false,
) where T<:Real

    len = length(sample)
    # sample must be sorted, but we can skip if sample is already sorted
    if sort_sample
        sample = sort(sample; alg=QuickSort)
    end

    grid_size = min_grid_pts + isqrt(len)  # isqrt = floor sqrt
    n_0 = 10  # determines how strongly to nudge ξ towards .5
    x_star = inv(3 * sample[(len + 2) ÷ 4])  # magic number. ¯\_(ツ)_/¯
    invmax = inv(sample[len])

    # build pointwise estimates of ξ and θ at each grid point
    θ_hats = similar(sample, grid_size)
    @fastmath @. θ_hats = invmax + (1 - sqrt((grid_size + 1) / $(1:grid_size))) * x_star
    @tullio threads=false ξ_hats[i] := log1p(-θ_hats[i] * sample[j])
    ξ_hats /= len

    log_like = ξ_hats  # Reuse preallocated array
    # Calculate profile log-likelihood at each estimate:
    @tullio threads=false ξ_hats[i] =
        len * (log(-θ_hats[i] / ξ_hats[i]) - ξ_hats[i] - 1)
    # Calculate weights from log-likelihood:
    weights = log_like  # Reuse preallocated array
    log_norm = logsumexp(log_like)
    @tullio threads=false log_like[x] = exp(log_like[x] - log_norm)
    # Take weighted mean:
    @tullio threads=false θ_hat := weights[x] * θ_hats[x]
    @tullio threads=false ξ := log1p(-θ_hat * sample[i])
    ξ /= len
    σ::T = -ξ / θ_hat

    # Drag towards .5 to reduce variance for small len
    if wip
        @fastmath ξ = (r_eff * ξ * len + 0.5 * n_0) / (r_eff * len + n_0)
    end

    return ξ, σ

end

"""
    gpd_quantile(p::T, k::T, sigma::T) where {T<:Real} -> T

Compute the `p` quantile of the Generalized Pareto Distribution (GPD).

# Arguments

  - `p`: A scalar between 0 and 1.
  - `ξ`: A scalar shape parameter.
  - `σ`: A scalar scale parameter.

# Returns

A quantile of the Generalized Pareto Distribution.
"""
function gpd_quantile(p, ξ::T, sigma::T) where {T <: Real}
    return sigma * expm1(-ξ * log1p(-p)) / ξ
end

max_stars_repo_name: ParadaCarleton/JuLOOa

max_issues_repo_name: ParadaCarleton/JuLOOa

max_forks_repo_name: ParadaCarleton/JuLOOa


text:
module sNMF_MC

include("../Misc.jl")
include("BackwardKernel.jl")
using .Misc, .BackwardKernel

using Distributions, SpecialFunctions

export Particle, smc_weight, particle_mcmc

mutable struct Particle{ℜ <: Real}
    S_R::Array{ℜ,1}
    S_IR::Array{ℜ,2}
    function Particle(I::Int, R::Int; γ::ℜ=0.1) where {ℜ<:Real}
        a::ℜ = I*I*γ
        return new{ℜ}(fill(a/R,R), fill(2.0 * a/(I*R),I,R))
    end
    function Particle(S_R::Array{ℜ,1}, S_IR::Array{ℜ,2}) where {ℜ<:Real}
        return new{ℜ}(S_R, S_IR)
    end
end

function smc_weight(X::Array{ℜ,2}, R::Ƶ, N::Ƶ=1; μ::ℜ=1.0, γ::ℜ=0.1, adaptive::Bool=true) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ = size(X,1)
    T::Ƶ = Ƶ(sum(X))

    ESS::Array{ℜ} = zeros(ℜ,T)
    ESS_min::Array{ℜ} = adaptive ? fill(N/2.0,T) : fill(N+1.0,T)

    P = [Particle(I,R; γ=γ) for n=1:N]
    P_temp = [Particle(I,R; γ=γ) for n=1:N]

    a::ℜ, b::ℜ = I*I*γ, γ/μ
    log_Zᵗ::ℜ = a*log(b) - (a+T)*log(b + 1.0) - sum(lgamma.(X .+ 1.0))

    log_wᵗ::Array{ℜ}, log_Wᵗ::Array{ℜ} = fill(log_Zᵗ,N), fill(-log(N),N)
    Wᵗ::Array{ℜ}, cum_Wᵗ::Array{ℜ} = fill(1.0/N,N), zeros(ℜ,N)

    log_νᵗ::ℜ, log_q_r::Array{ℜ}, q_r::Array{ℜ} = 0.0, zeros(ℜ,R), zeros(ℜ,R)
    rᵗ::Ƶ, i₁ᵗ::Ƶ, i₂ᵗ::Ƶ, uᵗ::ℜ = 0, 0, 0, 0.0

    for (t,(i₁ᵗ, i₂ᵗ)) in enumerate(EventQueue(X))
        for (n,p) in enumerate(P)

            log_q_r .= log.(p.S_IR[i₁ᵗ,:]) .+ log.(p.S_IR[i₂ᵗ,:] .+ (i₁ᵗ == i₂ᵗ)) .- log.(4.0 .* p.S_R .+ 2)
            log_νᵗ = logsumexp(log_q_r)
            log_q_r .-= log_νᵗ

            q_r .= exp.(log_q_r)
            rᵗ = rand(Categorical(q_r))

            p.S_R[rᵗ] += 1.0
            p.S_IR[i₁ᵗ,rᵗ] += 1.0
            p.S_IR[i₂ᵗ,rᵗ] += 1.0

            log_wᵗ[n] += log_νᵗ
        end

        log_Zᵗ = logsumexp(log_wᵗ)
        log_Wᵗ .= log_wᵗ .- log_Zᵗ
        log_Zᵗ -= log(N)
        Wᵗ .= exp.(log_Wᵗ)

        ESS[t] = 1.0/sum(Wᵗ .* Wᵗ)

        if ESS[t] < ESS_min[t]
            for n=1:N
                P_temp[n].S_R .= P[n].S_R
                P_temp[n].S_IR .= P[n].S_IR
            end

            cum_Wᵗ .= cumsum(Wᵗ)
            uᵗ = rand()/N
            p_id = 1

            for n=1:N # systematic resampling
                while cum_Wᵗ[p_id] < uᵗ
                    p_id += 1
                end
                P[n].S_R .= P_temp[p_id].S_R
                P[n].S_IR .= P_temp[p_id].S_IR

                uᵗ += 1.0/N
            end
            log_wᵗ .= log_Zᵗ
        end
    end
    p_id = rand(Categorical(Wᵗ))
    P[p_id].S_R .-= a/R
    P[p_id].S_IR .-= 2.0 .* a/I*R
    return log_Zᵗ, ESS, P[p_id]
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
# Taken from https://github.com/JuliaData/DataFrames.jl/blob/main/test/runtests.jl
using NighttimeLights
using Test

fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
quiet = length(ARGS) > 0 && ARGS[1] == "-q"
anyerrors = false

my_tests = ["aggregate.jl", "coordinate_system.jl", "data_io.jl", "f_apply.jl", "mask_area.jl", "polygons.jl", "sparse_cube.jl",

"data_cleaning/background_noise_removal.jl", "data_cleaning/bias_correction.jl", "data_cleaning/interpolation.jl", "data_cleaning/mark_nan.jl", "data_cleaning/outlier_removal.jl", "data_cleaning/full_procedures.jl",

"other/detrend.jl", "other/nan_functions.jl", "other/rank_correlation.jl", "other/weighted_mean.jl"]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch e
        global anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow(e)
        elseif !quiet
            showerror(stdout, e, backtrace())
            println()
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Documenter, ADCME
makedocs(sitename="ADCME", modules=[ADCME],
pages = Any[
    "index.md",
    "Tutorial"=>["tutorial.md","videos_and_slides.md", "tu_whatis.md", "tu_basic.md", "tu_optimization.md", "tu_sparse.md", "tu_fd.md", "tu_fem.md",
        "tu_inv.md", "tu_recipe.md", "tu_nn.md", "tu_implicit.md", "tu_customop.md","tu_debug.md", "exercise.md"],
    "Resources" => ["newton_raphson.md", "parallel.md", "optimizers.md", "optim.md", "ode.md", "global.md",
            "julia_customop.md", "nn.md", "ot.md", "resource_manager.md", "alphascheme.md", "factorization.md", "customopt.md",
            "options.md", "mcmc.md", "mpi.md", "mpi_benchmark.md", "multithreading.md", "rbf.md", "topopt.md", "quadrature.md",
            "sqlite3.md", "hessian.md", "plotly.md"],
    "Physics Informed Machine Learning" => ["fdtd.md"],
    "Deep Learning Schemes" => ["vae.md", "flow.md", "convnet.md", "bnn.md", "reinforcement_learning.md"],
    "Developer Guide" => ["designpattern.md", "toolchain.md", "installmpi.md", "windows_installation.md"],
    "Applications" => ["apps.md", "apps_ana.md", "apps_levy.md",
            "apps_constitutive_law.md", "apps_ad.md", "apps_adseismic.md", "apps_nnfem.md"],
    "api.md"
],
authors = "Kailai Xu")

deploydocs(
    repo = "github.com/kailaix/ADCME.jl.git",
)


text:
# An script to facilitate reading input concentrations,
# creating phreeqc file, running it, reding the output, and plotting it

# load the required packages
using CSV, DataFrames, LsqFit
# using LsqFit, Optim
include("SCM_functions.jl")

# read the data file
input_data = CSV.read("sample_input.csv")
capacitance = [1.3, 4.5]
zeta_exp = input_data[:zeta]

# calculate zeta potential for all data

function zeta_calc(chalk_eq, input_data, capacitance)
    (N, C) = size(input_data)
    MW_calcite = 100.09 # g/mol
    phqc_solution = Array{Any}(N)
    for i in 1:N
        phqc_solution[i] = """
        SOLUTION $i
            temp  $(input_data[:temp][i])
            pressure  $(input_data[:pressure][i])
            -units mol/kgw
            pH     $(input_data[:pH][i])   charge
            Ba     $(input_data[:Ba][i])
            Ca     $(input_data[:Ca][i])
            Mg     $(input_data[:Mg][i])
            Na     $(input_data[:Na][i])
            K      $(input_data[:K][i])
            Cl     $(input_data[:Cl][i])
            S      $(input_data[:S][i])
            C      $(input_data[:C][i]) as HCO3
        EQUILIBRIUM_PHASES $i
            CO2(g)   $(-input_data[:p_CO2][i])  1.0
            Calcite  0.0  $(input_data[:calcite][i]/MW_calcite)
        # END
        SURFACE $i
            -sites_units density
            -cd_music
            -equilibrate $i
            Chalk_a      $(input_data[:site_dens][i])      $(input_data[:sp_area][i])    $(input_data[:calcite][i])
            Chalk_c      $(input_data[:site_dens][i])
            -capacitance $(capacitance[1]) $(capacitance[2])
            #-donnan # uncomment, and see the effect on Charge Balance...
        """
    end

    # creae a selected output block
    sel_output = """
    SELECTED_OUTPUT 1
        -file phqc_out.csv
        -reset false
        USER_PUNCH
                -headings dielec sigma_chalk psi_chalk sigma_clay psi_clay psi_chalk1 psi_chalk2
                10 PUNCH EPS_R
                30 PUNCH EDL("sigma", "Chalk")
                50 PUNCH EDL("psi", "Chalk")
                51 PUNCH EDL("sigma", "Clay")
                52 PUNCH EDL("psi", "Clay")
                100 PUNCH EDL("psi1", "Chalk")
                110 PUNCH EDL("psi2", "Chalk")
        END
    END
    """
    # create the surface reactions
    surface_reactions=modified_so(log_k=chalk_eq)
    # write everything into a file
    zeta_pot = zeros(N)
    for i in 1:N
        ph_file=open("phreeqc_file", "w")
        write(ph_file, surface_reactions*phqc_solution[i]*sel_output)
        close(ph_file)
        run(`phreeqc phreeqc_file`)
        ph_res = CSV.read("phqc_out.csv"; delim = "\t  ")
        zeta_pot[i] = ph_res[:psi_chalk2][end]
    end
    return zeta_pot
end # end of zeta_calc function

complex_chalk_brine = param -> zeta_calc(param, input_data, capacitance)
model = (x, param) -> zeta_calc(param, input_data, capacitance)

function objective_function(param::Vector{Float64}, grad::Vector{Float64})
    zeta_calc = complex_chalk_brine(param)
    zeta_error=(zeta_calc-zeta_exp)
    if length(grad)>0
      dx=1e-5
      for i in 1:length(param)
        x=copy(param)
        x[i]+=dx
        grad_zeta=(complex_chalk_brine(x)-zeta_calc)/dx
        grad[i]=2.0*sum(grad_zeta.*zeta_error)
      end
    end
    return sum(zeta_error.*zeta_error)
end

function myconstraint(results::Vector{Float64}, param::Vector{Float64}, grad::Array{Float64,2})
    # zeta_exp.*zeta_calc > 0.0
    # -zeta_exp*zeta_calc<0
    # grad nxm [j,i] dci/dxj
    zeta_calc = complex_chalk_brine(param)
    zeta_mult = zeta_calc.*zeta_exp
    m=length(zeta_exp)
    n=length(param)
    if length(grad)>0
      dx=1e-5
      for i in 1:n
        x=copy(param)
        x[i]+=dx
        grad_zeta=(complex_chalk_brine(x)-zeta_calc)/dx
        grad[i,:]=-zeta_exp.*grad_zeta
      end
    end
    results[:]=-zeta_mult
end

  # LN_COBYLA works fine
  # LD_SLSQP
chalk_logk_init =[-3.58, -2.8, -2.2, 12.85, -24.73, 10.15, 1.55, 0.35, -9.8144, -7.5]

chalk_logk_lb=chalk_logk_init-0.5*abs(chalk_logk_init)
chalk_logk_ub=chalk_logk_init+0.5*abs(chalk_logk_init)
N_data = length(zeta_exp)
w = ones(N_data)
x_dummy = collect(linspace(1,length(zeta_exp)))
fit = curve_fit(model, x_dummy, zeta_exp, w, chalk_logk_init, lower = chalk_logk_lb, upper = chalk_logk_ub)

println(fit.param)

max_stars_repo_name: simulkade/calciteSCM

max_issues_repo_name: simulkade/calciteSCM

max_forks_repo_name: simulkade/calciteSCM


text:
function calc_wagons_time_from_first_to_last(s::Simulation, know_which_route::Bool = false)
	"""
	Description:
		For each TTC route calculates how many seconds it takes to get from the first station to the last on the route.
	Returns:
		Array with unique values or Dictionary with value per route.
	"""
	if know_which_route == true
		wagons_time_from_first_to_last = Dict{String,Float64}()
		for wag in s.wagons
			if wag.line ∉ keys(wagons_time_from_first_to_last)
				wagons_time_from_first_to_last[wag.line] = sum(values(wag.path_distances))
			end
		end
	elseif know_which_route == false
		wagons_time_from_first_to_last = Float64[]
		for wag in s.wagons
		   push!(wagons_time_from_first_to_last,sum(values(wag.path_distances)))
		end
	end

    return unique(wagons_time_from_first_to_last)
end

function calc_agents_time_to_work(s::Simulation)
	"""
	Description:
		For each agent calculates how many seconds it takes to get from home to work.
	Returns:
		Array with values per agent.
	"""
    agents_time_to_work = Float64[]
    for ag in s.agents
		sums = [sum(values(value)) for (key, value) in ag.distances_home_work]
		push!(agents_time_to_work, mean(sums))
    end
    return agents_time_to_work
end

function TTC_usage_by_agents(s::Simulation, orig_map_nodes_num::Int64)
	"""
	Description:
		For each agent calculates:
			1) how many seconds he/she spends in TTC on the way from home to work
			2) how many different TTC types he/she uses to get from home to work
			3) percent of time spent in TTC in relation to the total time from home to work
	Returns:
		time_spend_in_TTC, TTC_types_used, time_spend_in_TTC_comp_to_total_path_time
	"""
    time_spend_in_TTC = Float64[]
    TTC_types_used = Int64[]
    agents_time_to_work = calc_agents_time_to_work(s);

    for ag in s.agents
        ag_TTC_types_used = Int64[]
        ag_time_in_TTC = Dict()
		TTC_on_path = Dict()

		for (key, value) in ag.paths_home_work
			TTC_on_path[key] = [i for i in keys(value) if i > orig_map_nodes_num]
		end

		for (key, value) in TTC_on_path
			ag_time_in_TTC[key] = 0
			for ttc_node in value
            	ag_time_in_TTC[key] += ag.distances_home_work[key][ttc_node]
            	push!(ag_TTC_types_used, route_types[ttc_node-orig_map_nodes_num])
			end
        end
		ag_time_in_TTC = mean(values(ag_time_in_TTC))

        push!(time_spend_in_TTC, ag_time_in_TTC)
        push!(TTC_types_used, length(unique(ag_TTC_types_used)))
    end
    time_spend_in_TTC_comp_to_total_path_time  = time_spend_in_TTC ./ agents_time_to_work
    return time_spend_in_TTC, TTC_types_used, time_spend_in_TTC_comp_to_total_path_time
end

println("Gather statistics functions were read")

### Examples ###
# maximum(values(s.max_passengers_per_TTC_car))

# agents_time_to_work = calc_agents_time_to_work(s);
# describe(agents_time_to_work ./ 60)

# wagons_time_from_first_to_last = calc_wagons_time_from_first_to_last(s);
# describe(wagons_time_from_first_to_last ./ 60)

# time_spend_in_TTC, TTC_types_used, time_spend_in_TTC_comp_to_total_path_time = TTC_usage_by_agents(s,orig_map_nodes_num);
# describe(time_spend_in_TTC ./ 60)
# describe(time_spend_in_TTC_comp_to_total_path_time)
# describe(TTC_types_used)

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
import Gurobi
const GRBENV = Gurobi.Env()

function run_rmp_gurobi(finst)
    grb = Gurobi.Model(GRBENV, "")
    Gurobi.read_model(grb, finst)

    # Set parameters
    Gurobi.set_int_param!(grb, "OutputFlag", 1)
    Gurobi.set_int_param!(grb, "Threads", 1)
    Gurobi.set_int_param!(grb, "Method", 2)
    Gurobi.set_int_param!(grb, "Crossover", 0)
    Gurobi.set_int_param!(grb, "Presolve", 0)

    # Solve model
    t = @elapsed Gurobi.optimize(grb)

    println("Solve time: $t")
    println("Barrier iterations: ", Gurobi.get_intattr(grb, "BarIterCount"))
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_rmp_gurobi(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
using LeetCode
using Test