text:
macros = Dict()

defs = Dict{String, Union{JLPAbstValue, Vector}}()
#("Int" => JLPPrimitive("Int", 32),"Float" => JLPFloats("Float", 32))

asmcallfuncs = JLPAbstAsmFunc[]

function resetgv()
    macros = Dict()
    defs = Dict{String, Union{JLPAbstValue, Vector}}("Int" => JLPPrimitive("Int", 32)
                ,"Float" => JLPFloats("Float", 32))
    asmcallfuncs = JLPAbstAsmFunc[]
end

function _expandmacro(replacelists, l, newlist=JLPNil())
    if !(typeof(l) <: JLPList)
        return l
    else
        cdr = _expandmacro(replacelists, l.cdr, newlist)
        car = _expandmacro(replacelists, l.car, newlist)
        val = l.val
        if haskey(replacelists, l.car.val)
            if l.car.val[1:2] == "\$&"
                car = replacelists[l.car.val].car
                cdr = add_list(replacelists[l.car.val].cdr, cdr)
                val = replacelists[l.car.val].val
            else
                car = replacelists[l.car.val]
                val = JLPNil()
            end
        end
        return JLPList(val, car, cdr)
    end
end

function expandmacro(jlpmacroname, l)
    jlpmacro = macros[jlpmacroname]
    margs = jlpmacro.args
    sargs = Dict()
    m = margs.car
    tmplist = l
    while !(typeof(m) <: JLPNil)
        if m.car.val[1] == '&'
            sargs["\$" * m.car.val] = tmplist
            break
        end
        sargs["\$" * m.car.val] = tmplist.car
        m = m.cdr
        tmplist = tmplist.cdr
    end
    return _expandmacro(sargs, jlpmacro.returnlist)
end

function expandallmacro(lst)
    if !(typeof(lst) <: JLPList)
        return lst
    end
    if lst.car == jlpnil == lst.cdr
        return lst
    end

    newlst = if lst.cdr == jlpnil
                 jlpnil
             else
                 JLPList(lst.cdr.val, expandallmacro(lst.cdr.car), expandallmacro(lst.cdr.cdr))
             end
    newlst = expandallmacro(newlst)
    if haskey(macros, lst.car.val)
        return expandmacro(lst.car.val, newlst).car
    else
        return JLPList(lst.val, lst.car, newlst)
    end
end

function setmacro(l)
    macroname = l.car.val
    args = l.cdr
    returnlist = l.cdr.cdr
    macros[macroname] = JLPMacro(macroname, args, returnlist)
end

function getreturnfunc(l)
    lc = l.car
    tlc = typeof(l.car)
    returntype = if tlc <: JLPSymbol
                     getreturntype(lc)
                 else
                     JLPNil()
                 end
    args = JLPAbstValue[]
    arg = l.cdr.car
    while !(typeof(arg) <: JLPNil) && !(typeof(arg.cdr) <: JLPNil)
        tmp = getreturntype(arg.cdr.car)
        if typeof(tmp) <: JLPGenericStruct && tmp.name == "Vector"
            n = defs[arg.cdr.car.cdr.car.val]
            ln = "Vector___$(n.name)___"
            defs[ln] = JLPVectorType(ln, n)
            generatevectordefs(defs[ln])
            tmp = defs[ln]
        end
        push!(args, tmp)
        arg = arg.cdr.cdr
    end
    return JLPFn(returntype, args)
end

function getreturntype(l)
    tl = typeof(l)
    if tl <: JLPNil
        return JLPNil()
    elseif tl <: JLPSymbol
        if !haskey(defs, l.val)
            setdef(findfirst(x -> x.car.val == l.car.val, deflist))
        end
        @assert haskey(defs, l.val)
        return defs[l.val]
    elseif tl <: JLPInteger
        return defs["Int"]
    elseif tl <: JLPFloat
        return defs["Float"]
    elseif tl <: JLPBool
        return defs["Bool"]
    end
    if tl <: JLPList
        lc = l.car
        tlc = typeof(lc)
        if !(tlc <: JLPList)
            return getreturntype(lc)
        end
        lcc = lc.car
        if lcc.val == "fn"
            return getreturnfunc(lc.cdr)
        elseif lcc.val == "let"
            lets::JLPAbstList = lc.cdr
            arg = lets.car
            befdefs = Dict()
            letsyms = []
            while !(typeof(arg) <: JLPNil)
                @assert arg <: JLPSymbol "arg must be JLPSymbol"
                push!(letsyms, arg.val)
                if haskey(defs, arg.val)
                    befdefs[arg.val] = defs[arg.val]
                end
                setdef(lets)
                lets = lets.cdr.cdr
                arg = lets.cdr.cdr.car
            end
            map(x -> delete!(defs, x), letsyms)
            merge!(defs, befdefs)
            evall = l.cdr.cdr
            if !(typeof(evall) <: JLPNil)
                return JLPNil()
            end
            while !(typeof(evall.cdr) <: JLPNil)
                evall = evall.cdr
            end
            return getreturntype(evall)
        elseif lcc.val == "ret"
            return getreturntype(lc.cdr.car)
        elseif lcc.val == "if"
            ifcond::JLPAbstList = lc.cdr.car
            ift = getreturntype(lc.cdr.cdr.car)
            iff = getreturntype(lc.cdr.cdr.cdr.car)
            @assert typeof(ifcond) <: defs["Bool"]
            @assert ift == iff "!! true result != false result !!"
            return ift
        else
            funcs = if typeof(lc) <: JLPAbstList
                       getreturntype(lc)
                   else
                       defs[lc.val]
                   end
            if typeof(funcs) <: JLPGenericStruct
                name = funcs.name
                realtp = []
                ln = "$(name)___"
                for lst in l.cdr
                    tmp = getreturntype(lst)
                    push!(realtp, tmp)
                    ln *= tmp.name
                    ln *= "___"
                end
                tmp = get(defs, ln, nothing)
                if !isnothing(tmp)
                    return tmp
                end

                if name == "Vector"
                    defs[ln] = JLPVectorType(ln, realtp[1])
                    generatevectordefs(defs[ln])
                    return defs[ln]
                else
                    ftypeparams = funcs.typeparams
                    fparams = funcs.params
                    params = []
                    for i in 1:length(typeparams)
                        if typeparams[i].name
                        end
                    end
                end
            else
                args = JLPAbstValue[]
                arg = l.cdr
                while !(typeof(arg) <: JLPNil)
                    push!(args, getreturntype(arg.car))
                    arg = arg.cdr
                end
                f = findfirst(x -> x.args == args, funcs)
                return f.returntype
            end
        end
    end
end

function setdef(l)
    defname = l.car.val
    returntype = getreturntype(l.cdr)
    trt = typeof(returntype)
    if trt <: JLPFn
        if haskey(defs, defname)
            push!(defs[defname], returntype)
        else
            defs[defname] = [returntype]
        end
        asmname = returntype.returntype.name * "_" * defname * foldl((a,b) -> a * "_" * b.name, returntype.args, init="")
        push!(asmcallfuncs, JLPAsmFunc(defname, asmname, returntype))
    elseif trt <: JLPPrimitive
        defs[defname] = returntype
    elseif trt <: JLPStruct
        defs[defname] = returntype
    end
end

function generateprimnew(prim)
    pt = prim.second
end

function generateprimop(prim, name, asmname, option, rept = prim.second)
    pt = prim.second
    tmp = JLPFn(rept, [pt,pt])
    if haskey(defs, name)
        push!(defs[name], tmp)
    else
        defs[name] = [tmp]
    end
    push!(asmcallfuncs, JLPBinaryOp(name, tmp, asmname, option))
end

generateprimadd(prim) = generateprimop(prim, "+", "add", "")

generateprimsub(prim) = generateprimop(prim, "-", "sub", "")

generateprimmul(prim) = generateprimop(prim, "*", "mul", "")

generateprimdiv(prim) = generateprimop(prim, "/", "sdiv", "")

generateprimeq(prim) = generateprimop(prim, "=", "icmp", " eq", defs["Bool"])

generateprimne(prim) = generateprimop(prim, "/=", "icmp", " nq", defs["Bool"])

function generateprimfuncs(defprims)
    for prim in defprims
        generateprimnew(prim)
        generateprimadd(prim)
        generateprimsub(prim)
        generateprimmul(prim)
        generateprimdiv(prim)
        generateprimeq(prim)
        generateprimne(prim)
    end
end

function generatevectordefs(ve)
    returntype = ve.T
    asmname = ve.T.name
    setf = JLPFn(ve, [ve , defs["Int"], ve.T])
    setname = "set_$(ve.name)"
    getf = JLPFn(ve.T, [ve, defs["Int"]])
    getname = "get_$(ve.name)"
    if haskey(defs, "set")
        push!(defs["set"], setf)
    else
        defs["set"] = [setf]
    end
    if haskey(defs, "get")
        push!(defs["get"], getf)
    else
        defs["get"] = [getf]
    end
    push!(asmcallfuncs, JLPAsmStore("set", setf, setname))
    push!(asmcallfuncs, JLPAsmLoad("get", getf, getname))
    #push!(asmcallfuncs, )
end

function setprimitives(ls)
    defprims = Dict()
    for l in ls
        defprims[l.car.val] = JLPPrimitive(l.car.val, l.cdr.car.val)
    end
    merge!(defs, defprims)
    generateprimfuncs(defprims)
end

function generateasmcode(l::JLPSexp, asms::Vector, lets::Dict, count::Count = Count(1,1))
    if typeof(l) <: JLPList
        if typeof(l.car) <: JLPSymbol && l.car.val == "if"
            condition = generateasmcode(l.cdr.car, asms, lets, count)
            tl = JLPAsmLabel("tlabel_" * string(count.labelc))
            fl = JLPAsmLabel("flabel_" * string(count.labelc))
            el = JLPAsmLabel("elabel_" * string(count.labelc))
            count.labelc += 1
            sc = count.labelc
            push!(asms, JLPAsmBr(condition, JLPAsmLabel[tl, fl]))
            push!(asms, tl)
            truecode = generateasmcode(l.cdr.cdr.car, asms, lets, count)
            tec = count.labelc
            push!(asms, JLPAsmBr(el))
            push!(asms, fl)
            falsecode = generateasmcode(l.cdr.cdr.cdr.car, asms, lets, count)
            @assert truecode.ty == falsecode.ty "if method must return one type"
            fec = count.labelc
            push!(asms, JLPAsmBr(el))
            push!(asms, el)
            befl = if (tec == sc)
                       tl
                   else
                       JLPAsmLabel("elabel_" * string(sc))
                   end
            tvarlabel = (truecode, befl)
            aftl = if (fec == tec)
                       fl
                   else
                        JLPAsmLabel("elabel_" * string(tec))
                    end
            fvarlabel = (falsecode, aftl)
            varname = JLPAsmVar(count.c, truecode.ty)
            count.c += 1
            push!(asms, JLPAssignment(varname, JLPAsmPhi(truecode.ty, [tvarlabel, fvarlabel])))
            return varname
        elseif typeof(l.car) <: JLPSymbol && l.car.val == "ret"
            ret = JLPAsmRet(generateasmcode(l.cdr.car, asms, lets, count))
            push!(asms, ret)
            return ret
        elseif typeof(l.car) <: JLPSymbol && l.car.val == "let"
            lst = l.cdr.car
            while typeof(lst) <: JLPList && !(typeof(lst.car) <: JLPNil)
                varname = lst.car.val
                varval = generateasmcode(lst.cdr.car, asms, lets, count)
                lets[varname] = varval
                lst = lst.cdr.cdr
            end
            codelst = l.cdr.cdr
            target = generateasmcode(codelst.car, asms, lets, count)
            for lst in codelst.cdr
                target = generateasmcode(lst, asms, lets, count)
            end
            return target
        elseif typeof(l.car) <: JLPList && typeof(l.car.car) <: JLPSymbol && l.car.car.val == "Vector"
            #defs[l.car.cdr.car.val] = lets[l.car.cdr.car.val].ty
            t = getreturntype(l)
            int = generateasmcode(l.cdr.car, asms, lets, count)
            tmp = JLPVector(t.name, t.T, int)
            varname = JLPAsmVar(count.c, defs[tmp.name])
            count.c += 1
            lets[varname.name] = varname
            push!(asms, JLPAssignment(varname, tmp))
            return varname
        else
            cs = []
            tys = []
            for lst in l.cdr
                push!(cs, generateasmcode(lst, asms, lets, count))
                tmp = get(lets, cs[end].name, nothing)
                if !isnothing(tmp)
                    push!(tys, tmp)
                else
                    push!(tys, get(defs, cs[end].name, nothing))
                end
            end
            asmfunc = asmcallfuncs[findfirst(x -> x.name == l.car.val && x.fn.args == map(x -> x.ty, tys), asmcallfuncs)]
            varname = JLPAsmVar(count.c, asmfunc.fn.returntype)
            count.c += 1
            lets[varname.name] = varname
            push!(asms, JLPAssignment(varname, JLPAsmCallFunc(asmfunc, tys)))
            return varname
        end
    else
        if haskey(lets, l.val)
            return lets[l.val]
        elseif haskey(defs, l.val)
            return defs[l.val]
        elseif typeof(l) <: JLPInteger
            varname = JLPAsmVar(count.c, defs["Int"])
            count.c += 1
            push!(asms, JLPAssignment(varname, JLPPrimNew(defs["Int"], l.val)))
            lets[varname.name] = varname
            return varname
        else
            @assert false "$(l.val) is undefined"
        end
    end
