text:
# TODO:
# - function definition (wrap the function like a pycallback)
# - class definition (e.g. `struct User <: BaseModel; id::int=0; name::str=""; end`)
# - property syntax (e.g. `classmethod |> function foo(cls); cls(); end`)
# - with syntax (`@with`)
# - bytes literals (`b"..."`)
# - raw string literals (`r"..."` or raw"..."?)
# - regex literals (`re"..."` or `r"..."`?)
# - formatted string literals (`f"..."`)
# - comprehensions (tuples, lists, dicts, sets, generators)
# - generator functions??
# - splatting

const PY_MACRO_UNOPS = Dict(
    # operators
    :(+) => (pypos, true),
    :(-) => (pyneg, true),
    :(~) => (pyinv, true),
    # builtins
    :abs => (pyabs, true),
    :ascii => (pyascii, true),
    :bytes => (pybytes, true),
    :bool => (pybool, true),
    :complex => (pycomplex, true),
    :dict => (pydict, true),
    :dir => (pydir, true),
    :float => (pyfloat, true),
    :frozenset => (pyfrozenset, true),
    :int => (pyint, true),
    :iter => (pyiter, true),
    :list => (pylist, true),
    :range => (pyrange, true),
    :repr => (pyrepr, true),
    :set => (pyset, true),
    :slice => (pyslice, true),
    :str => (pystr, true),
    :tuple => (pytuple, true),
    :type => (pytype, true),
    # builtins converting to julia
    :jlascii => (x->pyascii(String,x), false),
    :jlbool => (pytruth, false),
    :jlbytes => (x->pybytes(Base.CodeUnits,x), false),
    :jlhash => (pyhash, false),
    :jllen => (pylen, false),
    :jlrepr => (x->pyrepr(String,x), false),
    :jlstr => (x->pystr(String,x), false),
    # jlcomplex
)

const PY_MACRO_BINOPS = Dict(
    # operators
    :(+) => (pyadd, true),
    :(-) => (pysub, true),
    :(*) => (pymul, true),
    :(/) => (pytruediv, true),
    :(÷) => (pyfloordiv, true),
    :(%) => (pymod, true),
    :(^) => (pypow, true),
    :(<<) => (pylshift, true),
    :(>>) => (pyrshift, true),
    :(&) => (pyand, true),
    :(|) => (pyor, true),
    :(⊻) => (pyxor, true),
    :(==) => (pyeq, true),
    :(!=) => (pyne, true),
    :(≠ ) => (pyne, true),
    :(<=) => (pyle, true),
    :(≤ ) => (pyle, true),
    :(< ) => (pylt, true),
    :(>=) => (pyge, true),
    :(≥ ) => (pyge, true),
    :(> ) => (pygt, true),
    :(===) => (pyis, false),
    :(≡) => (pyis, false),
    :(!==) => (pyisnot, false),
    :(≢) => (pyisnot, false),
    :in => (pyin, false),
    :∈ => (pyin, false),
    :∉ => (pynotin, false),
    # builtins
    :pow => (pypow, true),
    :complex => (pycomplex, true),
    :delattr => (pydelattr, false),
    :divmod => (pydivmod, true),
    :getattr => (pygetattr, true),
    :hasattr => (pyhasattr, false),
    :issubclass => (pyissubclass, false),
    :isinstance => (pyisinstance, false),
    :range => (pyrange, true),
    :slice => (pyslice, true),
    # builtins converting to julia
    # jlcomplex
)

const PY_MACRO_TERNOPS = Dict(
    :pow => (pypow, true),
    :range => (pyrange, true),
    :setattr => (pysetattr, false),
    :slice => (pyslice, true),
)

Base.@kwdef mutable struct PyMacroState
    mod :: Module
    src :: LineNumberNode
    consts :: IdDict{Any,Py} = IdDict{Any,Py}()
    inits :: Vector{Any} = []
end

function py_macro_err(st, ex, msg=nothing)
    ex2 = ex isa Expr ? Expr(ex.head, [arg isa Expr ? :__ : arg for arg in ex.args]...) : ex
    message = "@py syntax error: $ex2"
    if msg !== nothing
        message *= ": " * msg
    end
    message *= "\n  at $(st.src.file):$(st.src.line)"
    error(message)
end

py_macro_assign(body, ans, ex) = push!(body, :($ans = $ex))

py_macro_del(body, var, tmp) = if tmp; push!(body, :($pydel!($var))); end

function py_macro_lower(st, body, ans, ex; flavour=:expr)

    # scalar literals
    if ex isa Union{Nothing, String, Bool, Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128, BigInt, Float16, Float32, Float64}
        x = get!(pynew, st.consts, ex)
        py_macro_assign(body, ans, x)
        return false

    # Int128 literals
    elseif isexpr(ex, :macrocall) && ex.args[1] === GlobalRef(Core, Symbol("@int128_str"))
        value = parse(Int128, ex.args[3])
        x = get!(pynew, st.consts, value)
        py_macro_assign(body, ans, x)
        return false

    # UInt128 literals
    elseif isexpr(ex, :macrocall) && ex.args[1] === GlobalRef(Core, Symbol("@uint128_str"))
        value = parse(UInt128, ex.args[3])
        x = get!(pynew, st.consts, value)
        py_macro_assign(body, ans, x)
        return false

    # big integer literals
    elseif isexpr(ex, :macrocall) && ex.args[1] === GlobalRef(Core, Symbol("@big_str"))
        value = parse(BigInt, ex.args[3])
        x = get!(pynew, st.consts, value)
        py_macro_assign(body, ans, x)
        return false

    # __file__
    elseif ex === :__file__
        return py_macro_lower(st, body, ans, string(st.src.file))

    # __line__
    elseif ex === :__line__
        return py_macro_lower(st, body, ans, st.src.line)

    # x
    elseif ex isa Symbol
        if ex in BUILTINS
            py_macro_assign(body, ans, :($pybuiltins.$ex))
        else
            py_macro_assign(body, ans, ex)
        end
        return false

    # x:y:z
    elseif flavour==:index && @capture(ex, ax_:ay_:az_)
        @gensym x y z
        tx = py_macro_lower(st, body, x, ax===:_ ? :None : ax)
        ty = py_macro_lower(st, body, y, ay===:_ ? :None : ay)
        tz = py_macro_lower(st, body, z, az===:_ ? :None : az)
        py_macro_assign(body, ans, :($pyslice($x, $y, $z)))
        py_macro_del(body, x, tx)
        py_macro_del(body, y, ty)
        py_macro_del(body, z, tz)
        return true

    # x:y
    elseif flavour==:index && @capture(ex, ax_:ay_)
        @gensym x y
        tx = py_macro_lower(st, body, x, ax===:_ ? :None : ax)
        ty = py_macro_lower(st, body, y, ay===:_ ? :None : ay)
        py_macro_assign(body, ans, :($pyslice($x, $y)))
        py_macro_del(body, x, tx)
        py_macro_del(body, y, ty)
        return true

    # x + y + z + ...
    elseif @capture(ex, +(ax_, ay_, az_, args__))
        return py_macro_lower(st, body, ans, foldl((x, y)->:($x+$y), (ax, ay, az, args...)))

    # x * y * z * ...
    elseif @capture(ex, *(ax_, ay_, az_, args__))
        return py_macro_lower(st, body, ans, foldl((x, y)->:($x*$y), (ax, ay, az, args...)))

    # f(args...; kwargs...)
    elseif isexpr(ex, :call)
        af = ex.args[1]
        # is it a special operator?
        isop = false
        if haskey(PY_MACRO_UNOPS, af)
            isop = true
            op, t = PY_MACRO_UNOPS[af]
            if length(ex.args) == 2
                ax, = ex.args[2:end]
                @gensym x
                tx = py_macro_lower(st, body, x, ax)
                py_macro_assign(body, ans, :($op($x)))
                py_macro_del(body, x, tx)
                return t
            end
        end
        if haskey(PY_MACRO_BINOPS, af)
            isop = true
            op, t = PY_MACRO_BINOPS[af]
            if length(ex.args) == 3
                ax, ay = ex.args[2:end]
                @gensym x y
                tx = py_macro_lower(st, body, x, ax)
                ty = py_macro_lower(st, body, y, ay)
                py_macro_assign(body, ans, :($op($x, $y)))
                py_macro_del(body, x, tx)
                py_macro_del(body, y, ty)
                return t
            end
        end
        if haskey(PY_MACRO_TERNOPS, af)
            isop = true
            op, t = PY_MACRO_TERNOPS[af]
            if length(ex.args) == 4
                ax, ay, az = ex.args[2:end]
                @gensym x y z
                tx = py_macro_lower(st, body, x, ax)
                ty = py_macro_lower(st, body, y, ay)
                tz = py_macro_lower(st, body, z, az)
                py_macro_assign(body, ans, :($op($x, $y, $z)))
                py_macro_del(body, x, tx)
                py_macro_del(body, y, ty)
                py_macro_del(body, z, tz)
                return t
            end
        end
        # if it's an operator but we get this far, we used the operator wrong
        # if it's also a builtin, we allow calling it as a function
        # otherwise it's a syntax error
        isop && af ∉ BUILTINS && py_macro_err(st, ex)
        # disallow things like calling other julia infix operators like functions
        af isa Symbol && !Base.isidentifier(af) && py_macro_err(st, ex)
        # pycall
        aargs = []
        akwargs = []
        aparams = []
        for a in ex.args[2:end]
            if isexpr(a, :parameters)
                append!(aparams, a.args)
            elseif isexpr(a, :kw)
                push!(akwargs, a)
            else
                push!(aargs, a)
            end
        end
        append!(akwargs, aparams)
        if !isempty(akwargs)
            @gensym f args kwargs
            tf = py_macro_lower(st, body, f, af)
            targs = py_macro_lower(st, body, args, Expr(:tuple, aargs...))
            tkwargs = py_macro_lower(st, body, kwargs, Expr(:braces, Expr(:parameters), akwargs...))
            py_macro_assign(body, ans, :($pycallargs($f, $args, $kwargs)))
            py_macro_del(body, f, tf)
            py_macro_del(body, args, targs)
            py_macro_del(body, kwargs, tkwargs)
        elseif !isempty(aargs)
            @gensym f args
            tf = py_macro_lower(st, body, f, af)
            targs = py_macro_lower(st, body, args, Expr(:tuple, aargs...))
            py_macro_assign(body, ans, :($pycallargs($f, $args)))
            py_macro_del(body, f, tf)
            py_macro_del(body, args, targs)
        else
            @gensym f
            tf = py_macro_lower(st, body, f, af)
            py_macro_assign(body, ans, :($pycallargs($f)))
            py_macro_del(body, f, tf)
        end
        return true

    # (...)
    elseif isexpr(ex, :tuple)
        if any(isexpr(arg, :...) for arg in ex.args)
            py_macro_err(st, ex, "splatting into tuples not implemented")
        else
            py_macro_assign(body, ans, :($pynulltuple($(length(ex.args)))))
            @gensym a
            for (i, aa) in enumerate(ex.args)
                ta = py_macro_lower(st, body, a, aa, flavour = flavour==:index ? :index : :expr)
                push!(body, :($pytuple_setitem($ans, $(i-1), $a)))
                py_macro_del(body, a, ta)
            end
            return true
        end

    # [...]
    elseif isexpr(ex, :vect)
        if any(isexpr(arg, :...) for arg in ex.args)
            py_macro_err(st, ex, "splatting into tuples not implemented")
        else
            py_macro_assign(body, ans, :($pynulllist($(length(ex.args)))))
            @gensym a
            for (i, aa) in enumerate(ex.args)
                ta = py_macro_lower(st, body, a, aa)
                push!(body, :($pylist_setitem($ans, $(i-1), $a)))
                py_macro_del(body, a, ta)
            end
            return true
        end

    # {...}
    elseif isexpr(ex, :braces)
        # Like Python, we allow braces to be set or dict literals.
        #
        # The expression is a dict if either:
        # - it is empty ({})
        # - it has a parameters section (e.g. {x,y;})
        # - it has any assignment expressions (e.g. {k=v})
        # - it has any pair expressions (e.g. {k=>v})
        # and otherwise it is a set.
        #
        # e.g. these are sets: {x, y}, {x...}
        # e.g. these are dicts: {}, {x, y;}, {k=x}, {k=>x}
        if isempty(ex.args) || any(isexpr(arg, :parameters, :(=), :kw) || @capture(arg, k_:v_) for arg in ex.args)
            aargs = []
            for a in ex.args
                if isexpr(a, :parameters)
                    append!(aargs, a.args)
                else
                    push!(aargs, a)
                end
            end
            py_macro_assign(body, ans, :($pydict()))
            @gensym k v
            for aa in aargs
                if isexpr(aa, :...)
                    py_macro_err(st, aa, "splatting into dicts not implemented")
                elseif isexpr(aa, :kw, :(=))
                    ak, av = aa.args
                    ak isa Symbol || py_macro_err(st, aa, "key of `key=value` must be a symbol - did you mean `key:value`?")
                    tk = py_macro_lower(st, body, k, string(ak))
                    tv = py_macro_lower(st, body, v, av)
                    push!(body, :($pydict_setitem($ans, $k, $v)))
                    py_macro_del(body, k, tk)
                    py_macro_del(body, v, tv)
                elseif @capture(aa, ak_ : av_)
                    tk = py_macro_lower(st, body, k, ak)
                    tv = py_macro_lower(st, body, v, av)
                    push!(body, :($pydict_setitem($ans, $k, $v)))
                    py_macro_del(body, k, tk)
                    py_macro_del(body, v, tv)
                else
                    py_macro_err(st, aa, "this kind of dict entry not implemented")
                end
            end
            return true
        else
            py_macro_assign(body, ans, :($pyset()))
            @gensym a
            for aa in ex.args
                if isexpr(aa, :...)
                    py_macro_err(st, ex, "splatting into sets not implemented")
                end
                ta = py_macro_lower(st, body, a, aa)
                push!(body, :($pyset_add($ans, $a)))
                py_macro_del(body, a, ta)
            end
            return true
        end

    # x.k
    elseif isexpr(ex, :.)
        ax, ak = ex.args
        @gensym x k
        tx = py_macro_lower(st, body, x, ax)
        if ak isa QuoteNode && ak.value isa Symbol
            tk = py_macro_lower(st, body, k, string(ak.value))
        else
            tk = py_macro_lower(st, body, k, ak)
        end
        py_macro_assign(body, ans, :($pygetattr($x, $k)))
        py_macro_del(body, x, tx)
        py_macro_del(body, k, tk)
        return true

    # x[k]
    elseif @capture(ex, ax_[ak__])
        @gensym x k
        tx = py_macro_lower(st, body, x, ax)
        if length(ak) == 1
            tk = py_macro_lower(st, body, k, ak[1]; flavour=:index)
        else
            tk = py_macro_lower(st, body, k, Expr(:tuple, ak...); flavour=:index)
        end
        py_macro_assign(body, ans, :($pygetitem($x, $k)))
        py_macro_del(body, x, tx)
        py_macro_del(body, k, tk)
        return true

    # x = y
    elseif @capture(ex, ax_ = ay_)
        ty = py_macro_lower(st, body, ans, ay)
        py_macro_lower_assign(st, body, ax, ans)
        return ty

    # @del x, y, ...
    elseif @capture(ex, @del (args__,))

        for arg in args

            # @del x
            if arg isa Symbol
                if arg in BUILTINS
                    py_macro_err(st, ex, "can't delete a builtin")
                else
                    push!(body, :($pydel!($arg::$Py)))
                end

            # @del x.k
            elseif isexpr(arg, :.)
                ax, ak = ex.args
                @gensym x k
                tx = py_macro_lower(st, body, x, ax)
                if ak isa QuoteNode && ak.value isa Symbol
                    tk = py_macro_lower(st, body, k, string(ak.value))
                else
                    tk = py_macro_lower(st, body, k, ak)
                end
                push!(body, :($pydelattr($x, $k)))
                py_macro_del(body, x, tx)
                py_macro_del(body, k, tk)

            # @del x[k]
            elseif @capture(arg, ax_[ak__])
                @gensym x k
                tx = py_macro_lower(st, body, x, ax)
                if length(ak) == 1
                    tk = py_macro_lower(st, body, k, ak[1], flavour=:index)
                else
                    tk = py_macro_lower(st, body, k, Expr(:tuple, ak...), flavour=:index)
                end
                push!(body, :($pydelitem($x, $k)))
                py_macro_del(body, x, tx)
                py_macro_del(body, k, tk)

            else
                py_macro_err(st, ex, "@del argument must be a variable, reference or property")
            end
        end
        py_macro_assign(body, ans, nothing)
        return false

    # @del x
    elseif @capture(ex, @del arg_)
        return py_macro_lower(st, body, ans, :(@del ($arg,)))

    # @jl x
    elseif @capture(ex, @jl ax_)
        y = py_macro_lower_jl(st, ax)
        py_macro_assign(body, ans, y)
        return false

    # @compile code mode=mode ...
    elseif @capture(ex, @compile code_String mode=mode_String args__)
        x = pynew()
        args = [isexpr(arg, :(=)) ? Expr(:kw, arg.args...) : arg for arg in args]
        filename = "$(st.src.file):$(st.src.line)"
        push!(st.inits, :($pycopy!($x, $pybuiltins.compile($code, filename=$filename, mode=$mode, $(args...)))))
        py_macro_assign(body, ans, x)
        return false

    # @exec code ...
    elseif @capture(ex, @exec code_String args__)
        args = [isexpr(arg, :(=)) ? Expr(:kw, arg.args...) : arg for arg in args]
        ex2 = Expr(:macrocall, Symbol("@compile"), st.src, code, :(mode = "exec"))
        ex2 = Expr(:call, :exec, ex2, args...)
        @gensym a
        ta = py_macro_lower(st, body, a, ex2)
        py_macro_del(body, a, ta)
        py_macro_assign(body, ans, nothing)
        return false

    # @eval code ...
    elseif @capture(ex, @eval code_String args__)
        args = [isexpr(arg, :(=)) ? Expr(:kw, arg.args...) : arg for arg in args]
        ex2 = Expr(:macrocall, Symbol("@compile"), st.src, code, :(mode = "eval"))
        ex2 = Expr(:call, :eval, ex2, args...)
        return py_macro_lower(st, body, ans, ex2)

    # begin; ...; end
    elseif isexpr(ex, :block)
        if isempty(ex.args)
            py_macro_assign(body, ans, nothing)
            return false
        else
            @gensym a
            for (i, aa) in enumerate(ex.args)
                if aa isa LineNumberNode
                    st.src = aa
                    push!(body, aa)
                    if i == length(ex.args)
                        py_macro_assign(body, ans, nothing)
                        return false
                    end
                elseif i == length(ex.args)
                    ta = py_macro_lower(st, body, ans, aa)
                    return ta
                else
                    ta = py_macro_lower(st, body, a, aa)
                    py_macro_del(body, a, ta)
                end
            end
            @assert false
        end

    # if x; ...; end
    elseif isexpr(ex, :if, :elseif)
        if length(ex.args) == 2
            return py_macro_lower(st, body, ans, Expr(ex.head, ex.args..., nothing))
        elseif length(ex.args) == 3
            ax, ay, az = ex.args
            @gensym x
            py_macro_lower_bool(st, body, x, ax)
            body2 = []
            ty = py_macro_lower(st, body2, ans, ay)
            body3 = []
            tz = py_macro_lower(st, body3, ans, az)
            t = ty || tz
            if t
                ty || push!(body2, :($ans = $Py($ans)))
                tz || push!(body3, :($ans = $Py($ans)))
            end
            push!(body, Expr(:if, x, Expr(:block, body2...), Expr(:block, body3...)))
            return t
        end

    # x && y
    elseif isexpr(ex, :&&)
        ax, ay = ex.args
        tx = py_macro_lower(st, body, ans, ax)
        body2 = []
        body3 = []
        tx && push!(body2, :($pydel!($ans)))
        ty = py_macro_lower(st, body2, ans, ay)
        t = tx || ty
        if t
            tx || push!(body3, :($ans = $Py($ans)))
            ty || push!(body2, :($ans = $Py($ans)))
        end
        push!(body, Expr(:if, :($pytruth($ans)), Expr(:block, body2...), Expr(:block, body3...)))
        return t

    # x || y
    elseif isexpr(ex, :||)
        ax, ay = ex.args
        tx = py_macro_lower(st, body, ans, ax)
        body2 = []
        body3 = []
        tx && push!(body3, :($pydel!($ans)))
        ty = py_macro_lower(st, body3, ans, ay)
        t = tx || ty
        if t
            tx || push!(body2, :($ans = $Py($ans)))
            ty || push!(body3, :($ans = $Py($ans)))
        end
        push!(body, Expr(:if, :($pytruth($ans)), Expr(:block, body2...), Expr(:block, body3...)))
        return t

    # while x; ...; end
    elseif isexpr(ex, :while)
        ax, ay = ex.args
        @gensym x y
        body2 = []
        py_macro_lower_bool(st, body2, x, ax)
        push!(body2, Expr(:if, :(!$x), :(break)))
        ty = py_macro_lower(st, body2, y, ay)
        py_macro_del(body2, y, ty)
        push!(body, Expr(:while, true, Expr(:block, body2...)))
        py_macro_assign(body, ans, nothing)
        return false

    # for x in y; ...; end
    elseif @capture(ex, for ax_ in ay_; az_; end)
        @gensym y i v z
        ty = py_macro_lower(st, body, y, ay)
        push!(body, :($i = $pyiter($y)))
        py_macro_del(body, y, ty)
        body2 = []
        push!(body2, :($v = $pynext($i)))
        push!(body2, Expr(:if, :($ispynull($v)), Expr(:block, :($pydel!($v)), :(break))))
        py_macro_lower_assign(st, body2, ax, v)
        py_macro_del(body2, v, true)
        tz = py_macro_lower(st, body2, z, az)
        py_macro_del(body2, z, tz)
        push!(body, Expr(:while, true, Expr(:block, body2...)))
        py_macro_del(body, i, true)
        py_macro_assign(body, ans, nothing)
        return false

    # import ...
    elseif isexpr(ex, :import)
        for aa in ex.args
            if isexpr(aa, :as)
                am, av = aa.args
                av isa Symbol || py_macro_err(st, av)
                py_macro_lower_import(st, body, av, am)
            elseif isexpr(aa, :.)
                @gensym m
                if length(aa.args) > 1
                    tm = py_macro_lower_import(st, body, m, aa)
                    py_macro_del(body, m, tm)
                end
                py_macro_lower_import(st, body, aa.args[1], Expr(:., aa.args[1]))
            elseif isexpr(aa, :(:))
                @gensym m
                py_macro_lower_import(st, body, m, aa.args[1])
                for ak in aa.args[2:end]
                    if isexpr(ak, :as)
                        ak, av = ak.args
                        av isa Symbol || py_macro_err(st, av)
                        py_macro_lower_nested_getattr(st, body, av, m, ak)
                    elseif isexpr(ak, :.)
                        py_macro_lower_nested_getattr(st, body, ak.args[end], m, ak)
                    else
                        py_macro_err(st, ak)
                    end
                end
                py_macro_del(body, m, true)
            else
                py_macro_err(st, aa)
            end
        end
        py_macro_assign(body, ans, nothing)
        return false

    # "...$foo..."
    elseif isexpr(ex, :string)
        args = [a isa String ? a : :(str($a)) for a in ex.args]
        return py_macro_lower(st, body, ans, :("".join(($(args...),))))

    end

    py_macro_err(st, ex)
end

function py_macro_lower_nested_getattr(st, body, ans, x, ks)
    @gensym m v
    isexpr(ks, :.) || py_macro_err(st, ks)
    length(ks.args) > 0 || py_macro_err(st, ks)
    py_macro_lower(st, body, ans, :($x.$(ks.args[1])))
    for k in ks.args[2:end]
        py_macro_lower(st, body, v, :($ans.$k))
        py_macro_del(body, ans, true)
        py_macro_assign(body, ans, v)
    end
end

function py_macro_lower_import(st, body, ans, ex)
    @gensym m
    isexpr(ex, :.) || py_macro_err(st, ex)
    tm = py_macro_lower(st, body, m, join(ex.args, "."))
    py_macro_assign(body, ans, :($pyimport($m)))
    py_macro_del(body, m, tm)
end

function py_macro_lower_bool(st, body, ans, ex)
    # TODO: special cases to avoid round-tripping bools from Julia to Python to Julia
    @gensym x
    tx = py_macro_lower(st, body, x, ex)
    py_macro_assign(body, ans, :($pytruth($x)))
    py_macro_del(body, x, tx)
    return
end

function py_macro_lower_assign(st, body, lhs, rhs::Symbol)
    if lhs isa Symbol
        push!(body, :($lhs = $Py($rhs)))
    elseif @capture(lhs, ax_[ak__])
        @gensym x k
        tx = py_macro_lower(st, body, x, ax)
        if length(ak) == 1
            tk = py_macro_lower(st, body, k, ak[1], flavour=:index)
        else
            tk = py_macro_lower(st, body, k, Expr(:tuple, ak...), flavour=:index)
        end
        push!(body, :($pysetitem($x, $k, $rhs)))
        py_macro_del(body, x, tx)
        py_macro_del(body, k, tk)
    elseif isexpr(lhs, :.)
        ax, ak = lhs.args
        @gensym x k
        tx = py_macro_lower(st, body, x, ax)
        if ak isa QuoteNode && ak.value isa Symbol
            tk = py_macro_lower(st, body, k, string(ak.value))
        else
            tk = py_macro_lower(st, body, k, ak)
        end
        push!(body, :($pysetattr($x, $k, $rhs)))
        py_macro_del(body, x, tx)
        py_macro_del(body, k, tk)
    else
        py_macro_err(st, lhs, "assignment")
    end
    return
end

function py_macro_lower_jl(st, ex)
    # TODO: expand out macros
    if @capture(ex, @py ax_)
        @gensym x
        body = []
        py_macro_lower(st, body, x, ax)
        return Expr(:block, body..., x)
    elseif ex isa Expr
        return Expr(ex.head, [py_macro_lower_jl(st, aa) for aa in ex.args]...)
    else
        return ex
    end
end

function py_macro(ex, mod, src)
    body = []
    @gensym ans
    st = PyMacroState(mod=mod, src=src)
    py_macro_lower(st, body, ans, ex)
    if !isempty(st.consts) || !isempty(st.inits)
        doinit = Ref(true)
        inits = copy(st.inits)
        for (v, x) in st.consts
            if v isa String
                push!(inits, :($pycopy!($x, $pystr_intern!($pystr($v)))))
            else
                push!(inits, :($pycopy!($x, $Py($v))))
            end
        end
        push!(inits, :($doinit[] = false))
        pushfirst!(body, Expr(:if, :($doinit[]), Expr(:block, inits...)))
    end
    pushfirst!(body, src)
    Expr(:block, body..., ans)