end

function generateasmdef(l::JLPList)
    name::String = l.car.val
    lda = l.cdr.car
    if typeof(lda) <: JLPList
        ldaa = lda.car
        if ldaa.val == "fn"
            ldad = lda.cdr
            returntype = defs[ldad.car.val]
            ldadd = ldad.cdr
            arglst = ldadd.car
            args = Dict()
            argl = arglst.car
            argv = []
            argty = []
            while argl != jlpnil
                if typeof(arglst.cdr.car) <: JLPList
                    tmp = getreturntype(arglst.cdr.car)
                else
                    tmp = defs[arglst.cdr.car.val]
                end
                if typeof(tmp) <: JLPGenericStruct && tmp.name == "Vector"
                    n = defs[arglst.cdr.car.cdr.car.val]
                    ln = "Vector___$(n.name)___"
                    defs[ln] = JLPVectorType(ln, n)
                    generatevectordefs(defs[ln])
                    tmp = defs[ln]
                end
                push!(argty, tmp)
                push!(argv, JLPAsmVar(argl.val, tmp))
                args[argl.val] = JLPAsmVar(argl.val, tmp)
                arglst = arglst.cdr.cdr
                if arglst == jlpnil
                    break
                end
                argl = arglst.car
            end
            asms = []
            count = Count(1,1)
            rettarget = ""
            for lst in ldadd.cdr
                rettarget = generateasmcode(lst, asms, args, count)
            end
            if typeof(rettarget) <: JLPAsmVar
                push!(asms, JLPAsmRet(rettarget))
            elseif typeof(asms[end]) <: JLPAssignment
                target = asms[end].target
                push!(asms, JLPAsmRet(target))
            else
                @assert false "end must be assinment struct"
            end
            asmfunc = asmcallfuncs[findfirst(x -> x.name == l.car.val && x.fn.args == argty, asmcallfuncs)]
            return JLPLLVMFunc(asmfunc.asmname, returntype, argv, asms)
        end
    else
        return []
    end
end

emit(target::StringBuilder, str::String) = append!(target, str * "\n")

addcode(target::StringBuilder, str::String) = append!(target, "  $str\n")

function strllvmtype(ty::JLPAbstValue)
    tty = typeof(ty)
    if tty <: JLPPrimitive
        return "i" * string(ty.bitlen) * ty.option
    else tty <: JLPVectorType
        return (strllvmtype(ty.T) * "*")
    end
end

function generatellvmirfn(f)
    sb = StringBuilder()
    name = f.name
    returntype = f.returntype
    args = f.args
    append!(sb, "define ")
    if typeof(returntype) <: JLPPrimitive
        append!(sb, "i$(string(returntype.bitlen)) ")
    end
    append!(sb, "@$name(")
    for i in 1:length(args)
        append!(sb, "$(strllvmtype(args[i].ty)) ")
        append!(sb, "%$(args[i].name)")
        if i != length(args)
            append!(sb, ", ")
        end
    end
    append!(sb, ") {\n")
    asms = f.asms
    for asm in asms
        if typeof(asm) <: JLPAssignment
            varname = asm.target.name
            func = ""
            if typeof(asm.val) <: JLPAsmCallFunc
                if typeof(asm.val.fn) <: JLPBinaryOp
                    funcname = asm.val.fn.boname
                    returntype = asm.val.fn.fn.returntype
                    argtype = asm.val.fn.fn.args[1]
                    args = asm.val.args
                    op = asm.val.fn.op
                    argstr = ""
                    for i in 1:length(args)
                        argstr *= "%" * args[i].name
                        if i == length(args)
                            break
                        end
                        argstr *= ", "
                    end
                    func = "$funcname$(op) $(strllvmtype(argtype)) $argstr"
                elseif typeof(asm.val.fn) <: JLPAsmStore
                    returntype = asm.val.fn.fn.returntype
                    args = asm.val.args
                    argtys = asm.val.fn.fn.args
                    argstr = ""
                    tmpstr = strllvmtype(asm.val.fn.fn.returntype)
                    tmpTstr = strllvmtype(asm.val.fn.fn.returntype.T)
                    indexty = strllvmtype(asm.val.args[2].ty)
                    indexname = asm.val.args[2].name

                    addcode(sb, "%p$varname = getelementptr $(tmpTstr), $(tmpstr) %$(asm.val.args[1].name), $indexty %$indexname")
                    addcode(sb, "store $(strllvmtype(args[3].ty)) %$(args[3].name), $(strllvmtype(args[1].ty)) %p$(varname)")
                    func = "getelementptr $tmpTstr, $(tmpstr) %$(asm.val.args[1].name), $indexty 0"
                elseif typeof(asm.val.fn) <: JLPAsmLoad
                    returntype = asm.val.fn.fn.returntype
                    args = asm.val.args
                    argtys = asm.val.fn.fn.args
                    argstr = ""
                    tmpstr = strllvmtype(asm.val.fn.fn.returntype)
                    indexty = strllvmtype(asm.val.args[2].ty)
                    indexname = asm.val.args[2].name

                    addcode(sb, "%p$varname = getelementptr $tmpstr, $(tmpstr)* %$(asm.val.args[1].name), $indexty %$indexname")
                    func = "load $(tmpstr), $(tmpstr)* %p$varname"
                elseif typeof(asm.val.fn) <: JLPAsmFunc
                    funcname = asm.val.fn.asmname
                    returntype = asm.val.fn.fn.returntype
                    args = asm.val.args
                    argtys = asm.val.fn.fn.args
                    argstr = ""
                    for i in 1:length(args)
                        argstr *= strllvmtype(argtys[i]) * " %" * args[i].name
                        if i == length(args)
                            break;
                        end
                        argstr *= ", "
                    end
                    func = "call $(strllvmtype(returntype)) @$funcname($argstr)"
                end
            elseif typeof(asm.val) <: JLPAsmPhi
                func = "phi $(strllvmtype(asm.val.returntype)) "
                for i in 1:length(asm.val.vallabel)
                    (val, label) = asm.val.vallabel[i]
                    func *=  "[%$(val.name) , %$(label.label)]"
                    if i == length(asm.val.vallabel)
                        break
                    end
                    func *= " , "
                end
            elseif typeof(asm.val) <: JLPPrimNew
                returntype = asm.val.type
                func = "add i$(string(returntype.bitlen)) $(asm.val.val), 0"
            elseif typeof(asm.val) <: JLPVector
                len = asm.val.length
                addcode(sb, "%i$varname = sext $(strllvmtype(len.ty)) %$(len.name) to i64")
                addcode(sb, "%p$varname = call i8* @malloc(i64 %i$varname)")
                func = "bitcast i8* %p$varname to $(strllvmtype(asm.val.T))*"
            end
            addcode(sb, "%$varname = $func")
        elseif typeof(asm) <: JLPAsmBr
            if isnothing(asm.condition)
                addcode(sb, "br label %$(asm.labels.label)")
            else
                addcode(sb, "br i1 %$(asm.condition.name), label %$(asm.labels[1].label), label %$(asm.labels[2].label)")
            end
        elseif typeof(asm) <: JLPAsmLabel
            append!(sb, "$(asm.label):\n")
        elseif typeof(asm) <: JLPAsmRet
            reurntype = asm.val.ty
            name = asm.val.name
            rettypestr = ""
            tr = typeof(returntype)
            if tr <: JLPPrimitive
                rettypestr *= "i$(string(returntype.bitlen))"
            end
            addcode(sb, "ret $rettypestr %$name")
        end
    end
    append!(sb, "}")
    return String(sb)
end

function generatellvmir(asms)
    sb = StringBuilder()
    for asm in asms
        ta = typeof(asm)
        if ta <: JLPLLVMFunc
            if asm.name == "Int_getchar" || asm.name == "Int_putchar_Int"
                continue;
            end
            emit(sb, generatellvmirfn(asm))
        elseif ta <:JLPLLVMDef
        end
    end
    String(sb)
end

function setall(ls::Vector)
    macrolist = []
    deflist = []
    primitivelist = []
    structlist = []
    otherlist = []
    for l in ls
        if l.car.val == "macro"
            push!(macrolist, l.cdr)
        end
    end
    for l in macrolist
        setmacro(l)
    end
    ls = expandallmacro.(ls)
    for l in ls
        if l.car.val == "def"
            push!(deflist, l.cdr)
        elseif l.car.val == "primitive"
            push!(primitivelist, l.cdr)
        elseif l.car.val == "struct"
            push!(structlist, l.cdr)
        else
            push!(otherlist, l)
        end
    end
    setprimitives(primitivelist)
    defs["Vector"] = JLPGenericStruct("Vector", [JLPAnyValue("T")], ["T" => JLPAnyValue("T")])
    for l in deflist
        setdef(l)
    end
    asms = []
    for l in deflist
        tmp = generateasmdef(l)
        push!(asms, tmp)
    end
"""
declare i8* @malloc(i64)
declare i32 @getchar()
declare i32 @putchar(i32)
define i32 @Int_getchar() {
  %1 = call i32 @getchar()
  ret i32 %1
}
define i32 @Int_putchar_Int(i32 %a) {
  %1 = call i32 @putchar(i32 %a)
  ret i32 %1
}
""" * generatellvmir(asms) * """
define i32 @main() {
  %1 = call i32 @Int_main()
  ret i32 %1
}"""
end


text:
"""
    half-quadruped
    q = (x, z, t, q1, ..., q8)
        x - lateral position
        z - vertical position
        t - body orientation
        q1 - leg 1 hip angle  (absolute)
        q2 - leg 1 knee angle (absolute)
        q5 - leg 3 hip angle  (absolute)
        q6 - leg 3 knee angle (absolute)
"""
mutable struct HalfQuadruped{T} <: Model{T}
    # dimensions
	nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters
    nc::Int # contact points

    # parameters

    # torso
    l_torso1::T
    l_torso2::T
    m_torso::T
    J_torso::T

    # leg 1
        # thigh
    l_thigh1::T
    lc_thigh1::T
    m_thigh1::T
    J_thigh1::T
        # calf
    l_calf1::T
    lc_calf1::T
    m_calf1::T
    J_calf1::T

	# leg 3
        # thigh
    l_thigh3::T
    lc_thigh3::T
    m_thigh3::T
    J_thigh3::T
        # calf
    l_calf3::T
    lc_calf3::T
    m_calf3::T
    J_calf3::T

    # joint friction
	friction_joint::Vector{T}
	joint_spring_stiffness::Vector{T}
	joint_spring_offset::Vector{T}

    # environment
    friction_body_world::Vector{T}
    friction_foot_world::Vector{T}
    gravity::T
end