end

"""
    @py expr

Evaluate the given expression using Pythonic semantics.

For example:
- `f(x, y)` is translated to `pycall(f, x, y)`
- `x + y` is translated to `pyadd(x, y)`
- `x === y` is translated to `pyis(x, y)`
- `x.foo` is translated to `pygetattr(x, "foo")`

Compound statements such as `begin`, `if`, `while` and `for` are supported.

See the online documentation for more details.
"""
macro py(ex)
    esc(py_macro(ex, __module__, __source__))
end
export @py


text:
@testset "Test Saving and Loading PLS2 Models" begin

	Xtr        = [1 2;2 4;3 6;6 12;7 14.0]
	Ytr        = [2 2;4 4;6 6;12 12;14 14.0]
	Xt         = [4 8;5 10.0]
	model1    = PLSRegressor.fit(Xtr,Ytr,nfactors=2)
	pred1     = PLSRegressor.predict(model1,Xt)

	PLSRegressor.save(model1)
	model2    = PLSRegressor.load()

	pred2     = PLSRegressor.predict(model2,Xt)
    rm(PLSRegressor.MODEL_FILENAME)
	@test all(pred1 .== pred2)

end

@testset "PLS2 Prediction Tests (in sample)" begin

	@testset "Single Column Prediction Test" begin

		X        = [1; 2; 3.0]
		Y        = [1 1; 2 2; 3 3.0]
		model    = PLSRegressor.fit(X,Y,nfactors=1)
		pred     = PLSRegressor.predict(model,X)

		@test isequal(round.(pred),[1 1; 2 2; 3 3.0])

	end

	@testset "Constant Values Prediction Tests (Ax + b) | A=0, b=1 " begin

		X        = [1 3;2 1;3 2.0]
		Y        = [1 1; 1 1; 1 1.0]
		try
			PLSRegressor.fit(X,Y,nfactors=2)
		catch
			@test true
		end

	end

	@testset "Linear Prediction Tests " begin

		X        = [1 2; 2 4; 4 6.0]
		Y        = [4 2;6 4;8 6.0]
		model    = PLSRegressor.fit(X,Y,nfactors=2)
		pred     = PLSRegressor.predict(model,X)
		@test isequal(round.(pred),[4 2;6 4;8 6.0])

		X           = [1 -2; 2 -4; 4 -6.0]
		Y           = [-4 -2;-6 -4;-8 -6.0]
		model       = PLSRegressor.fit(X,Y,nfactors=2)
		pred        = PLSRegressor.predict(model,X)
		@test isequal(round.(pred),[-4 -2;-6 -4;-8 -6.0])

	end

end

@testset "PLS2 Pediction Tests (out of sample)" begin

	@testset "Linear Prediction Tests (Ax + b) | A>0" begin

		Xtr        = [1 2;2 4;3 6;6 12;7 14.0]
		Ytr        = [2 2;4 4;6 6;12 12;14 14.0]
		Xt         = [4 8;5 10.0]
		model    = PLSRegressor.fit(Xtr,Ytr,nfactors=2)
		pred     = PLSRegressor.predict(model,Xt)
		@test isequal(round.(pred),[8 8;10 10.0])

		Xtr        = [1 2;2 4;3 6;6 12;7 14.0]
		Ytr        = [2 4;4 6;6 8;12 14;14 16.0]
		Xt         = [4 8;5 10.0]
		model    = PLSRegressor.fit(Xtr,Ytr,nfactors=2)
		pred     = PLSRegressor.predict(model,Xt)
		@test isequal(round.(pred),[8 10;10 12.0])

	end

	@testset "Linear Prediction Tests (Ax + b) | A<0" begin

		Xtr        = [1 -2;2 -4;3 -6;6 -12;7 -14.0]
		Ytr        = [2 -2;4 -4;6 -6;12 -12;14 -14.0]
		Xt         = [4 -8;5 -10.0]
		model    = PLSRegressor.fit(Xtr,Ytr,nfactors=2)
		pred     = PLSRegressor.predict(model,Xt)
		@test isequal(round.(pred),[8 -8;10 -10.0])

		Xtr        = [1 -2;2 -4;3 -6;6 -12;7 -14.0]
		Ytr        = [2 -4;4 -6;6 -8;12 -14;14 -16.0]
		Xt         = [4 -8;5 -10.0]
		model    = PLSRegressor.fit(Xtr,Ytr,nfactors=2)
		pred     = PLSRegressor.predict(model,Xt)
		@test isequal(round.(pred),[8 -10;10 -12.0])

	end

end;

ext: jl


text:
# Computes classical Stein discrepancy based on the Langevin operator
# for univariate target distributions.
# Args:
#  points - n x d array of sample points
#  weights - n x 1 array of real-valued weights associated with sample points
#    (default: equal weights)
#  gradlogdensity - the gradlogdensity of the target distribution
#  solver - optimization program solver supported by JuMP that can solve
#   a quadratically constrained quadratic program (QCQP)
function langevin_classical_discrepancy(points=[],
                                        weights=fill(1/size(points,1), size(points,1)),
                                        gradlogdensity=nothing,
                                        solver=nothing,
                                        kwargs...)
    # we require that the solver is Gurobi or else the quadratic expressions
    # don't work
    try
        isa(solver, Main.Gurobi.GurobiSolver) ||
            error("The solver used must be Gurobi for classical discrepancies.")
    catch
        error("Gurobi must be used as the solver to solve the classical discrepancy.")
    end
    # get primary objects
    sample = SteinDiscrete(points, weights)
    points = sample.support
    weights = vec(sample.weights)
    n = length(weights)
    d = size(points,2)
    # check if d is 1
    (d == 1) || error("The classical discrepancy only works for univariate distributions.")
    # Objective coefficients for each g(x_i)
    gradlogdensities = zeros(n, d)
    for i in 1:n
        gradlogdensities[i,:] = gradlogdensity(points[i,:])
    end
    gobjcoefficients = vec(weights .* gradlogdensities)
    gradgobjcoefficients = weights
    # solve the optimization problem
    result = affine_classical_discrepancy(sample,
                                          gobjcoefficients,
                                          gradgobjcoefficients,
                                          solver,
                                          "langevin";
                                          kwargs...)
    result
end


text:
module Solaris

using CSV
using CUDA: CuArray
using DataFrames
using DiffEqFlux
using DiffEqFlux.GalacticOptim
using Flux
using JLD2
using Optim
using Plots
using PyCall
import Tracker

export device
export AbstractLayer,
       AbstractChain,
       Shortcut,
       Convex,
       ICNN,
       FastConvex,
       FastICNN
export sci_train,
       sci_train!,
       vis_train

include("widget.jl")
include("io.jl")
include("Layer/layer.jl")
include("train.jl")

const tf = PyNULL()

end


text:
@testset "Simulator: Flamingo" begin
    # Reference trajectory
	s = get_simulation("flamingo", "flat_2D_lc", "flat")
	model = s.model
	env = s.env

	model.μ_world = 0.1
	ref_traj = get_trajectory(s.model, s.env,
		joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait1.jld2"),
		load_type = :split_traj_alt)
	ContactImplicitMPC.update_friction_coefficient!(ref_traj, model, env)

	for t = 1:ref_traj.H
		r = ContactImplicitMPC.residual(model, env, ref_traj.z[t], ref_traj.θ[t], 0.0)
		@test norm(r) < 1.0e-4
	end

	# initial conditions
	q0 = SVector{model.dim.q}(ref_traj.q[1])
	q1 = SVector{model.dim.q}(ref_traj.q[2])

	h = ref_traj.h
	H = ref_traj.H
	model.μ_world = 0.9

	"""
	    PD tracking policy
	"""
	mutable struct PDPolicy <: ContactImplicitMPC.Policy
		model::ContactModel
	    traj::ContactTraj
		q0::AbstractVector
	    idx::Int
	    cnt::Int
	    N_sample::Int
	end

	function pd_policy(model, traj; N_sample = 1)
	    PDPolicy(model, traj, copy(traj.q[1]), 0, N_sample, N_sample)
	end

	function ContactImplicitMPC.policy(p::PDPolicy, x, traj, t)
	    # reset
	    if t == 1
	        p.idx = 0
	        p.cnt = p.N_sample
	    end

	    if p.cnt == p.N_sample
	        p.idx += 1
	        p.cnt = 0
			p.q0 .= copy(x)
	    end

	    p.cnt += 1

		u = p.traj.u[p.idx]
		# PD
		kp = 1.0 * ones(p.model.dim.u)
		kp[1] *= 100.0
		# kd = 10.0

		# u = Diagonal(kp) * B_func(p.model, x) * (x - p.traj.q[p.idx + 1])
		# @show u
		# u -= kd * B_func(p.model, x) * (x - traj.q[t]) / traj.h

	    return u ./ p.N_sample
	end

	p = pd_policy(model, ref_traj)

	# simulator
	sim = ContactImplicitMPC.simulator(s, q0, q1, h, 24,
	    p = p,
	    ip_opts = ContactImplicitMPC.InteriorPointOptions(
			r_tol = 1.0e-8,
			κ_tol = 1.0e-8),
	    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

	# simulate
	status = ContactImplicitMPC.simulate!(sim, verbose = false)
	@test status
end

# include(joinpath(module_dir(), "src", "dynamics", "flamingo", "visuals.jl"))
# vis = Visualizer()
# open(vis)
# anim = visualize_robot!(vis, model, sim.traj, name = :sim)
# anim = visualize_robot!(vis, model, ref_traj, anim = anim, name = :ref)


text:
module Run
using DirichletProcess
using Distributions
using Plots

function dp_samples(alpha, n)
    dp = DP(randn, alpha)
    samples = [dp_sample!(dp) for _ in 1:n]
    return samples
end

function hdp_samples(alpha1, alpha2, n)
    hdp = HDP(randn, alpha1, alpha2)
    samples = [hdp_sample!(hdp) for _ in 1:n]
    return samples
end

function dp()
    n = 10000
    fig = plot(; layout=(2, 2), legend=false)
    alphas = [1, 10, 100, 1000]
    for i in 1:4
        alpha = alphas[i]
        histogram!(dp_samples(alpha, n); subplot=i, normalize=true)
    end
    display(fig)

    return nothing
end

function hdp()
    n = 10000
    fig = plot(; layout=(2, 2), legend=false)
    # alphas = [1, 10, 100, 1000]
    alphas = [10, 10, 10, 10]
    for i in 1:4
        alpha = alphas[i]
        histogram!(hdp_samples(100, alpha, n); subplot=i, normalize=true)
    end
    display(fig)

    return nothing
end

end


text:
import Base: chop

# BLAS/linear algebra overrides

@inline dot(x...) = LinearAlgebra.dot(x...)
@inline dot(M::Int,a::Ptr{T},incx::Int,b::Ptr{T},incy::Int) where {T<:Union{Float64,Float32}} =
    BLAS.dot(M,a,incx,b,incy)
@inline dot(M::Int,a::Ptr{T},incx::Int,b::Ptr{T},incy::Int) where {T<:Union{ComplexF64,ComplexF32}} =
    BLAS.dotc(M,a,incx,b,incy)

dotu(f::StridedVector{T},g::StridedVector{T}) where {T<:Union{ComplexF32,ComplexF64}} =
    BLAS.dotu(f,g)
dotu(f::AbstractVector{Complex{Float64}},g::AbstractVector{N}) where {N<:Real} = dot(conj(f),g)
dotu(f::AbstractVector{N},g::AbstractVector{T}) where {N<:Real,T<:Number} = dot(f,g)

normalize!(w::AbstractVector) = rmul!(w,inv(norm(w)))
normalize!(w::Vector{T}) where {T<:BlasFloat} = normalize!(length(w),w)
normalize!(n,w::Union{Vector{T},Ptr{T}}) where {T<:Union{Float64,Float32}} =
    BLAS.scal!(n,inv(BLAS.nrm2(n,w,1)),w,1)
normalize!(n,w::Union{Vector{T},Ptr{T}}) where {T<:Union{ComplexF64,ComplexF32}} =
    BLAS.scal!(n,T(inv(BLAS.nrm2(n,w,1))),w,1)

flipsign(x,y) = Base.flipsign(x,y)
flipsign(x,y::Complex) = y==0 ? x : x*sign(y)

# Used for spaces not defined yet
struct UnsetNumber <: Number  end
promote_rule(::Type{UnsetNumber},::Type{N}) where {N<:Number} = N
promote_rule(::Type{Bool}, ::Type{UnsetNumber}) = Bool

# Test the number of arguments a function takes
hasnumargs(f,k) = applicable(f,zeros(k)...)

# fast implementation of isapprox with atol a non-keyword argument in most cases
isapprox_atol(a,b,atol;kwds...) = isapprox(a,b;atol=atol,kwds...)
isapprox_atol(a::Vec,b::Vec,atol::Real=0;kwds...) = isapprox_atol(collect(a),collect(b),atol;kwds...)
function isapprox_atol(x::Number, y::Number, atol::Real=0; rtol::Real=Base.rtoldefault(x,y))
    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y)))
end
function isapprox_atol(x::AbstractArray{T}, y::AbstractArray{S},atol::Real=0; rtol::Real=Base.rtoldefault(T,S), norm::Function=vecnorm) where {T<:Number,S<:Number}
    d = norm(x - y)
    if isfinite(d)
        return d <= atol + rtol*max(norm(x), norm(y))
    else
        # Fall back to a component-wise approximate comparison
        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol), zip(x, y))
    end
end

# The second case handles zero
isapproxinteger(x) = isapprox(x,round(Int,x))  || isapprox(x+1,round(Int,x+1))

# This creates ApproxFunBase.real, ApproxFunBase.eps and ApproxFunBase.dou
# which we override for default julia types
real(x...) = Base.real(x...)
real(::Type{UnsetNumber}) = UnsetNumber
real(::Type{Array{T,n}}) where {T<:Real,n} = Array{T,n}
real(::Type{Array{T,n}}) where {T<:Complex,n} = Array{real(T),n}
real(::Type{Vec{N,T}}) where {N,T<:Real} = Vec{N,T}
real(::Type{Vec{N,T}}) where {N,T<:Complex} = Vec{N,real(T)}

float(x) = Base.float(x)
Base.float(::UnsetNumber) = UnsetNumber()
Base.float(::Type{UnsetNumber}) = UnsetNumber
float(::Type{Array{T,N}}) where {T,N} = Array{float(T),N}
float(::Type{SVector{N,T}}) where {T,N} = SVector{N,float(T)}

eps(x...) = Base.eps(x...)
eps(x) = Base.eps(x)

eps(::Type{T}) where T<:Integer = zero(T)
eps(::Type{T}) where T<:Rational = zero(T)
eps(::T) where T<:Integer = eps(T)

eps(::Type{Complex{T}}) where {T<:Real} = eps(real(T))
eps(z::Complex{T}) where {T<:Real} = eps(abs(z))
eps(::Type{Dual{Complex{T}}}) where {T<:Real} = eps(real(T))
eps(z::Dual{Complex{T}}) where {T<:Real} = eps(abs(z))

eps(::Type{Vector{T}}) where {T<:Number} = eps(T)
eps(::Type{Vec{k,T}}) where {k,T<:Number} = eps(T)

isnan(x) = Base.isnan(x)
isnan(x::Vec) = map(isnan,x)

# BLAS

# implement muladd default
muladd(a,b,c) = a*b+c
muladd(a::Number,b::Number,c::Number) = Base.muladd(a,b,c)

for TYP in (:Float64,:Float32,:ComplexF64,:ComplexF32)
    @eval scal!(n::Integer,cst::$TYP,ret::DenseArray{T},k::Integer) where {T<:$TYP} =
            BLAS.scal!(n,cst,ret,k)
end

scal!(n::Integer,cst::BlasFloat,ret::DenseArray{T},k::Integer) where {T<:BlasFloat} =
    BLAS.scal!(n,convert(T,cst),ret,k)

function scal!(n::Integer,cst::Number,ret::AbstractArray,k::Integer)
    @assert k*n ≤ length(ret)
    @simd for j=1:k:k*(n-1)+1
        @inbounds ret[j] *= cst
    end
    ret
end

scal!(cst::Number,v::AbstractArray) = scal!(length(v),cst,v,1)

# Helper routines

function reverseeven!(x::AbstractVector)
    n = length(x)
    if iseven(n)
        @inbounds @simd for k=2:2:n÷2
            x[k],x[n+2-k] = x[n+2-k],x[k]
        end
    else
        @inbounds @simd for k=2:2:n÷2
            x[k],x[n+1-k] = x[n+1-k],x[k]
        end
    end
    x
end

function negateeven!(x::AbstractVector)
    @inbounds @simd for k = 2:2:length(x)
        x[k] *= -1
    end
    x
end

#checkerboard, same as applying negativeeven! to all rows then all columns
function negateeven!(X::AbstractMatrix)
    for j = 1:2:size(X,2)
        @inbounds @simd for k = 2:2:size(X,1)
            X[k,j] *= -1
        end
    end
    for j = 2:2:size(X,2)
        @inbounds @simd for k = 1:2:size(X,1)
            X[k,j] *= -1
        end
    end
    X
end

const alternatesign! = negateeven!

alternatesign(v::AbstractVector) = alternatesign!(copy(v))

alternatingvector(n::Integer) = 2*mod([1:n],2) .- 1

function alternatingsum(v::AbstractVector)
    ret = zero(eltype(v))
    s = 1
    @inbounds for k=1:length(v)
        ret+=s*v[k]
        s*=-1
    end

    ret
end

# Sum Hadamard product of vectors up to minimum over lengths
function mindotu(a::AbstractVector,b::AbstractVector)
    ret,m = zero(promote_type(eltype(a),eltype(b))),min(length(a),length(b))
    @inbounds @simd for i=m:-1:1 ret += a[i]*b[i] end
    ret
end

# efficiently resize a Matrix.  Note it doesn't change the input ptr
function unsafe_resize!(W::AbstractMatrix,::Colon,m::Integer)
    if m == size(W,2)
        W
    else
        n=size(W,1)
        reshape(resize!(vec(W),n*m),n,m)
    end
end

function unsafe_resize!(W::AbstractMatrix,n::Integer,::Colon)
    N=size(W,1)
    if n == N
        W
    elseif n < N
        W[1:n,:]
    else
        m=size(W,2)
        ret=Matrix{eltype(W)}(n,m)
        ret[1:N,:] = W
        ret
    end
end

function unsafe_resize!(W::AbstractMatrix,n::Integer,m::Integer)
    N=size(W,1)
    if n == N
        unsafe_resize!(W,:,m)
    else
        unsafe_resize!(unsafe_resize!(W,n,:),:,m)
    end
end

function pad!(f::AbstractVector{T},n::Integer) where T
	if n > length(f)
		append!(f,zeros(T,n - length(f)))
	else
		resize!(f,n)
	end
end

function pad(f::AbstractVector{T},n::Integer) where T
	if n > length(f)
	   ret=Vector{T}(undef, n)
	   ret[1:length(f)]=f
	   for j=length(f)+1:n
	       ret[j]=zero(T)
	   end
       ret
	else
        f[1:n]
	end
end

function pad(f::AbstractVector{Any},n::Integer)
	if n > length(f)
        Any[f...,zeros(n - length(f))...]
	else
        f[1:n]
	end
end

pad(x::Number, n::Int) = [x; zeros(typeof(x), n-1)]

function pad(v::AbstractVector,n::Integer,m::Integer)
    @assert m==1
    pad(v,n)
end

function pad(A::AbstractMatrix,n::Integer,m::Integer)
    T=eltype(A)
	if n <= size(A,1) && m <= size(A,2)
        A[1:n,1:m]
	elseif n==0 || m==0
	   Matrix{T}(undef,n,m)  #fixes weird julia bug when T==None
    else
        ret = Matrix{T}(undef,n,m)
        minn=min(n,size(A,1))
        minm=min(m,size(A,2))
        for k=1:minn,j=1:minm
            @inbounds ret[k,j]=A[k,j]
        end
        for k=minn+1:n,j=1:minm
            @inbounds ret[k,j]=zero(T)
        end
        for k=1:n,j=minm+1:m
            @inbounds ret[k,j]=zero(T)
        end
        for k=minn+1:n,j=minm+1:m
            @inbounds ret[k,j]=zero(T)
        end

        ret
	end
end

pad(A::AbstractMatrix,::Colon,m::Integer) = pad(A,size(A,1),m)
pad(A::AbstractMatrix,n::Integer,::Colon) = pad(A,n,size(A,2))

function pad(v, ::PosInfinity)
    if isinf(length(v))
        v
    else
        Vcat(v, Zeros{Int}(∞))
    end
end

function pad(v::AbstractVector{T}, ::PosInfinity) where T
    if isinf(length(v))
        v
    else
        Vcat(v, Zeros{T}(∞))
    end
end

#TODO:padleft!

function padleft(f::AbstractVector,n::Integer)
	if (n > length(f))
        [zeros(n - length(f)); f]
	else
        f[end-n+1:end]
	end
end

##chop!
function chop!(c::AbstractVector,tol::Real)
    @assert tol >= 0

    for k=length(c):-1:1
        if abs(c[k]) > tol
            resize!(c,k)
            return c
        end
    end

    resize!(c,0)
    c
end

chop(f::AbstractVector,tol) = chop!(copy(f),tol)
chop!(f::AbstractVector) = chop!(f,eps())

function chop!(A::AbstractArray,tol)
    for k=size(A,1):-1:1
        if norm(A[k,:])>tol
            A=A[1:k,:]
            break
        end
    end
    for k=size(A,2):-1:1
        if norm(A[:,k])>tol
            A=A[:,1:k]
            break
        end
    end
    return A
end
chop(A::AbstractArray,tol)=chop!(A,tol)#replace by chop!(copy(A),tol) when chop! is actually in-place.

## interlace

function interlace(v::Union{Vector{Any},Tuple})
    #determine type
    T=Float64
    for vk in v
        if isa(vk,Vector{Complex{Float64}})
            T=Complex{Float64}
        end
    end
    b=Vector{Vector{T}}(undef, length(v))
    for k=1:length(v)
        b[k]=v[k]
    end
    interlace(b)
end

function interlace(a::AbstractVector{S},b::AbstractVector{V}) where {S<:Number,V<:Number}
    na=length(a);nb=length(b)
    T=promote_type(S,V)
    if nb≥na
        ret=zeros(T,2nb)
        ret[1:2:1+2*(na-1)]=a
        ret[2:2:end]=b
        ret
    else
        ret=zeros(T,2na-1)
        ret[1:2:end]=a
        if !isempty(b)
            ret[2:2:2+2*(nb-1)]=b
        end
        ret
    end
end

function interlace(a::AbstractVector,b::AbstractVector)
    na=length(a);nb=length(b)
    T=promote_type(eltype(a),eltype(b))
    if nb≥na
        ret=Vector{T}(undef, 2nb)
        ret[1:2:1+2*(na-1)]=a
        ret[2:2:end]=b
        ret
    else
        ret=Vector{T}(undef, 2na-1)
        ret[1:2:end]=a
        if !isempty(b)
            ret[2:2:2+2*(nb-1)]=b
        end
        ret
    end
end

### In-place O(n) interlacing

function highestleader(n::Int)
    i = 1
    while 3i < n i *= 3 end
    i
end

function nextindex(i::Int,n::Int)
    i <<= 1
    while i > n
        i -= n + 1
    end
    i
end

function cycle_rotate!(v::AbstractVector, leader::Int, it::Int, twom::Int)
    i = nextindex(leader, twom)
    while i != leader
        idx1, idx2 = it + i - 1, it + leader - 1
        @inbounds v[idx1], v[idx2] = v[idx2], v[idx1]
        i = nextindex(i, twom)
    end
    v
end

function right_cyclic_shift!(v::AbstractVector, it::Int, m::Int, n::Int)
    itpm = it + m
    itpmm1 = itpm - 1
    itpmpnm1 = itpmm1 + n
    reverse!(v, itpm, itpmpnm1)
    reverse!(v, itpm, itpmm1 + m)
    reverse!(v, itpm + m, itpmpnm1)
    v
end

"""
This function implements the algorithm described in:

    P. Jain, "A simple in-place algorithm for in-shuffle," arXiv:0805.1598, 2008.
"""
function interlace!(v::AbstractVector,offset::Int)
    N = length(v)
    if N < 2 + offset
        return v
    end

    it = 1 + offset
    m = 0
    n = 1

    while m < n
        twom = N + 1 - it
        h = highestleader(twom)
        m = h > 1 ? h÷2 : 1
        n = twom÷2

        right_cyclic_shift!(v,it,m,n)

        leader = 1
        while leader < 2m
            cycle_rotate!(v, leader, it, 2m)
            leader *= 3
        end

        it += 2m
    end
    v
end

## slnorm gives the norm of a slice of a matrix

function slnorm(u::AbstractMatrix,r::AbstractRange,::Colon)
    ret = zero(real(eltype(u)))
    for k=r
        @simd for j=1:size(u,2)
            #@inbounds
            ret=max(norm(u[k,j]),ret)
        end
    end
    ret
end

function slnorm(m::AbstractMatrix,kr::AbstractRange,jr::AbstractRange)
    ret=zero(real(eltype(m)))
    for j=jr
        nrm=zero(real(eltype(m)))
        for k=kr
            @inbounds nrm+=abs2(m[k,j])
        end
        ret=max(sqrt(nrm),ret)
    end
    ret
end

slnorm(m::AbstractMatrix,kr::AbstractRange,jr::Integer) = slnorm(m,kr,jr:jr)
slnorm(m::AbstractMatrix,kr::Integer,jr::AbstractRange) = slnorm(m,kr:kr,jr)

function slnorm(B::BandedMatrix{T},r::AbstractRange,::Colon) where T
    ret = zero(real(T))
    m=size(B,2)
    for k=r
        @simd for j=max(1,k-B.l):min(k+B.u,m)
            #@inbounds
            ret=max(norm(B[k,j]),ret)
        end
    end
    ret
end

slnorm(m::AbstractMatrix,k::Integer,::Colon) = slnorm(m,k,1:size(m,2))
slnorm(m::AbstractMatrix,::Colon,j::Integer) = slnorm(m,1:size(m,1),j)

## Infinity

Base.isless(x::Block{1}, y::PosInfinity) = isless(Int(x), y)
Base.isless(x::PosInfinity, y::Block{1}) = isless(x, Int(y))

## BandedMatrix