function kinematics_hip(model::HalfQuadruped, q; hip=:none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if hip == :hip1
        l = model.l_torso1
    elseif hip == :hip2
        l = -model.l_torso2
    end

    return [
            x + l * cos(t_torso);
            z + l * sin(t_torso)
           ]
end

function kinematics_jacobian_hip(model::HalfQuadruped, q; hip=:none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if hip == :hip1
        l = model.l_torso1
    elseif hip == :hip2
        l = -model.l_torso2
    end

    jac = zeros(eltype(q), 2, model.nq)

    jac[1, 1] = 1.0
    jac[1, 3] = -l * sin(t_torso)
    jac[2, 2] = 1.0
    jac[2, 3] = l * cos(t_torso)

    return jac
end

function kinematics_thigh(model::HalfQuadruped, q; leg=:none, mode = :ee)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        t_shoulder = q[4]
        l_torso = model.l_torso1
        le_thigh = model.l_thigh1
        lc_thigh = model.lc_thigh1
        hip = :hip1
    elseif leg == :leg3
        t_shoulder = q[6]
        l_torso = -model.l_torso2
        le_thigh = model.l_thigh3
        lc_thigh = model.lc_thigh3
        hip = :hip2
    else
        @error "incorrect leg specified"
    end

    if mode == :ee
        l_thigh = le_thigh
    elseif mode == :com
        l_thigh = lc_thigh
    else
        @error "incorrect mode specified"
    end
    k_hip = kinematics_hip(model, q, hip=hip)

    return k_hip + [l_thigh * sin(t_shoulder);
                   -l_thigh * cos(t_shoulder)]
end

function kinematics_jacobian_thigh(model::HalfQuadruped, q; leg=:none, mode = :ee)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
        l_torso = model.l_torso1
        le_thigh = model.l_thigh1
        lc_thigh = model.lc_thigh1
        hip = :hip1
    elseif leg == :leg3
        idx = 6
        l_torso = -model.l_torso2
        le_thigh = model.l_thigh3
        lc_thigh = model.lc_thigh3
        hip = :hip2
    else
        @error "incorrect leg specified"
    end

    t_shoulder = q[idx]

    if mode == :ee
        l_thigh = le_thigh
    elseif mode == :com
        l_thigh = lc_thigh
    else
        @error "incorrect mode specified"
    end

    jac = kinematics_jacobian_hip(model, q, hip=hip)

    jac[1, idx] += l_thigh * cos(t_shoulder)
    jac[2, idx] += l_thigh * sin(t_shoulder)

    return jac
end

function kinematics_calf(model::HalfQuadruped, q; leg=:none, mode = :none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
        l_torso = model.l_torso1
        l_thigh = model.l_thigh1
        le_calf = model.l_calf1
        lc_calf = model.lc_calf1
    elseif leg == :leg3
        idx = 6
        l_torso = -model.l_torso2
        l_thigh = model.l_thigh3
        le_calf = model.l_calf3
        lc_calf = model.lc_calf3
    else
        @error "incorrect leg specified"
    end

    t_shoulder = q[idx]
    t_calf = q[idx + 1]

    if mode == :ee
        l_calf = le_calf
    elseif mode == :com
        l_calf = lc_calf
    else
        @error "incorrect mode specified"
    end

    k_thigh = kinematics_thigh(model, q, leg=leg, mode=:ee)

    return k_thigh + [l_calf * sin(t_calf);
                     -l_calf * cos(t_calf)]
end

function kinematics_jacobian_calf(model::HalfQuadruped, q; leg=:none, mode = :none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
        l_torso = model.l_torso1
        l_thigh = model.l_thigh1
        le_calf = model.l_calf1
        lc_calf = model.lc_calf1
    elseif leg == :leg3
        idx = 6
        l_torso = -model.l_torso2
        l_thigh = model.l_thigh3
        le_calf = model.l_calf3
        lc_calf = model.lc_calf3
    else
        @error "incorrect leg specified"
    end

    t_shoulder = q[idx]
    t_calf = q[idx + 1]

    if mode == :ee
        l_calf = le_calf
    elseif mode == :com
        l_calf = lc_calf
    else
        @error "incorrect mode specified"
    end

    jac = kinematics_jacobian_thigh(model, q, leg=leg, mode=:ee)
    jac[1, idx + 1] += l_calf * cos(t_calf)
    jac[2, idx + 1] += l_calf * sin(t_calf)

    return jac
end

function kinematics(model::HalfQuadruped, q)
	p_foot_1 = kinematics_calf(model, q, leg=:leg1, mode=:ee) # foot 1
	p_foot_3 = kinematics_calf(model, q, leg=:leg3, mode=:ee) # foot 3

    p_knee_1 = kinematics_thigh(model, q, leg=:leg1, mode=:ee) # knee 1
	p_knee_3 = kinematics_thigh(model, q, leg=:leg3, mode=:ee) # knee 3

    p_hip_1 = kinematics_hip(model, q, hip=:hip1) # hip 1
	p_hip_2 = kinematics_hip(model, q, hip=:hip2) # hip 2

	[
     p_foot_1; p_foot_3;
     p_knee_1; p_knee_3;
     p_hip_1; p_hip_2
    ]
end

# Lagrangian
function lagrangian(model::HalfQuadruped, q, q̇)
	L = 0.0

	# torso
	p_torso = q[1:2]
	v_torso = q̇[1:2]

	L += 0.5 * model.m_torso * transpose(v_torso) * v_torso
	L += 0.5 * model.J_torso * q̇[3]^2.0
	L -= model.m_torso * model.gravity * p_torso[2]

	# thigh 1
	p_thigh_1 = kinematics_thigh(model, q, leg=:leg1, mode=:com)
	J_thigh_1 = kinematics_jacobian_thigh(model, q, leg=:leg1, mode=:com)
	v_thigh_1 = J_thigh_1 * q̇

	L += 0.5 * model.m_thigh1 * transpose(v_thigh_1) * v_thigh_1
	L += 0.5 * model.J_thigh1 * q̇[4]^2.0
	L -= model.m_thigh1 * model.gravity * p_thigh_1[2]

	# calf 1
	p_calf_1 = kinematics_calf(model, q, leg=:leg1, mode=:com)
	J_calf_1 = kinematics_jacobian_calf(model, q, leg=:leg1, mode=:com)
	v_calf_1 = J_calf_1 * q̇

	L += 0.5 * model.m_calf1 * transpose(v_calf_1) * v_calf_1
	L += 0.5 * model.J_calf1 * q̇[5]^2.0
	L -= model.m_calf1 * model.gravity * p_calf_1[2]

	# thigh 3
	p_thigh_3 = kinematics_thigh(model, q, leg=:leg3, mode=:com)
	J_thigh_3 = kinematics_jacobian_thigh(model, q, leg=:leg3, mode=:com)
	v_thigh_3 = J_thigh_3 * q̇

	L += 0.5 * model.m_thigh3 * transpose(v_thigh_3) * v_thigh_3
	L += 0.5 * model.J_thigh3 * q̇[6]^2.0
	L -= model.m_thigh3 * model.gravity * p_thigh_3[2]

	# calf 3
	p_calf_3 = kinematics_calf(model, q, leg=:leg3, mode=:com)
	J_calf_3 = kinematics_jacobian_calf(model, q, leg=:leg3, mode=:com)
	v_calf_3 = J_calf_3 * q̇

	L += 0.5 * model.m_calf3 * transpose(v_calf_3) * v_calf_3
	L += 0.5 * model.J_calf3 * q̇[7]^2.0
	L -= model.m_calf3 * model.gravity * p_calf_3[2]

	t_torso = q[3]
	t_hip1 = q[4]
	t_knee1 = q[5]
	t_hip2 = q[6]
	t_knee2 = q[7]
	L -= 0.5 * joint_spring_stiffness[1] * (joint_spring_offset[1] + t_torso - t_hip1)^2
	L -= 0.5 * joint_spring_stiffness[2] * (joint_spring_offset[2] + t_torso - t_hip2)^2
	L -= 0.5 * joint_spring_stiffness[3] * (joint_spring_offset[3] + t_hip1 - t_knee1)^2
	L -= 0.5 * joint_spring_stiffness[4] * (joint_spring_offset[4] + t_hip2 - t_knee2)^2
	return L
end

function signed_distance(model::HalfQuadruped, q)
	k = kinematics(model, q)
    return [k[2], k[4], k[6], k[8], k[10], k[12]]
end

function input_jacobian(::HalfQuadruped, q)
	[1.0  0.0  0.0  0.0  0.0  0.0  0.0; # infeasible controls
	 0.0  1.0  0.0  0.0  0.0  0.0  0.0; # infeasible controls
	 0.0  0.0  1.0  0.0  0.0  0.0  0.0; # infeasible controls
     0.0  0.0 -1.0  1.0  0.0  0.0  0.0;
     0.0  0.0  0.0 -1.0  1.0  0.0  0.0;
     0.0  0.0 -1.0  0.0  0.0  1.0  0.0;
     0.0  0.0  0.0  0.0  0.0 -1.0  1.0]
end

function contact_jacobian(model::HalfQuadruped, q)
    J1 = kinematics_jacobian_calf(model, q, leg=:leg1, mode=:ee)
	J3 = kinematics_jacobian_calf(model, q, leg=:leg3, mode=:ee)

    J5 = kinematics_jacobian_thigh(model, q, leg=:leg1, mode=:ee)
	J7 = kinematics_jacobian_thigh(model, q, leg=:leg3, mode=:ee)

    J9  = kinematics_jacobian_hip(model, q, hip=:hip1)
    J10 = kinematics_jacobian_hip(model, q, hip=:hip2)

    [
     J1;
     J3;
     J5;
     J7;
     J9;
     J10
    ]
end

# nominal configuration
function nominal_configuration(model::HalfQuadruped)
    [0.0; 0.5; 0.0 * π; 0.25 * π; 0.5 * π; -0.25 * π; 0.1 * π]
end

# friction coefficients
function friction_coefficients(model::HalfQuadruped)
    [model.friction_foot_world; model.friction_body_world]
end

# dimensions
nq = 3 + 2 * 2            # generalized coordinates
nu = 3 + 2 * 2            # controls
nw = 0                    # parameters
nc = 6                    # contact points

# parameters
gravity = 9.81            # gravity
friction_body_world = [0.5; 0.5; 0.5; 0.5] # coefficient of friction
friction_foot_world = [0.5; 0.5] # coefficient of friction
friction_joint = [0.1; 0.1; 0.1; 0.1]
joint_spring_stiffness = [1.0, 1.0, 1.0, 1.0]
joint_spring_offset = [0.0, 0.0, 0.0, 0.0]

# similar to Unitree A1
m_torso = 4.713 + 4 * 0.696
m_thigh = 1.013
m_calf = 0.166

J_torso = 0.01683 + 4 * 0.696 * 0.183^2.0
J_thigh = 0.00552
J_calf = 0.00299

l_torso1 = 0.183
l_torso2 = 0.183
l_thigh = 0.2
l_calf = 0.2

lc_thigh = 0.5 * l_thigh - 0.00323
lc_calf = 0.5 * l_calf - 0.006435

halfquadruped = HalfQuadruped(nq, nu, nw, nc,
				l_torso1, l_torso2, m_torso, J_torso,

				l_thigh, lc_thigh, m_thigh, J_thigh,
				l_calf, lc_calf, m_calf, J_calf,

				l_thigh, lc_thigh, m_thigh, J_thigh,
				l_calf, lc_calf, m_calf, J_calf,

	            friction_joint,
				joint_spring_stiffness,
				joint_spring_offset,
                friction_body_world,
                friction_foot_world,
                gravity)

halfquadruped_contact_kinematics = [
	q -> kinematics_calf(halfquadruped, q, leg=:leg1, mode=:ee),
    q -> kinematics_calf(halfquadruped, q, leg=:leg3, mode=:ee),
    q -> kinematics_thigh(halfquadruped, q, leg=:leg1, mode=:ee),
    q -> kinematics_thigh(halfquadruped, q, leg=:leg3, mode=:ee),
    q -> kinematics_hip(halfquadruped, q, hip=:hip1),
    q -> kinematics_hip(halfquadruped, q, hip=:hip2)]

halfquadruped_contact_kinematics_jacobians = [
 	q -> kinematics_jacobian_calf(halfquadruped, q, leg=:leg1, mode=:ee),
    q -> kinematics_jacobian_calf(halfquadruped, q, leg=:leg3, mode=:ee),
    q -> kinematics_jacobian_thigh(halfquadruped, q, leg=:leg1, mode=:ee),
    q -> kinematics_jacobian_thigh(halfquadruped, q, leg=:leg3, mode=:ee),
    q -> kinematics_jacobian_hip(halfquadruped, q, hip=:hip1),
    q -> kinematics_jacobian_hip(halfquadruped, q, hip=:hip2)]

name(::HalfQuadruped) = :halfquadruped

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export checkbox

"""
    checkbox(;kwargs...)

A Checkbox component.
Checkbox Input
Keyword arguments:
- `id` (String; required): Provide an `id` to uniquely identify the Checkbox input
- `label` (a list of or a singular dash component, string or number; required): Provide a label to provide a description of the Checkbox input that you are
exposing to the user
- `hideLabel` (Bool; optional): Specify whether the label should be hidden, or not
- `disabled` (Bool; optional): Specify whether the Checkbox should be disabled
- `value` (Bool; optional): Specify whether the underlying input is checked
- `title` (String; optional): Specify a title for the <label> node for the Checkbox
- `style` (Dict; optional): The inline styles
- `className` (String; optional): Specify an optional className to be applied to the <label> node
"""
function checkbox(; kwargs...)
        available_props = Symbol[:id, :label, :hideLabel, :disabled, :value, :title, :style, :className]
        wild_props = Symbol[]
        return Component("checkbox", "Checkbox", "dash_carbon_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
# Simple macro
macro hello()
    return :(
    for i in 1:3
        println("hello world")
    end
    )
end
#=
julia> macro hello()
           return :(
           for i in 1:3
               println("hello world")
           end
           )
       end

julia> @hello()
hello world
hello world
hello world

julia> @hello
hello world
hello world
hello world
=#

# -----------------------------------------------------------------------------
# passing literal argument
macro hello(n)
    return :(
    for i in 1:$n
        println("hello world")
    end
    )
end
#=
julia> macro hello(n)
           return :(
           for i in 1:$n
               println("hello world")
           end
           )
       end
@hello (macro with 2 methods)

julia> @hello(2)
hello world
hello world

julia> @hello 2
hello world
hello world
=#

# -----------------------------------------------------------------------------
# Passing expression arguments
function showme(x)
    @show x
end
#=
julia> a = 1; b = "hello"; c = :hello;

julia> function showme(x)
           @show x
       end

julia> showme(a);
x = 1

julia> showme(b);
x = "hello"

julia> showme(c);
x = :hello
=#

macro showme(x)
    @show x
end
#=
julia> macro showme(x)
           @show x
       end

julia> @showme(a);
x = :a

julia> @showme(b);
x = :b

julia> @showme(c);
x = :c
=#

# -----------------------------------------------------------------------------

# macro expansion
#=
julia> @macroexpand @hello 2
:(for #4#i = 1:2
      #= REPL[3]:4 =#
      Main.println("hello world")
  end)
=#

# what happens?
macro identity(ex)
    dump(ex)
    return ex
end
#=
julia> macro identity(ex)
           dump(ex)
           return ex
       end
=#

function foo()
    return @identity 1 + 2 + 3
end

#=
julia> function foo()
           return @identity 1 + 2 + 3
       end
Expr
  head: Symbol call
  args: Array{Any}((4,))
    1: Symbol +
    2: Int64 1
    3: Int64 2
    4: Int64 3
foo (generic function with 1 method)
=#

@code_lowered foo()
#=
julia> @code_lowered foo()
CodeInfo(
1 ─ %1 = 1 + 2 + 3
└──      return %1
)
=#

# -----------------------------------------------------------------------------
# manipulating expressions

# example 1
macro squared(ex)
    return :($(ex) * $(ex))
end
#=
julia> macro squared(ex)
           return :($(ex) * $(ex))
       end
@squared (macro with 1 method)

julia> @squared 3
9
=#

# does not quite work properly
function foo()
    x = 2
    return @squared x
end
#=
julia> function foo()
           x = 2
           return @squared x
       end

julia> foo()
ERROR: UndefVarError: x not defined
=#

# why?
@code_lowered foo()

#=
julia> @code_lowered foo()
CodeInfo(
1 ─      x = 2
│   %2 = Main.x * Main.x
└──      return %2
)
=#

# let's fix it
macro squared(ex)
    return :($(esc(ex)) * $(esc(ex)))
end
#=
julia> macro squared(ex)
           return :($(esc(ex)) * $(esc(ex)))
       end
=#

function foo()
    x = 2
    return @squared x
end

#=
julia> function foo()
           x = 2
           return @squared x
       end

julia> foo()
4
=#

# example 2
#=
julia> :(sin(x)) |> dump
Expr
  head: Symbol call
  args: Array{Any}((2,))
    1: Symbol sin
    2: Symbol x

julia> :(sin(sin(x))) |> dump
Expr
  head: Symbol call
  args: Array{Any}((2,))
    1: Symbol sin
    2: Expr
      head: Symbol call
      args: Array{Any}((2,))
        1: Symbol sin
        2: Symbol x
=#

macro compose_twice(ex)
    @assert ex.head == :call
    @assert length(ex.args) == 2
    me = copy(ex)
    ex.args[2] = me
    return ex
end
@compose_twice(sin(1)) == sin(sin(1))

#=
julia> macro compose_twice(ex)
           @assert ex.head == :call
           @assert length(ex.args) == 2
           me = copy(ex)
           ex.args[2] = me
           return ex
       end

julia> @compose_twice(sin(1)) == sin(sin(1))
true
=#

# -----------------------------------------------------------------------------

# macro hygiene
macro ntimes(n, ex)
    quote
        times = $(esc(n))
        for i in 1:times
            $(esc(ex))
        end
    end
end
#=
julia> macro ntimes(n, ex)
           quote
               times = $(esc(n))
               for i in 1:times
                   $(esc(ex))
               end
           end
       end
=#

function foo()
    times = 0
    @ntimes 3 println("hello world")
    println("times = ", times)
end
#=
julia> function foo()
           times = 0
           @ntimes 3 println("hello world")
           println("times = ", times)
       end

julia> foo()
hello world
hello world
hello world
times = 0
=#

# let's expand the macro and see what happens
@macroexpand(@ntimes 3 println("hello world"))

#=
julia> @macroexpand(@ntimes 3 println("hello world"))
quote
    #= REPL[73]:3 =#
    #118#times = 3
    #= REPL[73]:4 =#
    for #119#i = 1:#118#times
        #= REPL[73]:5 =#
        println("hello world")
    end
end
=#

# -----------------------------------------------------------------------------
# non standard string literal

using DataFrames

# regex
#=
julia> typeof(r"^hello")
Regex
=#

# our own
#=
julia> DataFrame(x1 = rand(Float64, 100000), x2 = rand(Int16, 100000))
=#

# nmeric data frame non-standard string literal macro
macro ndf_str(s)
    nstr, spec = split(s, ":")
    n = parse(Int, nstr) # number of rows
    types = split(spec, ",") # column type specifications

    num_columns = length(types)

    mappings = Dict(
    "f64"=>Float64, "f32"=>Float32,
    "i64"=>Int64, "i32"=>Int32, "i16"=>Int16, "i8"=>Int8)

    column_types = [mappings[t] for t in types]
    column_names = [Symbol("x$i") for i in 1:num_columns]

    return DataFrame([column_names[i] => rand(column_types[i], n)
        for i in 1:num_columns]...)
end

ndf"100000:f64,f32,i16,i8"
#=
julia> ndf"100000:f64,f32,i16,i8"
100000×4 DataFrame
│ Row │ x1       │ x2        │ x3    │ x4   │
│     │ Float64  │ Float32   │ Int16 │ Int8 │
├─────┼──────────┼───────────┼───────┼──────┤
│ 1   │ 0.160496 │ 0.404637  │ 6805  │ 96   │
│ 2   │ 0.994738 │ 0.636778  │ 23910 │ -59  │
│ 3   │ 0.59626  │ 0.520906  │ 11857 │ -71  │
│ 4   │ 0.637577 │ 0.152872  │ -3842 │ 48   │
│ 5   │ 0.396112 │ 0.0520195 │ -2919 │ -111 │
=#


text:
module BlueMesh7
using Agents, AgentsPlots
using LightGraphs, GraphPlot
using Random
using StatsBase
include("MeshAgents.jl")
using .MeshAgents: Node, Source, Packet
using Distributions
using LoopVectorization
using SpecialFunctions

export generate_positions, initialize_mesh, run, plotgrid

const ble_37_channel_wavelength = Float64(0.12491352)
const part_of_path_loss_for_37 = Float64(20 * log(10, 4 * π / ble_37_channel_wavelength))
const Standard_Normal_d = Normal()

function agent_step!(source::Source, model::AgentBasedModel)
    if rand() < model.packet_emit_rate
        # move off to a random position
        move_agent!(source, model)

        source.packet = Packet(model.tick, source.id, rand(2:nagents(model)), model.ttl)
        source.event_start = model.tick

        model.produced += 1
    end

    source.packet === nothing && return
    model.tick < source.event_start && return

    d, r = divrem(model.tick - source.event_start, source.t_interpdu)
    source.transmitting = r == 0
    source.channel = d + 37

    if source.channel > 39
        if source.n_transmit_left > 0
            source.event_start = model.tick + source.t_og_transmit_delay + rand(source.rt_og_random_delay)
            source.n_transmit_left -= 1
            return
        end

        source.packet = nothing
        source.transmitting = false
        source.n_transmit_left = source.n_og_transmit_count
    end
end

function agent_step!(node::Node, model::AgentBasedModel)
    if node.state == :scanning
        # scanning happens every tick (assuming here scan_window == scan_interval)
        node.channel = div(model.tick - node.event_start, node.t_scan_interval) + 37

        # no packets acquired, repeat
        if node.channel > 39
            node.channel = 37
            node.event_start = model.tick
        end

        # instead of type instability here, zero(Packet) might work
        node.packet === nothing && return

        if node.packet.seq in node.received_packets
            node.packet = nothing

            return
        end

        push!(node.received_packets, node.packet.seq)

        if node.packet.dst == node.id
            model.received += 1
            node.packet = nothing

            return
        end

        if node.packet.ttl == 0
            node.packet = nothing

            return
        end

        node.packet.ttl -= 1

        if node.role == :relay
            # initiating back-off delay before advertising this packet
            node.event_start = model.tick + node.t_back_off_delay
            node.state = :advertising

            return
        end
    end

    node.role != :relay && return

    if node.state == :advertising
        model.tick < node.event_start && return

        # return to scanning if the packet is absent
        if node.packet === nothing
            node.state = :scanning
            node.event_start = model.tick

            return
        end

        d, r = divrem(model.tick - node.event_start, node.t_interpdu)

        # transmitting happens once for each channel
        node.transmitting = iszero(r)

        node.channel = d + 37

        # the end of the advertising
        if node.channel > 39
            # bonus retransmissions for the original/retx packet
            if node.n_transmit_left > 0
                node.event_start = model.tick + node.t_retx_transmit_delay + rand(node.rt_retx_random_delay)

                node.n_transmit_left -= 1
                node.channel = 37
                return
            end

            node.packet = nothing
            node.transmitting = false
            node.state = :scanning
            node.event_start = model.tick + 1
            node.n_transmit_left = node.n_retx_transmit_count
            node.channel = 37
        end
    end
end

distance(a::Tuple{Int, Int}, b::Tuple{Int, Int}) = sqrt((a[1] - b[1]) ^ 2 + (a[2] - b[2]) ^ 2)

function vmap!(arr::SubArray, func::Function, arr2::SubArray)
    length(arr) > 0 || return Array{typeof(func(arr[1])), 1}(undef, 0)
    @avx for i in 1:length(arr)
        arr2[i] = func(arr[i])
    end
    arr2
end

function vmap!(arr::SubArray, func::Function)
    length(arr) > 0 || return Array{typeof(func(arr[1])), 1}(undef, 0)
    arr2 = Array{typeof(func(arr[1])), 1}(undef, length(arr))
    @avx for i in 1:length(arr)
        arr2[i] = func(arr[i])
    end
    arr2
end

function calc_rssi!(agents::Array{Union{Node, Source}, 1}, P_tx::Array{Float64, 1})
    n = length(agents)
    rssi_map = Matrix{Float64}(undef, n, n)

    L_p = Matrix{Float64}(undef, n, n)

    for i in 1:n
        for j in 1:(i - 1) # more cache-friendly
            L_p[i, j] = L_p[j, i]
        end

        for j in (i + 1):n
            L_p[i, j] = 1 / distance(agents[i].pos, agents[j].pos) ^ 2
        end
    end

    @avx for i in 1:n
        for j in 1:n
            rssi_map[i, j] = P_tx[j] * L_p[i, j] * part_of_path_loss_for_37 # dB
        end

        rssi_map[i, i] = 0
    end

    rssi_map
end

function recalc_sources_rssi!(agents::Array{Union{Node, Source}, 1}, sources::Array{Int, 1},
        rssi_map:: Matrix{Float64}, P_tx::Array{Float64, 1})
    n = length(agents)

    for i in sources
        Lps = vmap!(view(agents, 1:n), x -> 1 / distance(agents[i].pos, x.pos) ^ 2)

        @avx for j in 1:n
            Lps[j] = Lps[j] * part_of_path_loss_for_37
            rssi_map[i, j] = P_tx[j] * Lps[j] # dB
        end

        @avx for j in 1:n
            rssi_map[j, i] = P_tx[i] * Lps[j] # dB
        end

        rssi_map[i, i] = 0

    end

    rssi_map
end

function model_step!(model::AgentBasedModel)
    all_agents = collect(allagents(model))

    transmitters = [
            [a.transmitting && a.channel == 37 for a in all_agents],
            [a.transmitting && a.channel == 38 for a in all_agents],
            [a.transmitting && a.channel == 39 for a in all_agents]
    ]

    recalc_sources_rssi!(all_agents, model.source_nodes, model.rssi_map, model.tx_powers)

    for (agent_i, agent) in enumerate(all_agents)
        # since we support just 1 message within scanning period
        (agent.state == :scanning && agent.packet == nothing) || continue

        neighbours = model.rssi_map[agent_i, transmitters[agent.channel - 36]]
        neighbours_agents = all_agents[transmitters[agent.channel - 36]]
        length(neighbours) > 0 || continue

        # add random noises to calculate RSSI
        neighbours = map(n -> n * rand(model.shadow_d) * rand(model.multipath_d), neighbours)

        visible_mask = [rssi > model.scanner_sensitivity for rssi in neighbours]
        neighbours = neighbours[visible_mask]
        neighbours_agents = neighbours_agents[visible_mask]

        length(neighbours) > 0 || continue

        # add some noises
        total = sum(neighbours) + max(rand(model.wifi_noise_d), 0) + max(rand(model.gaussian_noise), 0)
        #println("RSSI", neighbours, total)

        # calculate SINR and take sqrt
        neighbours = map(n -> n / (total - n), neighbours)

        # calculate BER
        neighbours = map(n -> 0.5 * erfc(sqrt(n / 2)), neighbours)

        # calculate PAR (Package Accept Rate) (PAR = 1 - PER)
        neighbours = map(n -> (1 - n)^model.msg_bit_length, neighbours)

        total_PAR = sum(neighbours)

        rand() <= total_PAR || continue # then we will receive message

        tx_i = StatsBase.sample(1:length(neighbours), Weights(neighbours))
        agent.packet = deepcopy(neighbours_agents[tx_i].packet)
    end

    model.tick += 1
end

function plotgrid(model::AgentBasedModel)
    g = SimpleGraph(nagents(model))

    colors = Dict(:relay => "orange", :sink => "darkblue", :source => "transparent")

    agents = sort(collect(allagents(model)), by=a -> a.id)
    xs = map(a -> a.pos[1], agents)
    ys = map(a -> a.pos[2], agents)
    cs = map(a -> colors[a.role], agents)

    gplot(g, xs, ys; nodefillc = cs)
end

function plotgraph(model::AgentBasedModel)
    g = SimpleGraph(nagents(model))

    for (src_i, src) in enumerate(allagents(model)), (dst_i, dst) in enumerate(allagents(model))
        src.role == :relay || continue
        dst.role != :source || continue

        if model.rssi_map[src_i, dst_i] * mean(model.multipath_d) * mean(model.shadow_d) > model.scanner_sensitivity
            add_edge!(g, src.id, dst.id)
        end
    end

    colors = Dict(:relay => "orange", :sink => "darkblue", :source => "transparent")

    agents = sort(collect(allagents(model)), by=a -> a.id)
    xs = map(a -> a.pos[1], agents)
    ys = map(a -> a.pos[2], agents)
    cs = map(a -> colors[a.role], agents)

    gplot(g, xs, ys; nodefillc = cs)
end

"""
    generate_positions(; dims = (200, 200), n_nodes = 100) → positions::Array{Tuple{Int,Int},1}

Generates random xy-coordinates of the nodes
"""
generate_positions(; dims = (100, 100), n_nodes = 64) = [(rand(1:dims[1]), rand(1:dims[2])) for _ in 1:n_nodes]

"""
    initialize_mesh(positions::Vector{Tuple{Int, Int}}, roles::Vector{Int}) → mesh::AgentBasedModel

Create a mesh by specifying roles of each node by `roles`=[0, 1, 0...] where 1 means that the node is a relay, sink otherwise, and including `positions`=[(Int, Int)...] as nodes' positions
"""
to_mW(dBm::Int64) = 10^(dBm / 10)
function initialize_mesh(positions::Vector{Tuple{Int, Int}}, roles::Vector{Int})
    length(roles) == length(positions) || throw(ArgumentError("Both arguments must be of equal length"))

    properties = Dict(
        :tick                   => 0,
        :packet_error_rate      => 0.05,
        :packet_emit_rate       => 10 / 1000,
        :ttl                    => 4,
        :msg_bit_length         => 312,
        :received               => 0,
        :produced               => 0,
        :rssi_map               => Matrix{Float64}(undef, 0, 0),
        :source_nodes           => Array{Int64, 1}(),
        :tx_powers              => Array{Float64, 1}(undef, length(roles)),
        :scanner_sensitivity    => to_mW(-95),
        :shadow_d               => LogNormal((to_mW(0) - to_mW(-4)) / 2 + to_mW(-4), (to_mW(0) - to_mW(-4)) / 2),
        :multipath_d            => Rayleigh(to_mW(-4)),
        :wifi_noise_d           => Normal((to_mW(-125) - to_mW(-135)) / 2 + to_mW(-135), (to_mW(-125) - to_mW(-135)) / 2),
        :gaussian_noise         => Normal((to_mW(-125) - to_mW(-135)) / 2 + to_mW(-135), (to_mW(-125) - to_mW(-135)) / 2)
    )

    space = GridSpace((maximum(first.(positions)), maximum(last.(positions))), moore = true)
    model = ABM(Union{Source, Node}, space; properties, warn = false)

    model.source_nodes = [1]
    source = Source(id = 1, pos = (1, 1))
    add_agent_single!(source, model)
    model.tx_powers[1] = source.tx_power

    for i in 2:length(roles)
        node = Node(id = i; pos = positions[i], role = roles[i] == 1 ? :relay : :sink)

        add_agent_pos!(node, model)

        model.tx_powers[i] = node.tx_power
    end

    model.rssi_map = calc_rssi!(collect(allagents(model)), model.tx_powers)

    #println("model.rssi_map: ", model.rssi_map)
    model
end

"""
    start(model::AgentBasedModel, minutes = 1) → (received, produced)

Start an experiment from `model` for the given number of minutes.
"""
function start(model::AgentBasedModel; minutes = 1)
    steps = minutes * 60 * 1000

    _, dfm = run!(model, agent_step!, model_step!, steps; when_model = [steps], mdata = [:produced, :received])

    dfm[:produced][1], dfm[:received][1]
end
end


text:
using LinearAlgebra
using Zygote
using ForwardDiff
using ReverseDiff
using DiffEqOperators
using StaticArrays
using Base.Iterators
using Plots

## Stencil generation
neighbours(CI, deg)        = CI .+ CartesianIndices(ntuple(x -> -deg:deg, length(CI)))

stencil_depth(CI, depth)   = [ CI + CartesianIndex(ntuple(x -> ifelse(x == i, depth, 0), length(CI))) for i in 1:length(CI) ]

stencils(CI, deg, poke = false) = ifelse(poke, [ stencil_depth(CI, -(deg + 1)); neighbours(CI, deg)[:]; stencil_depth(CI, deg + 1) ], neighbours(CI, deg)[:])

von_neumann_stencil(CI, n) = stencils(CI, n, true)
moore_stencil(CI, n)       = stencils(CI, n, false)
cs_stencil(φ, index)       = @view φ[von_neumann_stencil(index, 0)]

## Cells and boundaries
abstract type AbstractCell end

# This structure can be encoded generically using the Store comonad; here we define it explicitly for testing.
struct CartesianPoint{M,N,T} <: AbstractCell
    index  :: CartesianIndex{N}
    values :: SVector{M,T}
end

von_neumann_stencil(ci :: CartesianPoint, n) = von_neumann_stencil(ci.index, n)
moore_stencil(ci :: CartesianPoint, n)       = moore_stencil(ci.index, n)

struct CartesianBoundary{M,N,T} <: AbstractCell
    point      :: CartesianPoint{M,N,T}
    neighbours :: Vector{Tuple{CartesianIndex{N},SVector{M,T}}}
end

## Finite differencing operations on values and stencils

# First derivatives
Δ_central(φp, φm, Δx)     = (φp - φm) / 2Δx
Δ_forward(φ1, φ2, Δt)     = (φ2 - φ1) / Δt
central_2nd(φp, φ, φm)    = φm - 2φ + φp

# Stencil
Δx_central(φ, dx)  = Δ_central(φ[4], φ[1], dx)
Δy_central(φ, dy)  = Δ_central(φ[5], φ[2], dy)

Δx_backward(φ, dx) = Δ_backward(φ[3], φ[1], dx)
Δy_backward(φ, dy) = Δ_backward(φ[3], φ[2], dy)

Δx_forward(φ, dx)  = Δ_forward(φ[4], φ[3], dx)
Δy_forward(φ, dy)  = Δ_forward(φ[5], φ[3], dy)

# Second derivatives
Δ²_central(φp, φ, φm, Δx) = central_2nd(φp, φ, φm) / Δx^2
Δ²_central(φp, φ, φm, Δxp, Δxm) = Δ_forward(φp, φ, Δxp) - Δ_forward(φ, φm, Δxm)

# Stencil
Δ²x_central(φ, dx) = Δ²_central(φ[4], φ[3], φ[1], dx)
Δ²y_central(φ, dy) = Δ²_central(φ[5], φ[3], φ[2], dy)

Δ²x_central(φ, dxs) = Δ²_central(φ[4], φ[3], φ[1], dxs[1], dxs[2])
Δ²y_central(φ, dys) = Δ²_central(φ[5], φ[3], φ[2], dys[1], dys[2])

# Boundary (this needs type distinctions between boundary and non-boundary cells to be efficient)
function substitute(φ, index, val)
    if index ∈ CartesianIndices(φ)
        φ[index]
    else
        val
    end
end

## Governing equations and residual setup
temperature_residual(φ, k, q, ds) = -k * (Δ²x_central(φ, ds[1:2]) + Δ²y_central(φ, ds[3:4])) + q

# Ruleset for thermal diffusion with a source field
function ruleset(φs, σs, k, dxs, dys, index, φ_boundaries)
    # Boundary conditions check
    if any(index.I .<= 1) || any(index.I .>= size(φs))
        # There have to be more conditions generally...
        φs_vec = [ substitute(φs, local_index, φ_boundaries) for local_index in von_neumann_stencil(index, 0) ]
        ds_vec = [ substitute(dxs, local_index, 0.1)          for local_index in von_neumann_stencil(index, 0) ]
    else
        # Get neighbours of a cell and corresponding sizes
        φs_vec   = cs_stencil(φs, index)
        dxp, dxn = dxs[index], dxs[index + CartesianIndex(1,0)]
        dyp, dyn = dys[index], dys[index + CartesianIndex(0,1)]
        ds_vec   = [ dxp; dxn; dyp; dyn ]
    end

    σ = σs[index]

    # Evaluate residuals
    temperature_residual(φs_vec, k, σ, ds_vec)
end

# Compute residual equations
compute_residuals(p, σ, k, dxs, dys, φ_boundary) = map(index -> ruleset(p, σ, k, dxs, dys, index, φ_boundary), CartesianIndices(p))

newton_update(p, δp, α = 1.0) = p - α * δp

# Pre-allocated version
newton_update!(p, δp, α = 1.0) = @. p = p - α * δp

## Initial setup
nx, ny  = 50, 40

# Less allocations than broadcasting inside sum
polynomial(x, coeffs) = sum(k -> k[1] * x^k[2], zip(coeffs, 0:length(coeffs)-1))

# Grid generation
xs      = range(0, 1, length = nx)
ys      = polynomial.(xs, Ref(ones(5)))
grid    = reshape([ SVector(x, y) for (x, y_max) in zip(xs, ys) for y in range(y_max, -y_max, length = 2ny) ], 2ny, nx)

# Compute grid spacings
dxs = norm.(grid[:,2:end] .- grid[:,1:end-1])
dys = norm.(grid[2:end,:] .- grid[1:end-1,:])

##
α       = 1.0
k       = 1.
φ_bound = 200

φ0      = φ_bound .* ones(size(grid)...)
σ       = rand(size(grid)...)

p       = φ0
R       = similar(p[:])
# ∂R∂p = zeros(2 .* (prod(size(ωs)), prod((reverse ∘ size)(ωs)))...)

num_iters = 3
compute_residuals!(R, p) = R .= compute_residuals(p, σ, k, dxs, dys, φ_bound)[:]
compute_grad_fwd(R, p) = ForwardDiff.jacobian(compute_residuals!, R, p)

##
# using NLsolve
# prob = nlsolve(compute_residuals!, p)

## NEWTON ITERATIONS
#=================================================================#

function newton_solver!(R, p, num_iters = 3, α = 1.0)
    # Array to store errors
    ε = zeros(num_iters);

    # Newton iteration loop
    for i in 1:num_iters
        # Compute residuals
        R    = compute_residuals!(R, p)

        # Compute Jacobian
        ∂R∂p = compute_grad_fwd(R, p)

        # Compute Newton step
        Δp   = ∂R∂p \ -R

        # Update state with relaxation factor
        p    = newton_update!(p, reshape(Δp, size(p)), α)

        # Error processing
        ε[i] = maximum(abs.(Δp))
        println("Newton step error L² norm: $(ε[i])")
    end
    R, p, ε
end

##
R, p, ε = newton_solver!(R, p, num_iters, α)

max_stars_repo_name: GodotMisogi/AdjointCFD

max_issues_repo_name: GodotMisogi/AdjointCFD

max_forks_repo_name: GodotMisogi/AdjointCFD


text:
"""
# Arguments
- `S::Int` - Current Susceptible Population
- `I::Vector{Int}` - Current Infected Population
- `R::Int` - Current Recovered Population
- `T::Matrix{Int}` - Array for which people belonging to array element ``T_{i,j}`` are ``i-1`` days away
    from receiving positive test and have infection age ``j``
"""
struct SingleCovidState <: CovidPOMDPs.CovidState
    S::Int # Current Susceptible Population
    I::Vector{Int} # Current Infected Population
    R::Int # Current Recovered Population
    T::Matrix{Int} # Rows: Days from receiving test result; Columns: Infection Age
    params::InfParams
    prev_action::CovidAction
end

CovidPOMDPs.CovidState(S,I,R,T,p,a) = SingleCovidState(S,I,R,T,p,a)

CovidPOMDPs.statevars(s::SingleCovidState) = (s.S, s.I, s.R, s.T)

CovidPOMDPs.susceptible(s::SingleCovidState) = s.S
CovidPOMDPs.infected(s::SingleCovidState) = sum(s.I)
CovidPOMDPs.recovered(s::SingleCovidState) = s.R

Base.@kwdef struct SingleCovidPOMDP{A} <: POMDP{SingleCovidState, CovidAction, Int}
    "Delay (in days) between test being administered and result of test being received `(≥ 0)`"
    test_delay::Int = 1

    "Total population count `(> 0)`"
    N::Int = 10^6

    "POMDP discount factor `(γ ∈ [0,1])`"
    discount::Float64 = 0.95

    "Weight with which to penalize new infections `(≥ 0.0)`"
    inf_loss::Float64 = 1.0

    "Weight with which to penalize testing rate `(≥ 0.0)`"
    test_loss::Float64 = 1.0

    "Weight with which to penalize changes in testing rate `(≥ 0.0)`"
    testrate_loss::Float64 = 1.0

    "Number of days for which a testing policy must be held `(≥ 1)`"
    test_period::Int = 1

    "CovidPOMDP action space (default is continuous [0,1])"
    actions::A = CovidPOMDPs.CovidActionSpace()
end

"""
Take given CovidPOMDP obj and return same CovidPOMDP obj only with test_period changed to 1
"""
function CovidPOMDPs.unity_test_period(pomdp::SingleCovidPOMDP)
    return SingleCovidPOMDP(
        pomdp.test_delay,
        pomdp.N,
        pomdp.discount,
        pomdp.inf_loss,
        pomdp.test_loss,
        pomdp.testrate_loss,
        1,
        pomdp.actions
    )
end

function CovidPOMDPs.init_SIRT(pomdp::SingleCovidPOMDP, rng::AbstractRNG=Random.GLOBAL_RNG)
    N = pomdp.N
    S, inf, R = floor.(Int, simplex_sample(3, Float64(N), rng))

    horizon = INFECTION_HORIZON

    I = floor.(Int, simplex_sample(horizon, Float64(inf), rng))

    leftover = N - (S + sum(I) + R)
    R += leftover

    T = zeros(Int, pomdp.test_delay+1, horizon)

    return S, I, R, T
end

function CovidPOMDPs.init_SIRT(pomdp::SingleCovidPOMDP, inf::Int, rng::AbstractRNG=Random.GLOBAL_RNG)
    N = pomdp.N
    R = 0

    horizon = INFECTION_HORIZON

    I = floor.(Int, simplex_sample(horizon, Float64(inf), rng))

    leftover = N - sum(I)
    S = leftover

    T = zeros(Int, pomdp.test_delay+1, horizon)

    return S, I, R, T
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

using LinearAlgebra

function _last_primal_solution(model::Model)
    if !has_values(model)
        error(
            "No primal solution is available. You must provide a point at " *
            "which to check feasibility.",
        )
    end
    return Dict(v => value(v) for v in all_variables(model))
end

"""
    primal_feasibility_report(
        model::Model,
        point::AbstractDict{VariableRef,Float64} = _last_primal_solution(model),
        atol::Float64 = 0.0,
        skip_missing::Bool = false,
    )::Dict{Any,Float64}

Given a dictionary `point`, which maps variables to primal values, return a
dictionary whose keys are the constraints with an infeasibility greater than the
supplied tolerance `atol`. The value corresponding to each key is the respective
infeasibilility. Infeasibility is defined as the distance between the primal
value of the constraint (see `MOI.ConstraintPrimal`) and the nearest point by
Euclidean distance in the corresponding set.

## Notes

 * If `skip_missing = true`, constraints containing variables that are not in
   `point` will be ignored.
 * If `skip_missing = false` and a partial primal solution is provided, an error
   will be thrown.
 * If no point is provided, the primal solution from the last time the model was
   solved is used.

## Examples

```jldoctest; setup=:(using JuMP)
julia> model = Model();

julia> @variable(model, 0.5 <= x <= 1);

julia> primal_feasibility_report(model, Dict(x => 0.2))
Dict{Any,Float64} with 1 entry:
  x ≥ 0.5 => 0.3
```
"""
function primal_feasibility_report(
    model::Model,
    point::AbstractDict{VariableRef,Float64} = _last_primal_solution(model);
    atol::Float64 = 0.0,
    skip_missing::Bool = false,
)
    function point_f(x::VariableRef)
        fx = get(point, x, missing)
        if ismissing(fx) && !skip_missing
            error(
                "point does not contain a value for variable $x. Provide a " *
                "value, or pass `skip_missing = true`.",
            )
        end
        return fx
    end
    violated_constraints = Dict{Any,Float64}()
    for (F, S) in list_of_constraint_types(model)
        _add_infeasible_constraints(
            model,
            F,
            S,
            violated_constraints,
            point_f,
            atol,
        )
    end
    if num_nl_constraints(model) > 0
        if skip_missing
            error(
                "`skip_missing = true` is not allowed when nonlinear " *
                "constraints are present.",
            )
        end
        _add_infeasible_nonlinear_constraints(
            model,
            violated_constraints,
            point_f,
            atol,
        )
    end
    return violated_constraints
end

function _add_infeasible_constraints(
    model::Model,
    ::Type{F},
    ::Type{S},
    violated_constraints::Dict{Any,Float64},
    point_f::Function,
    atol::Float64,
) where {F,S}
    for con in all_constraints(model, F, S)
        obj = constraint_object(con)
        d = _distance_to_set(value.(obj.func, point_f), obj.set)
        if d > atol
            violated_constraints[con] = d
        end
    end
    return
end

function _add_infeasible_nonlinear_constraints(
    model::Model,
    violated_constraints::Dict{Any,Float64},
    point_f::Function,
    atol::Float64,
)
    evaluator = NLPEvaluator(model)
    MOI.initialize(evaluator, Symbol[])
    g = zeros(num_nl_constraints(model))
    MOI.eval_constraint(evaluator, g, point_f.(all_variables(model)))
    for (i, con) in enumerate(model.nlp_data.nlconstr)
        d = max(0.0, con.lb - g[i], g[i] - con.ub)
        if d > atol
            cref =
                ConstraintRef(model, NonlinearConstraintIndex(i), ScalarShape())
            violated_constraints[cref] = d
        end
    end
    return
end

function _distance_to_set(::Any, set::MOI.AbstractSet)
    return error(
        "Feasibility checker for set type $(typeof(set)) has not been " *
        "implemented yet.",
    )
end

_distance_to_set(::Missing, ::MOI.AbstractSet) = 0.0

###
### MOI.AbstractScalarSets
###

function _distance_to_set(x::T, set::MOI.LessThan{T}) where {T<:Real}
    return max(x - set.upper, zero(T))
end

function _distance_to_set(x::T, set::MOI.GreaterThan{T}) where {T<:Real}
    return max(set.lower - x, zero(T))
end

function _distance_to_set(x::T, set::MOI.EqualTo{T}) where {T<:Number}
    return abs(set.value - x)
end

function _distance_to_set(x::T, set::MOI.Interval{T}) where {T<:Real}
    return max(x - set.upper, set.lower - x, zero(T))
end

function _distance_to_set(x::T, ::MOI.ZeroOne) where {T<:Real}
    return min(abs(x - zero(T)), abs(x - one(T)))
end

function _distance_to_set(x::T, ::MOI.Integer) where {T<:Real}
    return abs(x - round(x))
end

function _distance_to_set(x::T, set::MOI.Semicontinuous{T}) where {T<:Real}
    return min(max(x - set.upper, set.lower - x, zero(T)), abs(x))
end

function _distance_to_set(x::T, set::MOI.Semiinteger{T}) where {T<:Real}
    d = max(ceil(set.lower) - x, x - floor(set.upper), abs(x - round(x)))
    return min(d, abs(x))
end

###
### MOI.AbstractVectorSets
###

function _check_dimension(v::AbstractVector, s)
    if length(v) != MOI.dimension(s)
        throw(DimensionMismatch("Mismatch between value and set"))
    end
    return
end

function _distance_to_set(x::Vector{T}, set::MOI.Nonnegatives) where {T<:Real}
    _check_dimension(x, set)
    return LinearAlgebra.norm(max(-xi, zero(T)) for xi in x)
end

function _distance_to_set(x::Vector{T}, set::MOI.Nonpositives) where {T<:Real}
    _check_dimension(x, set)
    return LinearAlgebra.norm(max(xi, zero(T)) for xi in x)
end

function _distance_to_set(x::Vector{T}, set::MOI.Zeros) where {T<:Number}
    _check_dimension(x, set)
    return LinearAlgebra.norm(x)
end

function _distance_to_set(x::Vector{T}, set::MOI.Reals) where {T<:Real}
    _check_dimension(x, set)
    return zero(T)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import Pkg; Pkg.activate("../..")
using StatsPlots
import Random

Random.seed!(0)
X = cumsum(randn(10, 4) .+ 1, dims=1) * .1 .+ [1 2 3 4]
plot(X, legendfontsize=10, m=:square, ms=1, msc=:white, legend=:outerright,
     lw=5, palette=:tab10, msw=0)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
import ASTInterpreter2
import DebuggerFramework: DebuggerState, execute_command, print_status, locinfo, eval_code, dummy_state
import ..Atom: fullpath, handle, @msg, wsitem, Inline, EvalError, Console
import Juno: Row, ProgressBar, Tree
import REPL
using Media
using MacroTools

chan = nothing
state = nothing

isdebugging() = chan ≠ nothing

# entrypoint
function enter(arg)
  quote
    # FIXME: should use `qualifyASTInterpreter(ASTInterpreter2._make_stack(arg))`
    # instead, but I can't get that to work.
    let stack = $(make_stack(arg))
      Atom.Debugger.startdebugging(stack)
    end
  end
end

function make_stack(arg)
    arg = expand(arg)
    @assert isa(arg, Expr) && arg.head == :call
    kws = collect(filter(x->isexpr(x,:kw),arg.args))
    if !isempty(kws)
      args = Expr(:tuple,:(Core.kwfunc($(args[1]))),
        Expr(:call,Base.vector_any,mapreduce(
          x->[QuoteNode(x.args[1]),x.args[2]],vcat,kws)...),
        map(x->isexpr(x,:parameters)?QuoteNode(x):x,
          filter(x->!isexpr(x,:kw),arg.args))...)
    else
      args = Expr(:tuple,
        map(x->isexpr(x,:parameters)?QuoteNode(x):x, arg.args)...)
    end
    quote
        theargs = $(esc(args))
        stack = [Atom.Debugger.ASTInterpreter2.enter_call_expr(Expr(:call,theargs...))]
        Atom.Debugger.ASTInterpreter2.maybe_step_through_wrapper!(stack)
        stack[1] = Atom.Debugger.ASTInterpreter2.JuliaStackFrame(stack[1], Atom.Debugger.ASTInterpreter2.maybe_next_call!(stack[1]))
        stack
    end
end

function qualifyASTInterpreter(expr)
  postwalk(expr) do ex
    if @capture(ex, f_(xs__)) && startswith(string(f), "ASTInterpreter2")
      :($(Symbol("Atom.Debugger.", f))($(xs...)))
    else
      ex
    end
  end
end

# setup interpreter
function startdebugging(stack)
  global state = dummy_state(stack)
  global chan = Channel(0)

  debugmode(true)

  if Atom.isREPL()
    # FIXME: Should get rid of the second code path (see comment in repl.jl).
    if Atom.repleval[]
      t = @async debugprompt()
    else
      Atom.changeREPLprompt("debug> ", color=:magenta)
    end
  end

  stepto(state)
  res = nothing

  try
    evalscope() do
      for val in chan
        if val == :nextline
          execute_command(state, state.stack[state.level], Val{:n}(), "n")
        elseif val == :stepexpr
            execute_command(state, state.stack[state.level], Val{:nc}(), "nc")
        elseif val == :stepin
          execute_command(state, state.stack[state.level], Val{:s}(), "s")
        elseif val == :finish
          execute_command(state, state.stack[state.level], Val{:finish}(), "finish")
        elseif typeof(val) <: Tuple && val[1] == :toline
          while locinfo(state.stack[state.level]).line < val[2]
            execute_command(state, state.stack[state.level], Val{:n}(), "n")
            length(state.stack) == 0 && break
          end
        else
          warn("Internal: Unknown debugger message $val.")
        end
        length(state.stack) == 0 && break
        stepto(state)
      end
    end
  catch e
    ee = EvalError(e, catch_stacktrace())
    if Atom.isREPL()
      print("\r                        \r")
      print_with_color(:red, stderr, "ERROR: ")
      Base.showerror(stderr, e, backtrace())
      println(stderr)
    else
      render(Console(), ee)
    end
  finally
    res = state.overall_result
    chan = nothing
    state = nothing
    debugmode(false)

    if Atom.isREPL()
      if Atom.repleval[]
        istaskdone(t) || schedule(t, InterruptException(); error=true)
        print("\r                        \r")
      else
        Atom.changeREPLprompt(Atom.current_prompt, color = :green)
      end
    end
  end
  res
end

function debugprompt()
  try
    panel = Base.LineEdit.Prompt("debug> ";
              prompt_prefix="\e[35m",
              prompt_suffix=Base.text_colors[:white],
              on_enter = s -> true)

    panel.on_done = (s, buf, ok) -> begin
      Atom.msg("working")

      line = String(take!(buf))

      isempty(line) && return true

      if !ok
        Base.LineEdit.transition(s, :abort)
        Base.LineEdit.reset_state(s)
        return false
      end

      try
        r = Atom.Debugger.interpret(line)
        r ≠ nothing && display(r)
        println()
      catch e
        print_with_color(:red, stderr, "ERROR: ")
        Base.showerror(stderr, e, backtrace())
        println(stderr)
      end

      Atom.msg("doneWorking")
      Atom.msg("updateWorkspace")

      return true
    end

    REPL.run_interface(Base.active_repl.t, Base.LineEdit.ModalInterface([panel]))
  catch e
    e isa InterruptException || rethrow(e) end
end

# setup handlers for stepper commands
for cmd in [:nextline, :stepin, :stepexpr, :finish]
  handle(()->put!(chan, cmd), string(cmd))
end

handle((line)->put!(chan, (:toline, line)), "toline")

# notify the frontend that we start debugging now
debugmode(on) = @msg debugmode(on)

## Stepping

# perform step in frontend
function stepto(state::DebuggerState)
  frame = state.stack[state.level]
  loc = locinfo(frame)
  if loc == nothing
    file, line = ASTInterpreter2.determine_line_and_file(frame, frame.pc.next_stmt)[end]
  else
    file, line = loc.filepath, loc.line
  end
  stepto(Atom.fullpath(string(file)), line, stepview(nextstate(state)))
end
stepto(file, line, text) = @msg stepto(file, line, text)
stepto(::Nothing) = debugmode(false)

# return expression that will be evaluated next
function nextstate(state)
  frame = state.stack[state.level]
  expr = ASTInterpreter2.pc_expr(frame, frame.pc)
  isa(expr, Expr) && (expr = copy(expr))
  if isexpr(expr, :(=))
      expr = expr.args[2]
  end
  if isexpr(expr, :call) || isexpr(expr, :return)
      expr.args = map(var -> ASTInterpreter2.lookup_var_if_var(frame, var), expr.args)
  end
  expr
end

function stepview(ex)
  out = if @capture(ex, f_(as__))
    Row(span(".syntax--support.syntax--function", string(typeof(f).name.mt.name)),
             text"(", interpose(as, text", ")..., text")")
  elseif @capture(ex, x_ = y_)
    Row(Text(string(x)), text" = ", y)
  elseif @capture(ex, return x_)
    Row(span(".syntax--support.syntax--keyword", "return "), x)
  else
    Text(string(ex))
  end
  render(Inline(), out)
end

function evalscope(f)
  try
    @msg doneWorking()
    unlock(Atom.evallock)
    f()
  finally
    lock(Atom.evallock)
    @msg working()
  end
end

## Workspace

function contexts(s::DebuggerState = state)
  cxt = []
  trace = ""
  for frame in reverse(state.stack)
    trace = string(trace, "/", frame.meth.name)
    c = Dict(:context => string("Debug: ", trace), :items => localvars(frame))
    push!(cxt, c)
  end
  cxt
end

function localvars(frame::ASTInterpreter2.JuliaStackFrame)
  items = []
  for i = 1:length(frame.locals)
    if !isnull(frame.locals[i])
      # #self# is only interesting if it has values inside of it. We already know
      # which function we're in otherwise.
      val = get(frame.locals[i])
      if frame.code.slotnames[i] == Symbol("#self#") && (isa(val, Type) || sizeof(val) == 0)
        continue
      end
      push!(items, wsitem(frame.code.slotnames[i], val))
    end
  end
  for i = 1:length(frame.sparams)
    push!(items, wsitem(frame.meth.sparam_syms[i], get(Nullable{Any}(frame.sparams[i]), Undefined())))
  end
  return items
end

type Undefined end
@render Inline u::Undefined span(".fade", "<undefined>")
Atom.wsicon(::Undefined) = "icon-circle-slash"

## Evaluation

function interpret(code::AbstractString, frame::DebuggerState = state)
  eval_code(state, state.stack[state.level], code)
end


text:
@get_units @with_kw struct SurfaceFlow{T,R,L}
    β::T | "-"                              # constant in Manning's equation
    sl::Vector{T} | "m m-1"                 # Slope [m m⁻¹]
    n::Vector{T} | "s m-1/3"                # Manning's roughness [s m⁻⅓]
    dl::Vector{T} | "m"                     # Drain length [m]
    q::Vector{T} | "m3 s-1"                 # Discharge [m³ s⁻¹]
    qin::Vector{T} | "m3 s-1"               # Inflow from upstream cells [m³ s⁻¹]
    q_av::Vector{T} | "m3 s-1"              # Average discharge [m³ s⁻¹]
    qlat::Vector{T} | "m2 s-1"              # Lateral inflow per unit length [m² s⁻¹]
    inwater::Vector{T} | "m3 s-1"           # Lateral inflow [m³ s⁻¹]
    inflow::Vector{T} | "m3 s-1"            # External inflow (abstraction/supply/demand) [m³ s⁻¹]
    volume::Vector{T} | "m3"                # Kinematic wave volume [m³] (based on water level h)
    h::Vector{T} | "m"                      # Water level [m]
    h_av::Vector{T} | "m"                   # Average water level [m]
    h_bankfull::Vector{T} | "m"             # Bankfull water level [m]
    Δt::T | "s"                             # Model time step [s]
    its::Int | "-"                          # Number of fixed iterations
    width::Vector{T} | "m"                  # Flow width [m]
    alpha_pow::T | "-"                      # Used in the power part of α
    alpha_term::Vector{T} | "-"             # Term used in computation of α
    α::Vector{T} | "s3/5 m1/5"              # Constant in momentum equation A = αQᵝ, based on Manning's equation
    cel::Vector{T} | "m s-1"                # Celerity of the kinematic wave
    to_river::Vector{T} | "m3 s-1"          # Part of overland flow [m³ s⁻¹] that flows to the river
    rivercells::Vector{Bool} | "-"          # Location of river cells (0 or 1)
    wb_pit::Vector{Bool} | "-"              # Boolean location (0 or 1) of a waterbody (wb, reservoir or lake).
    reservoir_index::Vector{Int} | "-"      # map cell to 0 (no reservoir) or i (pick reservoir i in reservoir field)
    lake_index::Vector{Int} | "-"           # map cell to 0 (no lake) or i (pick lake i in lake field)
    reservoir::R                            # Reservoir model struct of arrays
    lake::L                                 # Lake model struct of arrays
    kinwave_it::Bool                        # Boolean for iterations kinematic wave

    # TODO unclear why this causes a MethodError
    # function SurfaceFlow{T,R,L}(args...) where {T,R,L}
    #     equal_size_vectors(args)
    #     return new(args...)
    # end
end

statevars(::SurfaceFlow) = (:q, :h, :h_av)

function update(
    sf::SurfaceFlow,
    network;
    frac_toriver = nothing,
    inflow_wb = nothing,
    doy = 0,
)
    @unpack graph,
    order,
    subdomain_order,
    topo_subdomain,
    indices_subdomain,
    upstream_nodes = network

    n = length(order)
    ns = length(subdomain_order)

    @. sf.alpha_term = pow(sf.n / sqrt(sf.sl), sf.β)
    # use fixed alpha value based on 0.5 * h_bankfull
    @. sf.α = sf.alpha_term * pow(sf.width + sf.h_bankfull, sf.alpha_pow)

    # two options for iteration, fixed or based on courant number.
    if sf.kinwave_it
        if sf.its > 0
            its = sf.its
        else
            # calculate celerity
            courant = zeros(n)
            for k = 1:ns
                for m in subdomain_order[k]
                    for v in topo_subdomain[m]
                        if sf.q[v] > 0.0
                            sf.cel[v] = 1.0 / (sf.α[v] * sf.β * pow(sf.q[v], (sf.β - 1.0)))
                            courant[v] = (sf.Δt / sf.dl[v]) * sf.cel[v]
                        end
                    end
                end
            end
            filter!(x -> x ≠ 0.0, courant)
            its = isempty(courant) ? 1 : ceil(Int, (1.25 * quantile!(courant, 0.95)))
        end
    else
        its = 1
    end

    # sub time step
    adt = sf.Δt / its

    sf.q_av .= 0.0
    sf.h_av .= 0.0
    sf.to_river .= 0.0
    # because of possible iterations set reservoir and lake inflow and total outflow at
    # start to zero, the total sum of inflow and outflow at each sub time step is calculated
    if !isnothing(sf.reservoir)
        sf.reservoir.inflow .= 0.0
        sf.reservoir.totaloutflow .= 0.0
    end
    if !isnothing(sf.lake)
        sf.lake.inflow .= 0.0
        sf.lake.totaloutflow .= 0.0
    end

    for _ = 1:its
        sf.qin .= 0.0
        for k = 1:ns
            @threads for m in subdomain_order[k]
                for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])

                    # for overland flow frac_toriver needs to be defined
                    if frac_toriver !== nothing # run kinematic wave for land domain
                        # for a river cell without a reservoir or lake (wb_pit is false) part of the
                        # upstream surface flow goes to the river (frac_toriver) and part goes to
                        # the surface flow reservoir (1.0 - frac_toriver), upstream nodes with a
                        # reservoir or lake are excluded
                        sf.to_river[v] += sum_at(
                            i -> sf.q[i] * frac_toriver[i],
                            upstream_nodes[n],
                            eltype(sf.to_river),
                        )
                        if sf.width[v] > 0.0
                            sf.qin[v] = sum_at(
                                i -> sf.q[i] * (1.0 - frac_toriver[i]),
                                upstream_nodes[n],
                                eltype(sf.q),
                            )
                        end
                    else # run kinematic wave for river domain (including reservoirs and lakes)
                        # sf.qin by outflow from upstream reservoir or lake location is added
                        sf.qin[v] += sum_at(sf.q, upstream_nodes[n])
                    end

                    # Inflow supply/abstraction is added to qlat (divide by flow length)
                    # If inflow < 0, abstraction is limited
                    if sf.inflow[v] < 0.0
                        max_abstract = min(sf.qin[v] + sf.qlat[v] * sf.dl[v], -sf.inflow[v])
                        inflow = -max_abstract / sf.dl[v]
                    else
                        inflow = sf.inflow[v] / sf.dl[v]
                    end

                    sf.q[v] = kinematic_wave(
                        sf.qin[v],
                        sf.q[v],
                        sf.qlat[v] + inflow,
                        sf.α[v],
                        sf.β,
                        adt,
                        sf.dl[v],
                    )

                    if !isnothing(sf.reservoir) && sf.reservoir_index[v] != 0
                        # run reservoir model and copy reservoir outflow to inflow (qin) of
                        # downstream river cell
                        i = sf.reservoir_index[v]
                        update(sf.reservoir, i, sf.q[v] + inflow_wb[v], adt)

                        downstream_nodes = outneighbors(graph, v)
                        n_downstream = length(downstream_nodes)
                        if n_downstream == 1
                            j = only(downstream_nodes)
                            sf.qin[j] = sf.reservoir.outflow[i]
                        elseif n_downstream == 0
                            error(
                                """A reservoir without a downstream river node is not supported.
                                Add a downstream river node or move the reservoir to an upstream node (model schematization).
                                """)
                        else
                            error("bifurcations not supported")
                        end

                    elseif !isnothing(sf.lake) && sf.lake_index[v] != 0
                        # run lake model and copy lake outflow to inflow (qin) of downstream river
                        # cell
                        i = sf.lake_index[v]
                        update(sf.lake, i, sf.q[v] + inflow_wb[v], doy, adt)

                        downstream_nodes = outneighbors(graph, v)
                        n_downstream = length(downstream_nodes)
                        if n_downstream == 1
                            j = only(downstream_nodes)
                            sf.qin[j] = sf.lake.outflow[i]
                        elseif n_downstream == 0
                            error(
                                """A lake without a downstream river node is not supported.
                                Add a downstream river node or move the lake to an upstream node (model schematization).
                                """)
                        else
                            error("bifurcations not supported")
                        end
                    end

                    # update h
                    crossarea = sf.α[v] * pow(sf.q[v], sf.β)
                    sf.h[v] = crossarea / sf.width[v]

                    sf.q_av[v] += sf.q[v]
                    sf.h_av[v] += sf.h[v]
                end

            end
        end
    end
    sf.q_av ./= its
    sf.h_av ./= its
    sf.to_river ./= its
    sf.volume .= sf.dl .* sf.width .* sf.h
end

@get_units @with_kw struct LateralSSF{T}
    kh₀::Vector{T} | "m Δt-1"              # Horizontal hydraulic conductivity at soil surface [m Δt⁻¹]
    f::Vector{T} | "m-1"                   # A scaling parameter [m⁻¹] (controls exponential decline of kh₀)
    soilthickness::Vector{T} | "m"         # Soil thickness [m]
    θₛ::Vector{T} | "-"                     # Saturated water content (porosity) [-]
    θᵣ::Vector{T} | "-"                    # Residual water content [-]
    t::T | "Δt s"                          # time step [Δt s]
    Δt::T | "s"                            # model time step [s]
    βₗ::Vector{T} | "m m-1"                 # Slope [m m⁻¹]
    dl::Vector{T} | "m"                    # Drain length [m]
    dw::Vector{T} | "m"                    # Flow width [m]
    zi::Vector{T} | "m"                    # Pseudo-water table depth [m] (top of the saturated zone)
    exfiltwater::Vector{T} | "m Δt-1"      # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions)
    recharge::Vector{T} | "m Δt-1"         # Net recharge to saturated store [m Δt⁻¹]
    ssf::Vector{T} | "m3 Δt-1"             # Subsurface flow [m³ Δt⁻¹]
    ssfin::Vector{T} | "m3 Δt-1"           # Inflow from upstream cells [m³ Δt⁻¹]
    ssfmax::Vector{T} | "m2 Δt-1"          # Maximum subsurface flow [m² Δt⁻¹]
    to_river::Vector{T} | "m3 Δt-1"        # Part of subsurface flow [m³ Δt⁻¹] that flows to the river
    wb_pit::Vector{Bool} | "-"             # Boolean location (0 or 1) of a waterbody (wb, reservoir or lake).

    function LateralSSF{T}(args...) where {T}
        equal_size_vectors(args)
        return new(args...)
    end
end

statevars(::LateralSSF) = (:ssf,)

function update(ssf::LateralSSF, network, frac_toriver)
    @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network

    ns = length(subdomain_order)
    for k = 1:ns
        @threads for m in subdomain_order[k]
            for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])
                # for a river cell without a reservoir or lake (wb_pit is false) part of the
                # upstream subsurface flow goes to the river (frac_toriver) and part goes to the
                # subsurface flow reservoir (1.0 - frac_toriver) upstream nodes with a reservoir or
                # lake are excluded
                ssf.ssfin[v] = sum_at(
                    i -> ssf.ssf[i] * (1.0 - frac_toriver[i]),
                    upstream_nodes[n],
                    eltype(ssf.ssfin),
                )
                ssf.to_river[v] = sum_at(
                    i -> ssf.ssf[i] * frac_toriver[i],
                    upstream_nodes[n],
                    eltype(ssf.to_river),
                )

                ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
                    ssf.ssfin[v],
                    ssf.ssf[v],
                    ssf.zi[v],
                    ssf.recharge[v],
                    ssf.kh₀[v],
                    ssf.βₗ[v],
                    ssf.θₛ[v] - ssf.θᵣ[v],
                    ssf.f[v],
                    ssf.soilthickness[v],
                    ssf.t,
                    ssf.dl[v],
                    ssf.dw[v],
                    ssf.ssfmax[v],
                )
            end
        end
    end
end

@get_units @with_kw struct GroundwaterExchange{T}
    Δt::T | "s"                         # model time step [s]
    exfiltwater::Vector{T} | "m Δt-1"   # Exfiltration [m Δt⁻¹]  (groundwater above surface level, saturated excess conditions)
    zi::Vector{T} | "m"                 # Pseudo-water table depth [m] (top of the saturated zone)
    to_river::Vector{T} | "m3 Δt-1"     # Part of subsurface flow [m³ Δt⁻¹] that flows to the river
    ssf::Vector{T} | "m3 Δt-1"          # Subsurface flow [m³ Δt⁻¹]
end


text:
#=
By starting at the top of the triangle below and moving to adjacent numbers on the row below, the maximum total from top to bottom is 23.

3
7 4
2 4 6
8 5 9 3

That is, 3 + 7 + 4 + 9 = 23.

Find the maximum total from top to bottom of the triangle below:

75
95 64
17 47 82
18 35 87 10
20 04 82 47 65
19 01 23 75 03 34
88 02 77 73 07 63 67
99 65 04 28 06 16 70 92
41 41 26 56 83 40 80 70 33
41 48 72 33 47 32 37 16 94 29
53 71 44 65 25 43 91 52 97 51 14
70 11 33 28 77 73 17 78 39 68 17 57
91 71 52 38 17 14 91 43 58 50 27 29 48
63 66 04 68 89 53 67 30 73 16 69 87 40 31
04 62 98 27 23 09 70 98 73 93 38 53 60 04 23