pad!(A::BandedMatrix,n,::Colon) = pad!(A,n,n+A.u)  # Default is to get all columns
columnrange(A,row::Integer) = max(1,row-bandwidth(A,1)):row+bandwidth(A,2)

## Store iterator
mutable struct CachedIterator{T,IT}
    iterator::IT
    storage::Vector{T}
    state
    length::Int
end

CachedIterator{T,IT}(it::IT, state) where {T,IT} = CachedIterator{T,IT}(it,T[],state,0)
CachedIterator(it::IT) where IT = CachedIterator{eltype(it),IT}(it, ())

function resize!(it::CachedIterator,n::Integer)
    m = it.length
    if n > m
        if n > length(it.storage)
            resize!(it.storage,2n)
        end

        @inbounds for k = m+1:n
            xst = iterate(it.iterator,it.state...)
            if xst == nothing
                it.length = k-1
                return it
            end
            it.storage[k] = xst[1]
            it.state = (xst[2],)
        end

        it.length = n
    end
    it
end

eltype(it::CachedIterator{T}) where {T} = T

iterate(it::CachedIterator) = iterate(it,1)
function iterate(it::CachedIterator,st::Int)
    if  st == it.length + 1 && iterate(it.iterator,it.state...) == nothing
        nothing
    else
        (it[st],st+1)
    end
end

function getindex(it::CachedIterator, k)
    mx = maximum(k)
    if mx > length(it) || mx < 1
        throw(BoundsError(it,k))
    end
    resize!(it,isempty(k) ? 0 : mx).storage[k]
end

function findfirst(f::Function,A::CachedIterator)
    k=1
    for c in A
        if f(c)
            return k
        end
        k+=1
    end
    return 0
end

function findfirst(A::CachedIterator,x)
    k=1
    for c in A
        if c == x
            return k
        end
        k+=1
    end
    return 0
end

length(A::CachedIterator) = length(A.iterator)

## nocat
vnocat(A...) = Base.vect(A...)
hnocat(A...) = Base.typed_hcat(mapreduce(typeof,promote_type,A),A...)
hvnocat(rows,A...) = Base.typed_hvcat(mapreduce(typeof,promote_type,A),rows,A...)
macro nocat(x)
    ex = expand(x)
    if ex.args[1] == :vcat
        ex.args[1] = :(ApproxFunBase.vnocat)
    elseif ex.args[1] == :hcat
        ex.args[1] = :(ApproxFunBase.hnocat)
    else
        @assert ex.args[1] == :hvcat
        ex.args[1] = :(ApproxFunBase.hvnocat)
    end
    esc(ex)
end

## Dynamic functions

struct DFunction <: Function
    f
end
(f::DFunction)(args...) = f.f(args...)

hasnumargs(f::DFunction, k) = hasnumargs(f.f, k)

dynamic(f) = f
dynamic(f::Function) = DFunction(f) # Assume f has to compile every time

# Matrix inputs

## conv

conv(x::AbstractVector, y::AbstractVector) = DSP.conv(x, y)
@generated function conv(x::SVector{N}, y::SVector{M}) where {N,M}
    NM = N+M-1
    quote
        convert(SVector{$NM}, DSP.conv(Vector(x), Vector(y)))
    end
end

conv(x::SVector{1}, y::SVector{1}) = x.*y
conv(x::AbstractVector, y::SVector{1}) = x*y[1]
conv(y::SVector{1}, x::AbstractVector) = y[1]*x
conv(x::AbstractFill, y::SVector{1}) = x*y[1]
conv(y::SVector{1}, x::AbstractFill) = y[1]*x
conv(x::AbstractFill, y::AbstractFill) = DSP.conv(x, y)

## BlockInterlacer
# interlaces coefficients by blocks
# this has the property that all the coefficients of a block of a subspace
# are grouped together, starting with the first bloc
#
# TODO: cache sums

struct BlockInterlacer{DMS<:Tuple}
    blocks::DMS
end

const TrivialInterlacer{d} = BlockInterlacer{NTuple{d,<:Ones}}

BlockInterlacer(v::AbstractVector) = BlockInterlacer(tuple(v...))

Base.eltype(it::BlockInterlacer) = Tuple{Int,Int}

dimensions(b::BlockInterlacer) = map(sum,b.blocks)
dimension(b::BlockInterlacer,k) = sum(b.blocks[k])
Base.length(b::BlockInterlacer) = mapreduce(sum,+,b.blocks)

# the state is always (whichblock,curblock,cursubblock,curcoefficients)
# start(it::BlockInterlacer) = (1,1,map(start,it.blocks),ntuple(zero,length(it.blocks)))

# are all Ints, so finite dimensional
function done(it::BlockInterlacer,st)
    for k=1:length(it.blocks)
        if st[end][k] < sum(it.blocks[k])
            return false
        end
    end
    return true
end

iterate(it::BlockInterlacer) =
    iterate(it, (1,1,ntuple(_ -> tuple(), length(it.blocks)),
            ntuple(zero,length(it.blocks))))

function iterate(it::BlockInterlacer, (N,k,blkst,lngs))
    done(it, (N,k,blkst,lngs)) && return nothing

    if N > length(it.blocks)
        # increment to next block
        blkst = map(function(blit,blst)
                xblst = iterate(blit, blst...)
                xblst == nothing ? blst : (xblst[2],)
            end,it.blocks,blkst)
        return iterate(it,(1,1,blkst,lngs))
    end

    Bnxtb = iterate(it.blocks[N],blkst[N]...)  # B is block size

    if Bnxtb == nothing
        # increment to next N
        return iterate(it,(N+1,1,blkst,lngs))
    end

    B,nxtb = Bnxtb

    if k > B
        #increment to next N
        return iterate(it,(N+1,1,blkst,lngs))
    end

    lngs = tuple(lngs[1:N-1]...,lngs[N]+1,lngs[N+1:end]...)
    return (N,lngs[N]),(N,k+1,blkst,lngs)
end

cache(Q::BlockInterlacer) = CachedIterator(Q)


text:
# Gvec is of shape (3,Ng)
# the argument lmax is different from the corresponding
# Fortran subroutine
function ylmr2( lmax::Int, Gvec::Array{Float64,2} )

  Ng = size(Gvec)[2]
  lmax2 = (lmax+1)^2

  ylm = zeros( Float64, Ng, lmax2 )

  gg = zeros(Float64,Ng)

  for ig = 1:Ng
    gg[ig] = Gvec[1,ig]^2 + Gvec[2,ig]^2 + Gvec[3,ig]^2
  end

  const SMALL = 1.e-9
  const FPI = 4.0*pi

  if (Ng < 1) || (lmax2 < 1)
    return ylm
  end

  if lmax == 0
    ylm[:,1] = sqrt(1.0/FPI)
    return ylm
  end

  cost = zeros(Float64,Ng)
  sent = zeros(Float64,Ng)
  phi  = zeros(Float64,Ng)
  Q    = zeros(Float64,Ng,lmax+1,lmax+1)  # originally 0:lmax

  for ig = 1:Ng

    gmod = sqrt( gg[ig] )

    if gmod < SMALL
      cost[ig] = 0.d0
    else
      cost[ig] = Gvec[3,ig]/gmod
    end

    # beware the arc tan, it is defined modulo pi
    # NOTE: verify this in Julia
    if Gvec[1,ig] > SMALL
      phi[ig] = atan( Gvec[2,ig]/Gvec[1,ig] )
      #
    elseif Gvec[1,ig] < -SMALL
      phi[ig] = atan( Gvec[2,ig]/Gvec[1,ig] ) + pi
      #
    else
      phi(ig) = sign( pi/2., Gvec[2,ig] )  # XXX CHECK THIS
    end

    sent[ig] = sqrt( max( 0.0, 1.0 - cost[ig]^2) )

  end # ig

  #
  # Q(:,l,m) are defined as sqrt ((l-m)!/(l+m)!) * P(:,l,m) where
  # P(:,l,m) are the Legendre Polynomials (0 <= m <= l)
  #
  lm = 0
  for l = 0:lmax
    @printf("l = %d\n", l)
    c = sqrt( (2*l+1) / FPI )

    if l == 0
      for ig = 1:Ng
        Q[ig,1,1] = 1.0  # l=0 -> 1
      end

    elseif l == 1
      for ig = 1:Ng
        Q[ig,2,1] =  cost[ig]
        Q[ig,2,2] = -sent[ig]/sqrt(2.0)
      end

    else
      # recursion on l for Q(:,l,m)
      for m = 0:(l-2)
        println("Enter this ...")
        for ig = 1:Ng
          Q[ig,l+1,m+1] = cost[ig]*(2*l-1)/sqrt(l*l-m*m) * Q[ig,l,m+1]
                         - sqrt((l-1)*(l-1)-m*m)/sqrt(l*l-m*m) * Q[ig,l-1,m+1]
        end
        println("sum(Q) = ", sum(Q[:,l+1,m+1]))
      end
      for ig = 1:Ng
        Q[ig,l+1,l] = cost[ig] * sqrt(2*l-1) * Q[ig,l,l]
      end
      for ig = 1:Ng
        Q[ig,l+1,l+1] = - sqrt(2*l-1)/sqrt(2*l)*sent[ig]*Q[ig,l,l]
      end

    end

    #
    # Y_lm, m = 0
    #
    lm = lm + 1
    @printf("m = 0: lm = %d\n", lm)
    for ig = 1:Ng
      ylm[ig, lm] = c * Q[ig,l+1,1]
    end

    for m = 1:l
      #
      # Y_lm, m > 0
      #
      lm = lm + 1
      @printf("m > 0: lm = %d\n", lm)
      for ig = 1:Ng
        ylm[ig, lm] = c * sqrt(2.0) * Q[ig,l+1,m+1] * cos( m*phi[ig] )
      end
      #
      # Y_lm, m < 0
      #
      lm = lm + 1
      @printf("m < 0: lm = %d\n", lm)
      for ig = 1:Ng
        ylm[ig,lm] = c * sqrt(2.0) * Q[ig,l+1,m+1] * sin(m*phi[ig])
      end
    end

  end # lm

  return ylm

end


text:
##### version 2 -> include the leg first

type Classical2dFractalSimulator{T} <: Classical2dSimulator{T}
	dimM::Int
	wholeiteration::Int
	countiteration::Int
	normalizationFactor::Array{T,2}
	coefficients::Array{T,2}
	lattice::Classical2dFractalLattice{T}

	function Classical2dFractalSimulator{T}(lattice::Classical2dFractalLattice{T}, dimM::Int, wholeiteration::Int)
		this = new{T}()
		this.dimM = dimM
		this.wholeiteration = wholeiteration
		this.lattice = lattice
		this.normalizationFactor = Array{T,2}(wholeiteration+1,4)
		this.coefficients = Array{T,2}(wholeiteration+1,4)
		initializeCount!(this)
		return this
	end
end

Classical2dFractalSimulator{T}(lattice::Classical2dFractalLattice{T}, dimM::Int, wholeiteration::Int) =
Classical2dFractalSimulator{T}(lattice,dimM,wholeiteration)

function (simulator::Classical2dFractalSimulator){T}(temperature::T = getTemperature(simulator), externalfield::T = getExternalfield(simulator))
	setEnvParameters!(simulator, temperature, externalfield)
	initializeCount!(simulator)
	countUp!(simulator)
	normalizeTensor!(simulator)
	initializeCoefficients!(simulator)
	while true
		countUp!(simulator)
		renormalize!(simulator.lattice, getDimM(simulator))
		normalizeTensor!(simulator)
		updateCoefficients!(simulator)

		##### debug
		##### maybe the norm gets too small
		# println(maximum(abs(tenT)))
		# println(sum(isnan(tenT)))
		# norm = getNorm(simulator,3,getCount(simulator))
		# println(norm)
		##### debug/
		if getCount(simulator) > getWholeiteration(simulator)
			break
		end
	end
	freeenergy = getFreeEnergy(simulator)
	magnetization = getExpectationValue(simulator)
	if isPottsModel(simulator)
		NOS = convert(T,getStates(simulator))
		magnetization = (NOS * magnetization - one(T))/(NOS - one(T))
	end

	return freeenergy, magnetization
end

function getFreeEnergy{T}(simulator::Classical2dFractalSimulator{T})
	numberofsites = sum(getCoefficient(simulator, 0, 1))

	tensorT = getTensorT(simulator)[1]
	termcorrection = log(traceTensorTPeriodic(tensorT))

	freeenergy = - getFirstTerm(simulator) * (
											termFromNormFactor(simulator) +
										    termcorrection
										   ) / numberofsites
	return freeenergy
end

function termFromNormFactor{T}(simulator::Classical2dFractalSimulator{T})
	termNorm = zero(T)
	for i in 1:getCount(simulator)
		for j in 1:4
			termNorm += getCoefficient(simulator, j, i) * log(getNorm(simulator, j, i))
		end
	end
	return termNorm
end

function normalizeTensor!{T}(simulator::Classical2dFractalSimulator{T})
	tensorT, tensorTtilde = getTensorT(simulator)
	newTensorT, normT = normalizeTensor(tensorT)
	newTensorTtilde = tensorTtilde ./ normT
	setTensorT!(simulator, newTensorT)
	setTensorT!(simulator, newTensorTtilde; tilde = true)

	newTensorPx, normPx = normalizeTensor(getTensorP(simulator, 1))
	setTensorP!(simulator, 1, newTensorPx)
	newTensorPy, normPy = normalizeTensor(getTensorP(simulator, 2))
	setTensorP!(simulator, 2, newTensorPy)

	newTensorQ, normQ = normalizeTensor(getTensorQ(simulator))
	setTensorQ!(simulator, newTensorQ)

	setNorm!(simulator, normT, 1)
	setNorm!(simulator, normPx, 2)
	setNorm!(simulator, normPy, 3)
	setNorm!(simulator, normQ, 4)
end

function setNorm!{T}(simulator::Classical2dFractalSimulator{T}, norm::T, which::Int)
	iteration = getCount(simulator)
	setNorm!(simulator, norm, which, iteration)
end

function initializeCoefficients!{T}(simulator::Classical2dFractalSimulator{T})
	wholeiteration = getWholeiteration(simulator)
	setCoefficient!(simulator, 1.0, 1, wholeiteration + 1)
	setCoefficient!(simulator, 0.0, 2, wholeiteration + 1)
	setCoefficient!(simulator, 0.0, 3, wholeiteration + 1)
	setCoefficient!(simulator, 0.0, 4, wholeiteration + 1)
end

function updateCoefficients!{T}(simulator::Classical2dFractalSimulator{T})
	iteration = getWholeiteration(simulator) +2 -getCount(simulator)
	currentCoef = getCoefficient(simulator, 0, iteration+1)
	newCoefT = updateCoefT(currentCoef)
	newCoefPx = updateCoefPx(currentCoef, getLegextension(simulator))
	newCoefPy = updateCoefPy(currentCoef, getLegextension(simulator))
	newCoefQ = updateCoefQ(currentCoef)

	setCoefficient!(simulator, newCoefT, 1, iteration)
	setCoefficient!(simulator, newCoefPx, 2, iteration)
	setCoefficient!(simulator, newCoefPy, 3, iteration)
	setCoefficient!(simulator, newCoefQ, 4, iteration)
end

function updateCoefT{T}(previousCoef::Array{T,1})
	newCoefT = 4.0 * sum(previousCoef)
	return newCoefT
end

function updateCoefPx{T}(previousCoef::Array{T,1}, legExtension::Int)
	newCoefPx = 4.0 * legExtension * previousCoef[1] +
		(4.0 * legExtension - 1.0) * previousCoef[2] +
		(4.0 * legExtension - 1.0) * previousCoef[3] +
		(4.0 * legExtension -2.0) * previousCoef[4]
	return newCoefPx
end

function updateCoefPy{T}(previousCoef::Array{T,1}, legExtension::Int)
	return updateCoefPx(previousCoef, legExtension)
end

function updateCoefQ{T}(previousCoef::Array{T,1})
	newCoefQ = 2.0 * previousCoef[2] +
		2.0 * previousCoef[3] +
		4.0 * previousCoef[4]
	return newCoefQ
end

#---
# test Function
function testClassical2dFractalSimulator()
	lattice = testClassical2dFractalLattice()
	dimM = 5
	wholeiteration = 15
	simulator = Classical2dFractalSimulator(lattice, dimM, wholeiteration)
	return simulator
end


text:
##########
##########
# Create series of networks to run through simulations
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-05-29
##########
##########

# set up =======================================================================
using DrWatson
@quickactivate "trait-based-rewiring"

include(srcdir("00_species_list_creation.jl"))

using DelimitedFiles

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# create different combinations of networks ====================================

# the different values that have different options are:
#   richness shape
#   biomass shape

all_network_data = Array{Any}(undef, (450000, 6))
all_networks = Array{Int}(undef, (450000, 51))
all_species_lists = Array{Any}(undef, (450000, 11))

"""
    create_networks(j,
        network_id,
        all_network_data,
        all_networks,
        all_species_lists
    )

For simplicity, put the process of making all of the 9000 networks in a function
"""
function create_networks(
    reps::Int,
    j::Int,
    all_network_data::Array,
    all_networks::Array,
    all_species_lists::Array
    )

    # make 1000 networks for each of the 9 combinations (3x3 of richness and
    # biomass) to have multiple starting points
    network_id = 1
    for i in 1:reps
        for richness in ["top-heavy", "bottom-heavy", "uniform"]
            for biomass in ["pyramid", "uniform", "inverted"]
                temp = Array{Any}(undef, (50, 6))
                # get species list and network
                species_list, network = network_draw(50, richness)
                # save network
                all_networks[j:j+49, 1] = repeat([network_id], 50)
                all_networks[j:j+49, 2:51] = network
                # get biomass
                species_list = initialize_biomass(species_list, biomass)
                # append the network id to each species list item
                species_list = hcat(species_list, fill(network_id, 50))
                # put the species list in the all_species_lists
                all_species_lists[j:j+49, :] = species_list
                # put each species with the network data in all_network_data
                for k in 1:size(species_list,1)
                    temp[k, 1] = network_id
                    temp[k, 2] = richness
                    temp[k, 3] = biomass
                    temp[k, 4] = species_list[k, body_size_pos]
                    temp[k, 5] = species_list[k, biomass_pos]
                    temp[k, 6] = species_list[k, organism_type_pos]
                end
                all_network_data[(j:j+49), :] = temp
                # increment
                j = j + 50
                network_id = network_id + 1
            end
        end
        if (i % 5 == 0) print(i, " out of 1000 --") end
    end
    return all_network_data, all_networks, all_species_lists
end

# run function =================================================================
all_network_data_done, all_networks_done, all_species_lists_done = create_networks(
    1000, 1, all_network_data, all_networks, all_species_lists)
# add alive/dead
all_species_lists_done = hcat(all_species_lists_done, ones(size(all_species_lists_done, 1)))

# save files ===================================================================

writedlm(datadir("./intermediate-files/all_networks_for_simulations.csv"),
        all_networks_done, ",")
writedlm(datadir("./intermediate-files/species_lists_for_simulations.csv"),
        all_species_lists_done, ",")
writedlm(datadir("./intermediate-files/all_network_data_for_simulations.csv"),
        all_network_data_done, ",")

# test that files are equivalent when read back in
networks_test = readdlm(datadir("./intermediate-files/all_networks_for_simulations.csv"), ',', Any)
if networks_test == all_networks_done
    print("networks are equivalent")
else
    ArgumentError("networks are NOT equivalent")
end

network_data_test = readdlm(datadir("./intermediate-files/all_network_data_for_simulations.csv"), ',', Any)
if network_data_test == all_network_data_done
    print("network data are equivalent")
else
    ArgumentError("network data are NOT equivalent")
end

species_list_test = readdlm(datadir("./intermediate-files/species_lists_for_simulations.csv"), ',', Any)
i = 1
while i <= size(all_species_lists_done, 1)
    if species_list_test[i,:] === all_species_lists_done[i,:]
        error("species lists are NOT equivalent at ", i)
    end
    i = i + 1
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
@testset "Genetic Programming" begin
    Random.seed!(9874984737486);
    pop = 10
    terms = Terminal[:x, :y, rand]
    funcs = Function[+,-,*,/]

    t = TreeGP(pop, terms, funcs, maxdepth=2)
    @test Evolutionary.population_size(t) == pop
    @test sort(Evolutionary.terminals(t)) == [:x, :y]
    @testset for (term, dim) in t.terminals
        @test dim == 1
    end
    @testset for (func, arity) in t.functions
        @test arity == 2
    end
    @test summary(t) == "TreeGP[P=10,Parameter[x,y],Function[*, +, /, -]]"

    popexp = Evolutionary.initial_population(t);
    @test length(popexp) == pop

    Random.seed!(9874984737484)
    ft = rand(t, method=:full)
    @test Evolutionary.nodes(ft) == 7
    @test Evolutionary.height(ft) == 2
    gt = rand(t, method=:grow)
    @test Evolutionary.nodes(gt) == 5
    @test Evolutionary.height(gt) == 2
    @test gt[0] == gt
    @test gt[1] == :x
    @test gt[2] == :x
    @test gt[4] == Expr(:call, /, :x, :x)

    # simplification
    @test Expr(:call, -, :x, :x) |> Evolutionary.simplify! == 0
    @test Expr(:call, /, :x, :x) |> Evolutionary.simplify! == 1
    @test Expr(:call, *, 0, :x) |> Evolutionary.simplify! == 0
    @test Expr(:call, *, :x, 0) |> Evolutionary.simplify! == 0
    @test Expr(:call, /, 0, :x) |> Evolutionary.simplify! == 0
    @test Expr(:call, /, 1, 0) |> Evolutionary.simplify! == 1
    @test Expr(:call, +, 0, :x) |> Evolutionary.simplify! == :x
    @test Expr(:call, +, :x, 0) |> Evolutionary.simplify! == :x
    @test Expr(:call, -, :x, 0) |> Evolutionary.simplify! == :x
    @test Expr(:call, +, :x, :x) |> Evolutionary.simplify! == Expr(:call, *, 2, :x)
    @test Expr(:call, +, 1, Expr(:call, -, 1, :x) ) |> Evolutionary.simplify! == Expr(:call, -, 2, :x)
    @test Expr(:call, -, 1, Expr(:call, +, 1, :x) ) |> Evolutionary.simplify!  == Expr(:call, +, 0, :x)
    @test Expr(:call, *, Expr(:call, *, :x, :y), Expr(:call, /, :z, :x)) |> Evolutionary.simplify! == Expr(:call, *, :y, :z)
    @test Expr(:call, *, Expr(:call, *, :x, :y), Expr(:call, /, :z, :y)) |> Evolutionary.simplify! == Expr(:call, *, :x, :z)
    @test Expr(:call, *, Expr(:call, /, :z, :x), Expr(:call, *, :x, :y)) |> Evolutionary.simplify! == Expr(:call, *, :z, :y)
    @test Expr(:call, *, Expr(:call, /, :z, :y), Expr(:call, *, :x, :y)) |> Evolutionary.simplify! == Expr(:call, *, :z, :x)
    @test Expr(:call, *, Expr(:call, *, 2, :y), Expr(:call, /, :z, 2)) |> Evolutionary.simplify! == Expr(:call, *, :y, :z)
    @test Expr(:call, *, Expr(:call, *, :y, 2), Expr(:call, /, :z, 2)) |> Evolutionary.simplify! ==  Expr(:call, *, :y, :z)
    @test Expr(:call, *, Expr(:call, /, :z, 2), Expr(:call, *, 2, :y)) |> Evolutionary.simplify! == Expr(:call, *, :z, :y)
    @test Expr(:call, *, Expr(:call, /, :z, 2), Expr(:call, *, :y, 2)) |> Evolutionary.simplify! ==  Expr(:call, *, :z, :y)
    @test Expr(:call, +, Expr(:call, -, :x, 1), 1 ) |> Evolutionary.simplify! == Expr(:call, -, :x, 2)
    @test Expr(:call, -, Expr(:call, +, :x, 1), 1 ) |> Evolutionary.simplify! == Expr(:call, +, :x, 0)

    # evaluation
    ex = Expr(:call, +, 1, :x) |> Evolutionary.Expression
    xs = rand(10)
    @test ex(xs[1]) == xs[1]+1
    @test ex.(xs) == xs.+1

    depth = 5
    fitfun(x) = x*x + x + 1.0
    function fitobj(expr)
        rg = -5.0:0.5:5.0
        ex = Evolutionary.Expression(expr)
        sum(v->isnan(v) ? 1.0 : v, abs2.(fitfun.(rg) - ex.(rg)) )/length(rg) |> sqrt
    end

    Random.seed!(9874984737426);
    res = Evolutionary.optimize(fitobj,
        TreeGP(50, Terminal[:x, randn], Function[+,-,*,/],
            mindepth=1,
            maxdepth=depth,
            optimizer = GA(
                selection = uniformranking(5),
                ɛ = 0.1,
                mutationRate = 0.8,
                crossoverRate = 0.2,
            ),
        )
    )
    @test minimum(res) < 1

end


text:
import Base.Order: lt, Ordering
export LenLex, WreathOrder, RecursivePathOrder

"""
    WordOrdering <: Ordering
Abstract type representing word orderings.

The subtypes of `WordOrdering` should contain a field `A` storing the `Alphabet`
over which a particular order is defined. Moreover, an `Base.lt` method should be
defined to compare whether one word is less than the other (in the ordering
defined).
"""
abstract type WordOrdering <: Ordering end

alphabet(o::WordOrdering) = o.A
Base.:(==)(o1::T, o2::T) where {T<:WordOrdering} = alphabet(o1) == alphabet(o2)

"""
    struct LenLex{T} <: WordOrdering

Basic structure representing Length+Lexicographic (left-to-right) ordering of
the words over given Alphabet. Lexicographic ordering of an Alphabet is
implicitly specified inside Alphabet struct.
"""
struct LenLex{T} <: WordOrdering
    A::Alphabet{T}
end

Base.hash(o::LenLex, h::UInt) = hash(o.A, hash(h, hash(LenLex)))

"""
    lt(o::LenLex, p::AbstractWord, q::AbstractWord)

Return whether the first word is less then the other one in a given LenLex ordering.
"""
function lt(o::LenLex, p::AbstractWord, q::AbstractWord)
    if length(p) == length(q)
        for (a, b) in zip(p,q)
            a == b || return isless(a, b)  # comparing only on positive pointer values
        end
        return false
    else
        return isless(length(p), length(q))
    end
end