NOTE: As there are only 16384 routes, it is possible to solve this problem by trying every route. However, Problem 67, is the same challenge with a triangle containing one-hundred rows; it cannot be solved by brute force, and requires a clever method! ;o)
=#
a = Array[
[75],
[95, 64],
[17, 47, 82],
[18, 35, 87, 10],
[20, 04, 82, 47, 65],
[19, 01, 23, 75, 03, 34],
[88, 02, 77, 73, 07, 63, 67],
[99, 65, 04, 28, 06, 16, 70, 92],
[41, 41, 26, 56, 83, 40, 80, 70, 33],
[41, 48, 72, 33, 47, 32, 37, 16, 94, 29],
[53, 71, 44, 65, 25, 43, 91, 52, 97, 51, 14],
[70, 11, 33, 28, 77, 73, 17, 78, 39, 68, 17, 57],
[91, 71, 52, 38, 17, 14, 91, 43, 58, 50, 27, 29, 48],
[63, 66, 04, 68, 89, 53, 67, 30, 73, 16, 69, 87, 40, 31],
[04, 62, 98, 27, 23, 09, 70, 98, 73, 93, 38, 53, 60, 04, 23]]

function calc(arr)
  result = Array[]
  push!(result, arr[length(arr)])
  for row in reverse(1:length(arr)-1)
    nextrow = Int64[]
    for column in 1:length(arr[row])
      lowleft = last(result)[column]
      lowright = last(result)[column+1]
      center = arr[row][column]
      push!(nextrow, max(center+lowleft, center+lowright))
    end
    push!(result, nextrow)
  end
  last(result)[1]
end
@time println(calc(a))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
@testset "utils" begin
    @printline
    model2 = PCSAFT(["water","ethanol"])
    model4 = SAFTgammaMie(["methane","butane","isobutane","pentane"])
    gc3 = UNIFAC(["propane","butane","isobutane"])
    gc2 = SAFTgammaMie([
        "ethanol",
        ("ibuprofen", ["CH3"=>3, "COOH"=>1, "aCCH"=>1, "aCCH2"=>1, "aCH"=>4])])

    ideal1 = WalkerIdeal(["hexane"])
    noparam1 = gc3.puremodel[1].translation
    simple1 = gc3.puremodel[1].alpha
    @testset "split_model" begin
        models2 = split_model(model2)
        @test models2[1].components[1] == model2.components[1]
        @test models2[2].components[1] == model2.components[2]
        @test models2[1].icomponents == models2[2].icomponents == 1:1

        model2_unsplit = only(split_model(model2,[[1,2]]))
        @test model2_unsplit.icomponents == model2.icomponents
        @test model2_unsplit.components == model2.components

        model4_split = Clapeyron.split_model(model4)
        @test model4_split[4].groups.n_groups[1][2] == 3
        @test model4_split[1].components[1] == "methane"
        @test all(isone(length(model4_split[i].components)) for i in 1:4)

        gc3_split = Clapeyron.split_model(gc3)
        @test all(isone(length(gc3_split[i].components)) for i in 1:3)
        @test all(isone(length(gc3_split[i].puremodel)) for i in 1:3)
    end

    @testset "macros" begin
        comps(model) = Clapeyron.@comps
        groups(model) = Clapeyron.@groups
        groups(model,i) = Clapeyron.@groups(i)
        sites(model,i) = Clapeyron.@sites(i)
        f_eos(model,V,T,z,c) = 2 + c
        model = V = T = z = nothing

        @test groups(gc2) == 1:6
        @test comps(gc2) == 1:2
        @test groups(gc2,2) == [2,3,4,5,6]
        @test groups(gc3,1) == groups(gc3,2) #propane and butane has the same amount of groups
        @test sites(gc2.vrmodel,2) == [3,4,5]
        @test sites(model2,1) == [1,2]
        @test Clapeyron.@f(f_eos,pi) == 2+pi
        @test Clapeyron.@nan(Base.log(-1),3) == 3
        @test_throws MethodError Clapeyron.@nan(Base.log("s"),3)
    end
    using Clapeyron: has_sites,has_groups
    @testset "has_sites-has_groups" begin
        @test has_sites(typeof(gc3)) == false
        @test has_sites(typeof(model2)) == has_sites(typeof(model4)) == has_sites(typeof(gc2)) == true
        @test has_groups(typeof(gc2)) == has_groups(typeof(gc3)) == has_groups(typeof(ideal1)) == true
        @test has_groups(typeof(simple1)) == has_groups(typeof(model2)) == false
    end

    @testset "eosshow" begin
        #@newmodelgc
        @test repr(ideal1) == "WalkerIdeal{BasicIdeal}(\"hexane\")"
        @test repr("text/plain",ideal1) == "WalkerIdeal{BasicIdeal} with 1 component:\n \"hexane\": \"CH3\" => 2, \"CH2\" => 4\nContains parameters: Mw, Nrot, theta1, theta2, theta3, theta4, deg1, deg2, deg3, deg4"
        #@newmodel
        @test repr(model2) == "PCSAFT{BasicIdeal}(\"water\", \"ethanol\")"
        @test repr("text/plain",model2) == "PCSAFT{BasicIdeal} with 2 components:\n \"water\"\n \"ethanol\"\nContains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol"
        #@newmodelsimple
        @test repr(noparam1) == "NoTranslation()"
        @test repr("text/plain",noparam1) == "NoTranslation\n"
        @test repr(simple1) == "PRAlpha(\"propane\")"
        @test repr("text/plain",simple1) == "PRAlpha with 1 component:\n \"propane\"\nContains parameters: acentricfactor"
    end
    @printline