"""
    struct WreathOrder{T} <: WordOrdering

Structure representing Basic Wreath-Product ordering (determined by the Lexicographic
ordering of the Alphabet) of the words over given Alphabet. This Lexicographic
ordering of an Alphabet is implicitly specified inside Alphabet struct.
"""
struct WreathOrder{T} <: WordOrdering
    A::Alphabet{T}
end

Base.hash(o::WreathOrder, h::UInt) = hash(o.A, hash(h, hash(WreathOrder)))

"""
    lt(o::WreathOrder, p::AbstractWord, q::AbstractWord)

Return whether the first word is less then the other one in a given WreathOrder ordering.
"""
function lt(o::WreathOrder, p::AbstractWord, q::AbstractWord)
    iprefix = lcp(p, q) + 1

    @views pp = p[iprefix:end]
    @views qq = q[iprefix:end]

    isone(pp) && return !isone(qq)
    isone(qq) && return false

    # i.e. pp, qq ≂̸ ε
    max_p = maximum(pp)
    max_q = maximum(qq)

    max_p < max_q && return true
    max_p > max_q && return false

    # i.e. max_p == max_q
    head_p_len = count(isequal(max_p), pp)
    head_q_len = count(isequal(max_p), qq)

    head_p_len < head_q_len && return true
    head_p_len > head_q_len && return false

    # i.e. head_p_len == head_q_len
    first_pp = findfirst(isequal(max_p), pp)
    first_qq = findfirst(isequal(max_p), qq)
    return @views lt(o, pp[1:first_pp - 1], qq[1:first_qq - 1])
end

"""
    struct RecursivePathOrder{T} <: WordOrdering

Structure representing Recursive Path Ordering (determined by the Lexicographic
ordering of the Alphabet) of the words over given Alphabet. This Lexicographic
ordering of an Alphabet is implicitly specified inside Alphabet struct.
For the definition see
> Susan M. Hermiller, Rewriting systems for Coxeter groups
> _Journal of Pure and Applied Algebra_
> Volume 92, Issue 2, 7 March 1994, Pages 137-148.
"""
struct RecursivePathOrder{T} <: WordOrdering
    A::Alphabet{T}
end

Base.hash(o::RecursivePathOrder, h::UInt) = hash(o.A, hash(h, hash(RecursivePathOrder)))

"""
    lt(o::RecursivePathOrder, p::AbstractWord, q::AbstractWord)

Return whether the first word is less then the other one in a given
RecursivePathOrder ordering.
"""
function lt(o::RecursivePathOrder, p::AbstractWord, q::AbstractWord)
    isone(p) && return !isone(q)
    isone(q) && return false

    # i.e. p, q ≂̸ ε
    length(p) == 1 && length(q) == 1 && return (first(p) < first(q))

    # i.e. we are not comparing single letter words
    @views begin
        p == q[2:end] && return true
        lt(o, p, q[2:end]) && return true
        first(p) < first(q) && lt(o, p[2:end], q) && return true
        first(p) == first(q) && lt(o, p[2:end], q[2:end]) && return true
    end
    return false
end


text:
function setnumbeats!(system::LegSystem,n::Int64)
    # total time needed for all past/current cardiac cycles
    ttotal = system.solverparams.th*(system.solverparams.numbeats+1);

    # increment number of cycles if new cycle is starting
    if system.t[n+2] >= ttotal
        println("Old number of completed cardiac cycles: $(system.solverparams.numbeats)")
        println("New number of completed cardiac cycles: $(system.solverparams.numbeats+1)")
        println("Total number of cycles to complete: $(system.solverparams.numbeatstotal)")
        system.solverparams.numbeats+=1;
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
using jInv.Mesh
using Test

# setup 3D problem
nc = [6, 7, 8]
x0 = rand(3)
domain = [x0[1], 4, x0[2], 2, x0[3], 6]
h   = (domain[2:2:end]-domain[1:2:end])./nc
h1  = h[1]*ones(nc[1])
h2  = h[2]*ones(nc[2])
h3  = h[3]*ones(nc[3])

Mt            = getTensorMesh3D(h1,h2,h3,x0)
l             = diag(getLength(Mt))
nex           = Mt.ne[1]
ne            = sum(Mt.ne)
Xe1,Xe2,Xe3   = getEdgeGrids(h1,h2,h3)
Xe1           = Xe1 .+ x0'

xt1 = x0[1] +   h[1]
xt2 = x0[1] + 5*h[1]
yt1 = x0[2] + 3*h[2]
yt2 = x0[2] + 4*h[2]
zt1 = x0[3] + 4*h[3]
zt2 = x0[3] + 5*h[3]

# test polygonal path directly on edges
idx     = findall( (xt1 .< Xe1[:,1] .< xt2) .& (Xe1[:,2] .== yt1) .& (Xe1[:,3] .== zt1) )
s1      = zeros(ne)
s1[idx] = l[idx]
p1      = [xt1 yt1 zt1
           xt2 yt1 zt1]
s1Comp  = getEdgeIntegralOfPolygonalChain(Mt,p1)

@test norm(s1-s1Comp) < 1e-10

s1Norm = getEdgeIntegralOfPolygonalChain(Mt,p1,normalize=true)
@test abs(sum(s1Norm) - 1.0) < 1e-10

# test path offset by half a cell in the y direction
idx1     = findall( (xt1 .< Xe1[:,1] .< xt2) .& (Xe1[:,2] .== yt1) .& (Xe1[:,3] .== zt1) )
idx2     = findall( (xt1 .< Xe1[:,1] .< xt2) .& (Xe1[:,2] .== yt2) .& (Xe1[:,3] .== zt1) )
s2       = zeros(ne)
s2[idx1] = l[idx1]./2
s2[idx2] = l[idx2]./2
p2       = [xt1 (yt1+yt2)/2 zt1
            xt2 (yt1+yt2)/2 zt1]

s2Comp   = getEdgeIntegralOfPolygonalChain(Mt,p2)
@test norm(s2-s2Comp) < 1e-10

# test path offset by half a cell in the y and z directions
idx3     = findall( (xt1 .< Xe1[:,1] .< xt2) .& (Xe1[:,2] .== yt1) .& (Xe1[:,3] .== zt2) )
idx4     = findall( (xt1 .< Xe1[:,1] .< xt2) .& (Xe1[:,2] .== yt2) .& (Xe1[:,3] .== zt2) )
s3       = zeros(ne)
s3[idx1] = l[idx1]./4
s3[idx2] = l[idx2]./4
s3[idx3] = l[idx3]./4
s3[idx4] = l[idx4]./4
p3       = [xt1 (yt1+yt2)/2 (zt1 + zt2)/2
            xt2 (yt1+yt2)/2 (zt1 + zt2)/2]

s3Comp   = getEdgeIntegralOfPolygonalChain(Mt,p3)
@test norm(s3-s3Comp) < 1e-10


text:
using Kalman
using Base.Test
using StaticArrays
using LaTeXStrings
include("testsystem1.jl")
srand(11)
using Plots