end


text:
#!/usr/bin/env julia

using Pkg
using Dates
using Distributed
using MIToS.Utils.Scripts

Args = parse_commandline(
    # TO DO ----------------------------------------------------------------------
    ["--format", "-f"],
    Dict(
        :help => "Format of the MSA: Stockholm, Raw or FASTA",
        :arg_type => String,
        :default => "Stockholm"
    ),
    ["--clustering", "-c"],
    Dict(
        :help => "Sequence clustering (Hobohm I)",
        :action => :store_false
    ),
    ["--threshold", "-i"],
    Dict(
        :help => "Percent identity threshold for sequence clustering (Hobohm I)",
        :arg_type => Float64,
        :default => 62.0
    ),
    # Keywords...
    description="""
    This takes a MSA file as input and it calculates and saves on *.conservation.csv the
    Shannon entropy (H) and Kullback-Leibler divergence (KL) values for each column
    (Johansson and Toh 2010).\n
    It is possible to do a sequence clustering using the Hobohm I algorithm to avoid the
    effect of sequence redundancy in the conservation scores. Each sequence in a cluster
    is weighted using the inverse of the number of elements in that cluster.\n
    Shannon entropy is a common measure of the residue variability of a particular MSA column.
    For each column, we consider the frequency of the 20 natural protein residues. This uses
    the Euler's number e as the base of the logarithm, so the entropy is measured in nats.\n
    The Kullback-Leibler divergence, also called relative entropy, is a measure of residue
    conservation. It measures how much a probability distribution differs from a background
    distribution. In particular, this implementation measures the divergence between the
    residue distribution of an MSA column and the probabilities derived from the BLOSUM62
    substitution matrix.\n
    Johansson, F., Toh, H., 2010.
    A comparative study of conservation and variation scores.
    BMC Bioinformatics 11, 388.
    """,
    output=".conservation.csv"
    # ----------------------------------------------------------------------------
    )