Y, X = sample(20, 1, M)
Xhat = kalmanrts(Y, M)[2]
scatter(Y[:], label=L"Y (= X_1 + \epsilon = obs)")
plot!(reshape(X,2,20)', color = [:red :blue], label=[L"X_1" L"X_2"] )# label=[L"X_1", L"X_2" ])
plot!(reshape(Xhat,2,20)', line= :dot, color = [:red :blue], label=[L"X_1 (est.)" L"X_2 (est.)"] )# label=[L"X_1", L"X_2" ])

savefig("test/testplot.png")


text:
##############################################################################
##
## Type with stored cholesky
##
##############################################################################

type DenseCholeskyAllocatedSolver{Tc <: StridedMatrix} <: AbstractAllocatedSolver
    chol::Tc
    function DenseCholeskyAllocatedSolver(chol)
        size(chol, 1) == size(chol, 2) || throw(DimensionMismatch("chol must be square"))
        new(chol)
    end
end

DenseCholeskyAllocatedSolver{Tc}(chol::Tc) = DenseCholeskyAllocatedSolver{Tc}(chol)
function AbstractAllocatedSolver{Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer, ::Cholesky)
    return DenseCholeskyAllocatedSolver(Array(eltype(nls.J), length(nls.x), length(nls.x)))
end

##############################################################################
##
## solve J'J \ J'y by Cholesky
##
##############################################################################

function A_ldiv_B!(x::AbstractVector, J::StridedMatrix, y::AbstractVector, A::DenseCholeskyAllocatedSolver)
    chol = A.chol
    Ac_mul_B!(chol, J,  J)
    Ac_mul_B!(x, J,  y)
    A_ldiv_B!(cholfact!(chol, :U, Val{true}), x)
    return x, 1
end

##############################################################################
##
## solve (J'J + diagm(damp)) \ J'y by Cholesky
##
##############################################################################

function A_ldiv_B!(x::AbstractVector, J::StridedMatrix, y::AbstractVector,
            damp::AbstractVector, A::DenseCholeskyAllocatedSolver)
    chol = A.chol

    # update chol as J'J + λdtd
    Ac_mul_B!(chol, J, J)
    # transform dammp
    size(chol, 1) == length(damp) || throw(DimensionMismatch("size(chol, 1) should equal length(damp)"))
   for i in 1:size(chol, 1)
        chol[i, i] += damp[i]
    end

    # solve
    Ac_mul_B!(x, J, y)
    A_ldiv_B!(cholfact!(chol), x)
    return x, 1
end


text:
"""
            Arrays with particle specific fields.
            Needed if no info block present.

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

gas_arr = ["U", "RHO", "HSML", "BFLD", "DIVB", "ABVC", "NE", "NH", "SFR",
             "VRMS", "VBLK", "Z", "ACVC"]

bh_arr = ["BHMA", "BHMD", "BHPC", "AGE", "ACRB"]

gas_bh = "TNGB"

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
"""
    initial_state(instance)

Initial state for a given instance (used for root node).
"""
function initial_state end

"""
    domain_type(instance)

The type used as domain for the decision variables.
"""
function domain_type end

"""
    value_type(instance)

The (numeric) type used in the objective function.
"""
function value_type end

"""
    next_variable(instance, diagram::Diagram{S,D,V}, variable_order)::Int

The variable to build the next layer, as an integer index.

Optional: Defaults to order `1:length(instance)`.

Multiple strategies can be implemented for an instance type by defining and
passing objects as `variable_order`.

The (work-in-progress) diagram is also passed and can be queried about already
assigned variables or the nodes & states of the last layer.
"""
function next_variable end

"""
    transitions(instance, state, variable::Int)::Dict{Arc{S,D,V},S}

All feasible transitions from the given state by any assignment in the
variable's domain.

Arcs specify the original state, the variable assignment and the contribution
the objective and are mapped to the new state.
"""
function transitions end

"""
    process(processing, layer::Layer{S,D,V})::Layer{S,D,V}

Build a new layer by processing the nodes of the given layer.

For restrictions, this could mean simply selecting a subset of the nodes, to
stay within the width limit.

For relaxations, new nodes are created by merging existing nodes (and setting
`exact=false`).

Optional: Defaults to `identity`.
"""
function process end


text:
"""
    NAGARCH{p, q, T<:AbstractFloat} <: UnivariateVolatilitySpec{T}
"""
struct NAGARCH{p, q, T<:AbstractFloat} <: UnivariateVolatilitySpec{T}
    coefs::Vector{T}
    function NAGARCH{p, q, T}(coefs::Vector{T}) where {p, q, T}
        length(coefs) == nparams(NAGARCH{p, q})  || throw(NumParamError(nparams(NAGARCH{p, q}), length(coefs)))
        new{p, q, T}(coefs)
    end
end

NAGARCH{p, q}(coefs::Vector{T}) where {p, q, T}  = NAGARCH{p, q, T}(coefs)

@inline nparams(::Type{<:NAGARCH{p, q}}) where {p, q} = p+2q+1
@inline nparams(::Type{<:NAGARCH{p, q}}, subset) where {p, q} = isempty(subset) ? 1 : sum(subset) + 1

@inline presample(::Type{<:NAGARCH{p, q}}) where {p, q} = max(p, q)

Base.@propagate_inbounds @inline function update!(
            ht, lht, zt, at, ::Type{<:NAGARCH{p, q}}, garchcoefs
            ) where {p, q}
	mht = garchcoefs[1]
    for i = 1:p
        mht += garchcoefs[i+1]*ht[end-i+1]
    end
    for i = 1:q
        mht += garchcoefs[i+1+p] * (at[end-i+1] + garchcoefs[i+1+p+q] * sqrt(ht[end-i+1]))^2
    end
    push!(ht, mht)
    push!(lht, (mht > 0) ? log(mht) : -mht)
    return nothing
end

# using ω / (1-α-αγ²-β)

@inline function uncond(::Type{<:NAGARCH{p, q}}, coefs::Vector{T}) where {p, q, T}
    den=one(T)
    for i = 1:p
        den -= coefs[i+1]
    end
    for i = 1:q
        den -= coefs[1+p+i]
    end
    for i = 1:q
        den -= coefs[1+p+i] * coefs[1+p+q+i]^2
    end
    h0 = coefs[1]/den
end

function startingvals(spec::Type{<:NAGARCH{p, q}}, data::Array{T}) where {p, q, T}
    x0 = zeros(T, p+2q+1)
    x0[1]=1
    x0[2:p+1] .= 0.9/p
    x0[p+2:end] .= 0.025/q
    x0[1] = var(data)/uncond(spec, x0)
    return x0
end

# using ω / (1-α-αγ²-β)

function startingvals(TT::Type{<:NAGARCH}, data::Array{T} , subset::Tuple) where {T}
	p, q = subsettuple(TT, subsetmask(TT, subset))
	x0 = zeros(T, p+2*q+1)
    x0[2:p+1] .= 0.9/p
    x0[p+2:end] .= 0.025/q
    x0[1] = var(data)*(one(T)-sum(x0[2:p+q+1])- sum(x0[p+q+2:end].^2 .* x0[p+2:p+q+1]))
	mask = subsetmask(TT, subset)
	x0long = zeros(T, length(mask))
	x0long[mask] .= x0
    return x0long
end

function constraints(::Type{<: NAGARCH{p,q}}, ::Type{T}) where {p, q, T}
    lower = zeros(T, p+2q+1)
    upper = ones(T, p+2q+1)
    upper[1] = T(Inf)
	return lower, upper
end

function coefnames(::Type{<:NAGARCH{p,q}}) where {p, q}
    names = Array{String, 1}(undef, p+2q+1)
    names[1] = "ω"
    names[2:p+1] .= (i -> "β"*subscript(i)).([1:p...])
    names[2+p:p+q+1] .= (i -> "α"*subscript(i)).([1:q...])
    names[p+q+2:p+2*q+1] .= (i -> "γ"*subscript(i)).([1:q...])
    return names
end

@inline function subsetmask(VS_large::Union{Type{NAGARCH{p, q}}, Type{NAGARCH{p, q, T}}}, subs) where {p, q, T}
	ind = falses(nparams(VS_large))
	subset = zeros(Int, 2)
	subset[3-length(subs):end] .= subs
	ind[1] = true
	ps = subset[1]
	qs = subset[2]
	@assert ps <= p
	@assert qs <= q
	ind[2:2+ps-1] .= true
	ind[2+p:2+p+qs-1] .= true
	ind[2+p+q:2+p+qs+q-1] .= true
	ind
end

@inline function subsettuple(VS_large::Union{Type{NAGARCH{p, q}}, Type{NAGARCH{p, q, T}}}, subsetmask) where {p, q, T}
	ps = 0
	qs = 0
	@inbounds @simd ivdep for i = 2 : p + 1
		ps += subsetmask[i]
	end
	@inbounds @simd ivdep for i = p + 2 : p + q + 1
		qs += subsetmask[i]
	end
	(ps, qs)
end


text:
"""
results(FastNet)

Return a refernce to the results of *sim* as a DataFrame

# Example
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[]);

julia> makenodes!(net,12,1)

julia> makenodes!(net,34,2)

julia> sim=FastSim(net,(r,t)->nothing,[])
Simulation run, currently at time 0.0

julia> runsim!(sim,100,25)
Time     Node state 1    Node state 2
  0.0              12              34
 25.0              12              34
 50.0              12              34
 75.0              12              34
100.0              12              34

julia> results(sim)
5×3 DataFrame
 Row │ Time     Node state 1  Node state 2
     │ Float64  Int64         Int64
─────┼─────────────────────────────────────
   1 │     0.0            12            34
   2 │    25.0            12            34
   3 │    50.0            12            34
   4 │    75.0            12            34
   5 │   100.0            12            34
```
"""
function results(sim::FastSim)
    sim.results
end

"""
    listnodes(FastNet)
    listnodes(FastNet,state)

Return a vector of the IDs of all nodes in FastNet *net* or all nodes in state *state* in *net*.

The one-argument method of this function returns a Vector containing the IDs of all nodes that
are in the network. The two-argument method creates a vector of all nodes in state *state*.

See also [listneighbors](#Fastnet.listneighbors)

# Example
```jldoctest
julia> using Fastnet

julia> net=FastNet(10,10,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,3,1)

julia> makenodes!(net,2,2)

julia> listnodes(net)
5-element Vector{Int64}:
 1
 2
 3
 4
 5

julia> listnodes(net,2)
2-element Vector{Int64}:
 4
 5
```
"""
function listnodes(net::FastNet)
    n=countnodes_f(net)
    ret=Vector{Int}(undef,n)
    cur=1
    nc=net.C-1
    for c=1:nc
        @inbounds nclen=net.nclen[c]
        for r=1:nclen
            @inbounds ret[cur]=node(net,c,r)
            cur+=1
        end
    end
    ret
end

function listnodes(net::FastNet,state::Int)
    nclen=net.nclen[state]
    ret=Vector{Int}(undef,nclen)
    cur=1
    for r=1:nclen
        ret[cur]=node_f(net,state,r)
        cur+=1
    end
    ret
end

"""
    listneighbors(FastNet,nid)

Return a vector of the IDs of all nodes that are asjacent to node *nid* in FastNet *net*.

This function is comparatively slow as it needs to allocate the vector. In your *rates!*
and process functions it is preferable to iterate over the neighbors using firstlinkout,
firstlinkin, nexlinkout, nextlinkin.

See also [listnodes](#Fastnet.listnodes)

# Example
```jldoctest
julia> using Fastnet

julia> net=FastNet(10,10,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,3,1)

julia> makenodes!(net,2,2)

julia> listnodes(net)
5-element Vector{Int64}:
 1
 2
 3
 4
 5

 julia> listnodes(net,2)
 2-element Vector{Int64}:
 4
 5
```
"""
function listneighbors(net::FastNet,node)
    nid=0
    try
        nid=convert(Int,node)
    catch e
        throw(ArgumentError("listneighbors expects its second argument (node) to be an integer"))
    end
    ret=Int[]
    cur=net.nout[nid]
    while cur!=0
        push!(ret,net.kdst[cur])
        cur=net.knexts[cur]
    end
    cur=net.nin[nid]
    while cur!=0
        push!(ret,net.ksrc[cur])
        cur=net.knextd[cur]
    end
    ret
end

"""
    listlinks(net)

Return an array that contains the source and destination of all links in *net*.

The return value will be a two-dimensional array of dimensions (K,3), where
K is the number of links in the network. Each row corresponds to one link.
The contests of the array are

- First column: Identical to the linkid of the respective link
- Second column: Source node of the link
- Third column: Destination node of the link

Warning: Do not rely on the row index to be identical to the link ID

See also [savelinklist](#Fastnet.savelinklist)

# Example
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,5,1)

julia> makelink!(net,1,2)
1

julia> makelink!(net,1,3)
2

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,5,1)

julia> makelink!(net,1,2)
1

julia> makelink!(net,2,3)
2

julia> makelink!(net,3,4)
3

julia> makelink!(net,4,5)
4

julia> makelink!(net,5,1)
5

julia> listlinks(net)
5×3 Matrix{Int64}:
 1  1  2
 2  2  3
 3  3  4
 4  4  5
 5  5  1
```
"""
function listlinks(net::FastNet)
    K=countlinks(net)
    res=Array{Int,2}(undef,K,3)
    for i=1:K
        kid=link_f(net,i)
        res[i,1]=kid
        res[i,2]=linksrc_f(net,kid)
        res[i,3]=linkdst_f(net,kid)
    end
    return res
end

"""
    listnodestates(net)

Return an array that contains the ids and states of all nodes in *net*.

The return value will be a two-dimensional array of dimension (K,2), where
K is the number of nodes in the network. Each row corresponds to one node.
The contests of the array ar

- First Column: ID of the respective node
- Second column: State of the node

Warning: Do not rely on the row index to be identical to the node ID

See also [listlinks](#Fastnet.listlinks)

# Example
```jldoctest
julia> using Fastnet

julia> net=FastNet(100,100,3,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,2,1)

julia> makenodes!(net,2,2)

julia> makenodes!(net,3,3)

julia> listnodestates(net)
7×2 Matrix{Int64}:
 1  1
 2  1
 3  2
 4  2
 5  3
 6  3
 7  3
```
"""
function listnodestates(net::FastNet)
    N=countnodes(net)
    res=Array{Int,2}(undef,N,2)
    for i=1:N
        nid=node_f(net,i)
        res[i,1]=nid
        res[i,2]=nodestate_f(net,nid)
    end
    res
end

"""
    savelinklist(net,filename)

Save network link list of network *net* to file *filename*.

The file is written in text mode. Each line that is written correosponds to
one link. The lines have the form

    LINKID SOURCE DESTINATION

Where LINKID is the respective link, SOURCE is the node id of the source node and DESTINATION
is the node id  of the destination node. The elemnts are separated by space. The line is ended by a line feed '\\n'.

See also [linklist](#Fastnet.linklist)
"""
function savelinklist(net::FastNet,filename::String)
    fl=0
    try
        fl=open(filename,"w+")
    catch e
        throw(ArgumentError("savelinklist was unable to open '$filename' for writing"))
    end
    K=countlinks(net)
    for i=1:K
        kid=link_f(net,i)
        src=linksrc_f(net,kid)
        dst=linkdst_f(net,kid)
        print(fl,"$kid $src $dst \n")
    end
    close(fl)
    nothing
end

"""
    savenodeinfo(net,filename)

Save information about the nodes to file *filename*.

The file is written in text mode. Each line that is written correosponds to
one nodes. The lines have the form

    NODEID STATE INDEGREE OUTDEGREE IN-NEIGHBORS OUT-NEIGBORS

Where NODEID is the ID of the node, STATE is the state of the node, in degree is the number of links that have the node as the
destination, OUTDEGREE is the number of links that have the node as the source, IN-NEIGHBORS is a list of all nodes from which the
focal node receives an incoming link, OUTNEIGHBORS is a list of the IDs of all nodes from which the focal node casts an outgoing link.
All elements are separated by spaces The line is ended by a line feed '\\n'.

See also [linklist](#Fastnet.linklist)
"""
function savenodeinfo(net::FastNet,filename)
    fl=0
    try
        fl=open(filename,"w+")
    catch e
        throw(ArgumentError("savenodeinfo was unable to open '$filename' for writing"))
    end
    N=countnodes(net)
    for i=1:N
        nid=node_f(net,i)
        state=nodestate_f(net,nid)
        indeg=indegree_f(net,nid)
        outdeg=outdegree_f(net,nid)
        print(fl,"$nid $state $indeg $outdeg")
        curlink=firstlinkin_f(net,nid)
        while curlink!=0
            src=linksrc_f(net,curlink)
            print(fl," $src")
            curlink=nextlinkin_f(net,curlink)
        end
        curlink=firstlinkout_f(net,nid)
        while curlink!=0
            dst=linkdst_f(net,curlink)
            print(fl," $dst")
            curlink=nextlinkout_f(net,curlink)
        end
        print(fl,"\n")
    end
    close(fl)
    nothing
end

"""
    degreedist(net::FastNet)

Return a vector of Float64 that specifies the networks degree distribution.

The element k of the returned vector is the probability that the probability that the a randomly drawn node from the
network has degree k.

If the elements of the vector do not add up to 1.0, the remainder is the probability that a randomly drawn node has
degree zero.

# Example
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,4,1)

julia> makelink!(net,node(net,1),node(net,2));

julia> makelink!(net,node(net,2),node(net,3));

julia> degreedist(net)
2-element Vector{Float64}:
 0.5
 0.25
```

"""
function degreedist(net::FastNet)
    mxdegree=0
    N=countnodes_f(net)
    for i=1:N
        mxdegree=max(mxdegree,degree_f(net,node_f(net,i)))
    end
    dd=zeros(Float64,mxdegree)
    for i=1:N
        k=degree_f(net,node_f(net,i))
        if k>0
          dd[k]+=1.0
        end
    end
    for i=1:mxdegree
        dd[i]/=N
    end
    dd
end

"""
    quickplotnodes(sim)

Return a UnicodePlot of the node timeseries stored in FastSim *sim*.

"""
function quickplotnodes(sim::FastSim)
    res=sim.results
    plt=lineplot(res[!,1],res[!,2],ylim=[0,sim.net.N],name=sim.net.nodealias[1],width=60,xlabel="Time t")
    for i=2:sim.net.C-1
        plt=lineplot!(plt,res[!,1],res[!,i+1],name=sim.net.nodealias[i])
    end
    plt
end

"""
    quickplotlink(sim)

Return a UnicodePlot of the link timeseries from the results stored in FastSim *sim*.

"""
function quickplotlinks(sim::FastSim)
    res=sim.results
    plt=lineplot(res[!,1],res[!,sim.net.C+1],ylim=[0,sim.net.K],name=sim.net.linkalias[1],width=60,xlabel="Time t")
    for i=2:sim.net.L-2
        plt=lineplot!(plt,res[!,1],res[!,sim.net.C+i],name=sim.net.linkalias[i])
    end
    plt
end


text:
"""
    load_json_model(filename::String)::JSONModel

Load and return a JSON-formatted model that is stored in `file_name`.
"""
function load_json_model(filename::String)::JSONModel
    return JSONModel(JSON.parsefile(filename))
end

"""
    save_json_model(model::MetabolicModel, file_name::String)

Save a [`JSONModel`](@ref) in `model` to a JSON file `file_name`.

In case the `model` is not `JSONModel`, it will be converted automatically.
"""
function save_json_model(model::MetabolicModel, file_name::String)
    m =
        typeof(model) == JSONModel ? model :
        begin
            @_io_log @warn "Automatically converting $(typeof(model)) to JSONModel for saving, information may be lost."
            convert(JSONModel, model)
        end

    JSON.print(open(file_name, "w"), m.json)
end


text:
@testset "537.complex-number-multiplication.jl" begin
    @test complex_number_multiply("1+1i", "1+1i") == "0+2i"
    @test complex_number_multiply("1+-1i", "1+-1i") == "0+-2i"
end