set_parallel(Args["parallel"])

@everywhere begin

    const args = remotecall_fetch(()->Args,1)

    import MIToS.Utils.Scripts: script

    # TO DO ----------------------------------------------------------------------
    using MIToS.MSA
    using MIToS.Information
    using PairwiseListMatrices
    # ----------------------------------------------------------------------------

    function script(input::Union{Base.LibuvStream,  AbstractString},
                    args,
                    fh_out::Union{Base.LibuvStream, IO})
        # TO DO ------------------------------------------------------------------
        println(fh_out, "# MIToS ", Pkg.installed()["MIToS"], " Conservation.jl ", now())
        println(fh_out, "# used arguments:")
        for (key, value) in args
            println(fh_out, "# \t", key, "\t\t", value)
        end
        form = string(args["format"])
        if form == "Stockholm"
            msa = readorparse(input, Stockholm)
        elseif form == "FASTA"
            msa = readorparse(input, FASTA)
        elseif form == "Raw"
            msa = readorparse(input, Raw)
        else
            throw(ErrorException("--format should be Stockholm, Raw or FASTA."))
        end
        println(fh_out, "# Number of sequences: ", nsequences(msa))
        if args["clustering"]
            clusters = hobohmI(msa, args["threshold"])
            println(fh_out, "# Number of sequence clusters: ", nclusters(clusters))
        else
            clusters = NoClustering()
        end
        probability_table = Probabilities(ContingencyTable(Float64, Val{1}, UngappedAlphabet()))

        KL = mapcolfreq!(kullback_leibler, msa, probability_table, weights = clusters)
        H  = mapcolfreq!(entropy, msa, probability_table, weights = clusters)

        col_names = names(KL,2)
        println(fh_out, "i,H,KL")
        for i in 1:length(col_names)
            println(fh_out, col_names[i], ",", H[i], ",", KL[i])
        end
        # ------------------------------------------------------------------------
    end

end

runscript(args)