text:
###########################################################
## Canonical ordering of elements in orderless Mxpr       #
###########################################################

### !! This file has other canonicalizing code.: for Power

# There are four (more or less) steps:
# 1. Reduce sequences of numbers, leaving perhaps singletons. (We could choose to get singletons as well)
#    Doing this first is much,much faster for sums and prods with lots of numbers. May slow down
#    canoning other expressions.
# 2. Sort the terms and factors into a canonical order.
# 3. Evaluate all + and * between singleton numbers from 1., which have been sorted to beginning.
# 4. Combine
#    A. terms with same numeric coefficient
#    B. terms with same numeric power

## Comparison functions for sorting expressions to canonical order.

# Mathematica and Maxima sort terms and factors before evaluating-- it
# is difficult or impossible to hold an expression in non-canonical
# form (well, in Symata and Maxima, you have the source code.) Maple
# apparently does not sort expressions to canoncial order. But Maple
# (so they say) essentially never stores multiple copies of
# expressions that are the same: a hash key is computed for each
# expression, and the key is looked-up in a data structure. If the
# expression exists, a pointer to the original expression is used in
# place of the new one. In Maple, the hash key for for AC expressions
# is invariant under permutation of terms(factors).

# We try to follow the Mma order for Orderless Heads, only Plus and Times here.

# Anything with a Blank is greater than anything without a Blank.
# This is the natural order for pattern matching.
# Blank, BlankSequence, BlankNullSequence are not less than one another

#const _jstypeorder = Dict{DataType,Int}()
# In v0.6.0, Rational is of type UnionAll
const _jstypeorder = Dict{Any,Int}()

# orderless (commutative) functions will have terms ordered from first
# to last according to this order of types. Then lex within types.
# we'd like to put PyCall object in here, but we can't. Unless we generate the table at init time.
function _mklexorder()
    i = 1   ## CAREFULL! We use literal index of Any below
    for typ in (Float64,BigFloat,Int,BigInt,Any,Rational,Symbol,SJSym,Expr, Mxpr, Qsym)
        _jstypeorder[typ] = i
        i += 1
    end
    i = 1
end

_mklexorder()

# returns ordering precedence of Type, typ
function mxtypeorder(typ::DataType)
    ! haskey(_jstypeorder,typ)  && return 5  # Any
    return _jstypeorder[typ]
end

# Return the name of the blank or ""
# For sorting, no name becomes a zero-length string.
# The "blankname" is Blank[Name], this name is the type qualifier for the match
blankname(b) = length(b) == 0 ? "" : b[1]

# Order of Blankxxx's when they have no name (i.e. do not capture a match)
lessblanknoname(x::Mxpr{:Blank},y::Mxpr{:BlankSequence}) = true
lessblanknoname(x::Mxpr{:BlankSequence},y::Mxpr{:Blank}) = false
lessblanknoname(x::Mxpr{:Blank},y::Mxpr{:BlankNullSequence}) = true
lessblanknoname(x::Mxpr{:BlankNullSequence},y::Mxpr{:Blank}) = false
lessblanknoname(x::Mxpr{:BlankNullSequence},y::Mxpr{:BlankSequence}) = true
lessblanknoname(x::Mxpr{:BlankSequence},y::Mxpr{:BlankNullSequence}) = false

jslexless(x::Mxpr{:Pattern}, y::Mxpr{:Pattern}) = jslexless(x[1],y[1])

for b in ("Blank","BlankSequence","BlankNullSequence")
    for b1 in ("Blank","BlankSequence","BlankNullSequence")
        @eval begin
            function jslexless(x::Mxpr{Symbol($b)},y::Mxpr{Symbol($b1)})
                blankname(x) == blankname(y) == "" && return lessblanknoname(x,y)
                jslexless(blankname(x),blankname(y))
            end
        end
    end
    @eval begin
        jslexless(x::Mxpr{Symbol($b)},y::Mxpr{:Pattern}) = true
        jslexless(x::Mxpr{:Pattern},y::Mxpr{Symbol($b)}) = false
    end
end

# This works, but we did not need it
# is_nameless_blankxxx(mx::BlankXXX) = blankname(mx) == ""
# # Return true if mx contains a BlankXXX at any depth
# function expr_has_blank(mx::Mxpr)
#     is_blankxxx(mx) && return true
#     for i in 1:length(mx)
#         expr_has_blank(mx[i]) && return true
#     end
#     return false
# end
# expr_has_blank(x) = false

for b in ("Blank","BlankSequence","BlankNullSequence","Pattern")
    for o in ("Times","Plus","Power")
        @eval begin
            jslexless(x::Mxpr{Symbol($o)},y::Mxpr{Symbol($b)}) = false
            jslexless(x::Mxpr{Symbol($b)},y::Mxpr{Symbol($o)}) = true
        end
    end
    @eval begin
        jslexless(x::Mxpr{Symbol($b)},y::SJSym) = false
        jslexless(x::SJSym, y::Mxpr{Symbol($b)}) = true
    end
end

# Try using Union type BlankXXX here. No, this causes ambiguity errors.
for b in ("Blank","BlankSequence","BlankNullSequence")
    @eval begin
        function jslexless(x::Mxpr, y::Mxpr{Symbol($b)})
            return false
        end
        function  jslexless(x::Mxpr{Symbol($b)},y::Mxpr)
            return true
        end
    end
end

jslexless(x::Mxpr, y::Mxpr{:Pattern}) = true
jslexless(x::Mxpr{:Pattern}, y::Mxpr) = false

function jslexless(x::Mxpr{:Power}, y::Mxpr{:Power})
    jslexless(base(x),base(y)) || (!jslexless(base(y),base(x))) &&  jslexless(expt(x),expt(y))
end
jslexless(x::Mxpr{:Times}, y::Mxpr{:Times}) = jslexless(x[end],y[end])
jslexless(x::Mxpr{:Times}, y::Mxpr{:Power}) = jslexless(x[end],y)
jslexless(x::Mxpr{:Power}, y::Mxpr{:Times}) = jslexless(x,y[end])

function jslexless(x::Mxpr{:Power}, y::Mxpr)
    if y == base(x)
        return isa(expt(x),Real) && expt(x) < 0
    end
    jslexless(base(x),y)
end
function jslexless(x::Mxpr, y::Mxpr{:Power})
    if x == base(y)
        return isa(expt(y),Real) && expt(y) > 1
    end
    jslexless(x,base(y))
end

function jslexless(x::Mxpr{:Power}, y::SJSym)
    if y == base(x)
        return isa(expt(x),Real) && expt(x) < 0
    end
    jslexless(base(x),y)
end

function jslexless(x::SJSym, y::Mxpr{:Power})
    if x == base(y)
        return isa(expt(y),Real) && expt(y) > 1
    end
    jslexless(x,base(y))
end

function islesscomp(x, y)
    return !(jslexless(y, x)) || x == y
end

# Looks like symbol is less than times if it occurs among
# the arguments to times.
# Last element in a sorted expression is most significant

jslexless(x::Mxpr{:Times}, y::SJSym) = jslexless(x[end],y)
jslexless(x::SJSym, y::Mxpr{:Times}) = islesscomp(x,y)
#jslexless(x::SJSym, y::Mxpr{:Times}) = jslexless(x,y[end])

jslexless(x::Mxpr{:Plus}, y::SJSym) = jslexless(x[end],y)
jslexless(x::SJSym, y::Mxpr{:Plus}) = jslexless(x,y[end])

jslexless(x::Mxpr{:Times}, y::Mxpr) = jslexless(x[end],y)
jslexless(x::Mxpr, y::Mxpr{:Times}) = jslexless(x,y[end])

jslexless(x::Qsym, y::Qsym) = (x.context == y.context) ? jslexless(x.name, y.name) : jslexless(x.context, y.context)
function jslexless(x::T, y::Qsym) where T
    true
end
function jslexless(x::Qsym, y::T) where T
    false
end

const istrcmp =  @static Sys.iswindows() ? :_strnicmp : :strcasecmp

# cmp(a::Symbol, b::Symbol) is in base/string.jl
# Mma does case insensitive ordering of symbols, we do too
casecmp(a::Symbol, b::Symbol) = round(Int,sign(ccall(istrcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}), a, b)))
# Same as isless(a,b) for symbols, but we do case insensitive comparison first
# *and* we invert order of upper to lower case.  we want  a < A.
jsisless(a::Symbol, b::Symbol) = (r = casecmp(a,b); r == 0 ? cmp(a,b) > 0 : r < 0)
jslexless(x::SJSym, y::SJSym) = jsisless(x,y)
jslexless(x::Number, y::SJSym) = true
jslexless(x::Number, y::Mxpr) = true
jslexless(x::Mxpr, y::Number) = false
jslexless(x::SJSym, y::Mxpr) = true
jslexless(x::Mxpr, y::SJSym) = false

# Assume args of x and y are sorted. Compare
# last (most significant) element in expressions,
# x[i], y[j] (i and j are in general different becuase
# x and y are of different lengths)
# If x[i] != y[j], return true or false.
# Otherwise repeat, comparing next-to-last elements, etc.
# If all terms compared are equal, the shorter of
# x and y is the least.
function jslexless(x::Mxpr{T},y::Mxpr{T}) where T
    x === y && return false
    ax = margs(x)
    ay = margs(y)
    lx = length(ax)
    ly = length(ay)
    ix = lx
    iy = ly
    n = min(lx,ly)
    @inbounds for i in 1:n
        jslexless(ax[ix],ay[iy]) && return true
        if jslexless(ay[iy],ax[ix])
            return false
        end
        ix -= 1
        iy -= 1
    end
    return lx < ly
end

jslexless(x::Mxpr{T},y::Mxpr{V}) where {T,V} = T < V
# Following methods will only be called on all non-Symbolic types. (not Symbol or Mxpr)
_jslexless(x::DataType, y::DataType) = x <: y
# FIXE: Delete after testing on v0.7
_jslexless(x::T,y::T) where {T} = symlexless(x,y)
#_jslexless(x::T,y::T) where {T} = isless(x, y)  # use Julia definitions
jslexless(x::T,y::T) where {T} =  _jslexless(x, y)
# We have defined an order for different types
jslexless(x::T,y::V) where {T,V} = mxtypeorder(T) < mxtypeorder(V)

## Step 2. Sort expression according to jslexless.

# Order the args in orderless Mxpr.
function orderexpr!(mx::Mxpr)
    ar = mx.args
    sort!(ar,lt=jslexless)
    mx
end

# Step 1. Reduce sequences of explicit numbers.
# Sum (multiply) sequence of numbers in expression before sorting.
# This only removes one consecutive run of numbers, but we repeat.
for (op,name,id) in  ((:Plus,:plusfirst!,0),(:Times,:mulfirst!,1))
    @eval begin
        function ($name)(mx::Mxpr, n0::Int)
            args = margs(mx)
            len = length(args)
            n = 0
            @inbounds for i in n0:len
                if isa(args[i],Number)
                    n = i
                    break
                end
            end
            n == 0 && return (mx,len)
            s = $(op == :Plus ? :(zero(args[n]))  :  :(one(args[n])))
            m = 0
            @inbounds for i in n:len
                x = args[i]
                if isa(x,Number)
                    s = $(op == :Plus ? :(mplus(s,x)) : :(mmul(s,x)) )
                else
                    m = i
                    break
                end
            end
            m == 0 && n == 1 &&  return (s,len)
            if m == 0 m = len + 1 end
            splice!(args,n:m-1,[s])
            return (mx,m)
        end
    end
end

# Reduce all sequences of numbers in sums and products before sorting.
# Summing numbers in sums before sorting makes a large (eg 1000x) difference
# for large sums with many (or all) numbers.
# Singleton sequences may remain and they are sorted and removed with
# compactsumlike!. Probably could analyze sum and choose a strategy.
function loopnumsfirst!(mx::Orderless)
    len = length(mx)
    m = 1
    while true
        (mx,m) =  numsfirst!(mx,m)
        (isa(mx,Number) || m >= len) && return mx
        nlen = length(mx)
        p = len - m
        m = nlen - p
        len = nlen
    end
end

numsfirst!(mx::Mxpr{:Times},n) = mulfirst!(mx,n)
numsfirst!(mx::Mxpr{:Plus},n) = plusfirst!(mx,n)

# Try: Assume we can mutate and reuse the same Mxpr.
# (x^n)^m --> x^(n*m) , only for real n,m
function mulpowers(mx::Mxpr{:Power})
    (e,b) = numeric_expt_and_base(mx)
    (e1,b1) = numeric_expt_and_base(base(mx))
    e != 1 && e1 != 1 && return mxpr(:Power,b1,mmul(e,e1))
    return mx
end
mulpowers(x) = x

## Distribute -1 * sum

function distribute_minus_one_b(n,terms)
    newterms = newargs(length(terms))
    @inbounds for i in 1:length(terms)
        newterms[i] = mxpr(:Times,n,terms[i])
    end
    return newterms
end
function distribute_minus_one_a(mx, n::Integer, sum0::Mxpr{:Plus})
    n == -1 || return mx
    return mxpr(:Plus,distribute_minus_one_b(n,margs(sum0)))
end
distribute_minus_one_a(mx, n, m) = mx

function distribute_minus_one(mx::Mxpr{:Times})
    length(mx) != 2 && return mx
    distribute_minus_one_a(mx,mx[1],mx[2])
end
distribute_minus_one(x) = x

#### canonexpr!

struct Times_analysis
    gotFloat::Bool
    isZero::Bool
    isIndeterminate::Bool
    isInfinity::Bool
    isComplexInfinity::Bool
end

# For Plus we don't need to search for zero. Fix this.
function analyze_operands(mx::Orderless)
    ma = margs(mx)
    gotFloat::Bool = false
    gotZero::Bool = false
    gotIndeterminate::Bool = false
    gotInfinity::Bool = false
    gotComplexInfinity::Bool = false
    for i in 1:length(ma)
        a = ma[i]
        if is_Indeterminate(a)
            gotIndeterminate = true
            return  Times_analysis(gotFloat,gotZero,gotIndeterminate,gotInfinity,gotComplexInfinity)
        end
        if ! gotInfinity && is_Infintity(a) gotInfinity = true
        elseif ! gotComplexInfinity && is_ComplexInfinity(a) gotComplexInfinity = true
        elseif ! gotFloat && isa(a,AbstractFloat) gotFloat = true
        elseif ! gotZero && a == 0 gotZero = true
        end
    end
    if gotZero
        if gotInfinity || gotComplexInfinity
            gotIndeterminate = true
            return  Times_analysis(false,false,gotIndeterminate,false,false)
        else
            isZero = true
            return Times_analysis(gotFloat,isZero,false,false,false)
        end
    elseif gotInfinity
        if gotComplexInfinity
            isComplexInfinity = true
            return Times_analysis(false,false,false,false,isComplexInfinity)
        else
            isInfinity = true
            return Times_analysis(false,false,false,isInfinity,false)
        end
    elseif gotComplexInfinity
        isComplexInfinity = true
        return Times_analysis(false,false,false,false,isComplexInfinity)
    end
    return Times_analysis(gotFloat,gotZero,gotIndeterminate,gotInfinity,gotComplexInfinity)
end

# TODO: What to do about overflow ?
# maybe do the julia way. caveat emptor
function whichinfinity(mx::Mxpr{:Times})
    ma = margs(mx)
    prod = BigInt(1)  # This won't be type stable :(
    nonum = newargs()
    for i in 1:length(ma)
        x = ma[i]
        if isa(x,Number)
            prod *= x
        elseif mhead(x) != :DirectedInfinity
            push!(nonum,x)
        end
    end
    prodsq = real(prod * conj(prod))
    fac = mpow(prodsq,-1//2)
    direction = doeval(prod*mpow(prodsq,-1//2)) # TODO: optimize
    if length(nonum) > 0
        mxpr(:Times, nonum..., mxpr(:DirectedInfinity, direction))
    else
        return mxprcf(:DirectedInfinity, direction)
    end
end

# This is not complete
function whichinfinity(mx::Mxpr{:Plus})
    return mxprcf(:DirectedInfinity, 1)
end

function canonexpr_orderless!(mx::Mxpr{:Times})
    length(mx) == 0 && return 1
    length(mx) == 1 && return mx[1]
    analysis = analyze_operands(mx)
    if analysis.isIndeterminate
        return Indeterminate
    end
    if analysis.isComplexInfinity return ComplexInfinity end
    if analysis.isInfinity return whichinfinity(mx) end
    if analysis.isZero return 0 end # Need to fix the type of zero
    if analysis.gotFloat orderless_convert_to_float!(mx) end
    _canonexpr_orderless!(mx)
end

function canonexpr_orderless!(mx::Mxpr{:Plus})
    length(mx) == 0 && return 0
    length(mx) == 1 && return mx[1]
    analysis = analyze_operands(mx)
    if analysis.isIndeterminate return Indeterminate end
    if analysis.isComplexInfinity return ComplexInfinity end
    if analysis.isInfinity return whichinfinity(mx) end
    if analysis.gotFloat orderless_convert_to_float!(mx) end
    _canonexpr_orderless!(mx)
end

function orderless_convert_to_float!(mx)
    ma = margs(mx)
    for i in 1:length(ma)
        x = ma[i]
        if (! isa(x,Number)) && is_Numeric(x)
            ma[i] = do_N(x)
        end
    end
    nothing
end

### This is not a "toplevel" function.
## Apply all steps listed at top of this file.
# We say 'sum', but this applies to Times as well
function _canonexpr_orderless!(mx)
    @mdebug(2,"entering: _canonexpr_orderless! ", mx)
    mx = loopnumsfirst!(mx)  # remove sequences of numbers
    isa(mx,Number) && return mx
    mx = distribute_minus_one(mx)
    orderexpr!(mx)  # sort terms
    if isa(mx,Mxpr)
        mx = compactsumlike!(mx) # sum numbers not gotten by loopnumsfirst.
        @mdebug(2,"_canonexpr_orderless!, returned from compactsumlik1!: ", mx)
        if isa(mx,Mxpr)
            mx = collectordered!(mx)  # collect terms differing by numeric coefficients
            @mdebug(2,"_canonexpr_orderless!, returned from collectordered!: ", mx)
            # following is rarely if ever used.
            if isa(mx,PowerT)
                mx = mulpowers(mx)
            end  # add numeric exponents when base is same
        end
    end
    setcanon(mx)
    @mdebug(2,"exiting: _canonexpr_orderless! ", mx)
    mx
end

canonexpr!(mx::Orderless) = canonexpr_orderless!(mx)

function canonexpr!(mx::Mxpr)
    if get_attribute(mx,:Orderless)
        orderexpr!(mx)
    end
    mx
end

# We set fixed here, but not for orderless above. Which one is correct ?
function canonexpr!(mx::Mxpr{:Power})
    mx = do_canon_power!(mx,base(mx),expt(mx)) # Don't want "!" here
    if isa(mx, TimesT)
        mx = canonexpr!(mx)
    end
    setcanon(mx)
    # m = (a * b * c * d * e *f * g * h)^2  takes 2e-4 s if we don't do setfixed below
    # But, we can't do this: ex = 102^(1/2); ex^2 does not return 102 if we setfixed.
    if ! isa(mx, PowerT)
#        setfixed(mx)
    end
    return mx
end

# (expr1*expr2*....)^n --> expr1^n * expr2^n * .... for numeric n  *NO* expr1,.... must be positive as well
# Assume the product prod has already been sorted;
# the number, if present is first.
# Input is power of a product
function do_canon_power!(mx::PowerT,prod::TimesT, expt::Integer)
    len = length(prod)
    args = margs(prod)
    nargs = newargs(len) # we must use new args.
    if isa(args[1],Number)  # No more than the first item will be a number
        nargs[1] = mpow(args[1],expt)
    else
        nargs[1] = mxpr(:Power,args[1],expt)
    end
    @inbounds for i in 2:len
        nargs[i] = mxpr(:Power,args[i],expt)
    end
    @inbounds for i in 1:len
        setcanon(nargs[i])
    end
    mx = mxpra(:Times,nargs)
    return mx
end

do_canon_power!(mx,base,expt) = mx

## For expressions/objects that are not canonicalized
canonexpr!(x) = x

flatcanon!(x) = canonexpr!(maybeflatten(x))
#flatcanon!(x) = canonexpr!(x)  # These are for testing speed
#flatcanon!(x) = flatten(x)
#flatcanon!(x) = x

## Step 3.
# TODO  'compact' not a good name for this function
#################################################################
## Sum collected numerical args in :mplus, (or same for :Times)  #
#################################################################
# + and * are nary. Replace all numbers in the list of args, by one sum or product

compactsumlike!(mx::Mxpr{:Plus}) = compactplus!(mx)
compactsumlike!(mx::Mxpr{:Times}) = compactmul!(mx)

for (op,name,id) in  ((:Plus,:compactplus!,0),(:Times,:compactmul!,1))
    if op == :Times
        fop = :mmul
    else
        fop = :mplus
    end
    @eval begin
        function ($name)(mx::Mxpr)
            @mdebug(2,"In ", $name, "  ", mx)
            length(mx) < 2 && return mx
            @mdebug(3, $name, ": length(mx)>1 ", mx)
            a = margs(mx)
            @mdebug(3, $name, ": got margs: ", mx)
            isa(a[1],Number) || return mx
            sum0 = a[1]
            while length(a) > 1
                popfirst!(a)
                isa(a[1],Number) || break
                sum0 = ($fop)(sum0,a[1])
            end
            @mdebug(3, $name, ": done while loop, a=$a, sum0=$sum0")
            (length(a) == 0 || isa(a[1],Number)) && return sum0
            $(fop == :mmul ? :(sum0 == 0 && return handle_coefficient_is_zero(mx,sum0)) : :())
            sum0 != $id && pushfirst!(a, sum0)
            @mdebug(3, $name, ": checking length: ",a)
            length(a) == 1 && return a[1]
            @mdebug(3, $name, ": returning at end")
            return mx
        end
    end
end

# This is now handled earlier
function handle_coefficient_is_zero(mx,sum0)
    return sum0
end

# Get numeric coefficient of expr. It is 1 if there is none.
function numeric_coefficient(x::Mxpr{:Times})
    c = isa(x[1],Number) ? x[1] : 1
end
numeric_coefficient(x::Number) = x
numeric_coefficient(x) = 1

# Get numeric exponent; may be 1 (zero never encountered, I hope)
function numeric_expt(x::Mxpr{:Power})
    c = isa(expt(x),Number) ? expt(x) : 1
end
numeric_expt(x::Number) = 1
numeric_expt(x) = 1

# TODO: Tests fail unless we copy. But, there may be a way around this
function _rest!(mx::Mxpr)
    res=copy(mx)  # could be slow
    popfirst!(margs(res))
    return length(res) == 1 ? res[1] : res
end

# split product n*expr into (n,expr) with numeric n. n may be 1
function numeric_coefficient_and_factor(a)
    n = numeric_coefficient(a)
    res = n == 1 ? (n,a) : (n,_rest!(a))
    return res
end

function numeric_expt_and_base(a)
    n = numeric_expt(a)
    return n == 1 ? (n,a) : (n,base(a))
end

# Test if two terms differ only by numeric coeffcient,
# and return coefficient and common sub-expression.
# name is too generic.
function _matchterms(a,b)
    (na,a1) = numeric_coefficient_and_factor(a)
    (nb,b1) = numeric_coefficient_and_factor(b)
    return a1 == b1 ? (true,na+nb,a1) : (false,0,a1)
end

# Same as above, but for products
function _matchfacs(a,b)
    (na,a1) = numeric_expt_and_base(a)
    (nb,b1) = numeric_expt_and_base(b)
    return a1 == b1 ? (true,na+nb,a1) : (false,0,a1)
end

## Step 4.
# Replace n repeated terms expr by expr*x, and factors expr by expr^n
# This assumes that terms have been ordered so that collectable
# terms are neighboring.
# We now used newcollectmmul! below which collects exponents
# with m and or n symbolic (symbols or expressions)
# expr^n * expr^m  --> expr^(n+m)
collectordered!(x) = x
collectordered!(mx::Mxpr{:Plus}) = collectmplus!(mx)
collectordered!(mx::Mxpr{:Times}) = newcollectmmul!(mx)
for (op,name,matchf) in  ((:mplus,:collectmplus!, :_matchterms),
                          (:mmul,:collectmmul!,:_matchfacs))
    @eval begin
        function ($name)(mx::Mxpr)
            length(mx) < 2 && return mx
            a = margs(mx)
            n = 1
            count = 0
            coeffcount = 0
            while n < length(a)
                isa(a[n],Number) && (n += 1; continue)
                (success,coeffsum,fac) = ($matchf)(a[n],a[n+1])
                if success
                    count = 1
                    coeffcount = coeffsum
                    @inbounds for i in (n+1):(length(a)-1)
                        (success1,coeffsum1,fac1) = ($matchf)(fac,a[i+1])
                        #                        if a[i] == a[i+1]
                        if success1
                            count += 1
                            coeffcount += coeffsum1 - 1
                        else
                            break
                        end
                    end
                    newex = $(op== :mplus ?
                    :(coeffcount == 1 ? fac : mxpr(:Times,coeffcount,fac)) :
                    :(coeffcount == 0 ? 0 : coeffcount == 1 ? fac : mxpr(:Power,fac,coeffcount)))
                    if coeffcount == 0
                        splice!(a,n:n+count)
                    else
                        splice!(a,n:n+count,[newex])
                    end
                end
                n += 1
            end
            if  length(a) == 1
                return a[1]
            end
            if length(a) == 0
                $(op== :mplus ? :(return 0) : :(return 1))
            end
            return mx
        end
    end
end

function any_expt(x::Mxpr{:Power})
    expt(x)
end
any_expt(x::Number) = 1
any_expt(x) = 1

function any_expt_and_base(a)
    n = any_expt(a)
    return n == 1 ? (n,a) : (n,base(a))
end

# Same as above, but for products
function _matchfacs2a(a,b,cumterms)
    (na,a1) = any_expt_and_base(a)
    (nb,b1) = any_expt_and_base(b)
    if a1 == b1
        push!(cumterms,na)
        push!(cumterms,nb)
        return(true,a1)
    else
        return (false,a1)
    end
end

function _matchfacs2b(a,b,cumterms)
    (na,a1) = any_expt_and_base(a)
    (nb,b1) = any_expt_and_base(b)
    if a1 == b1
        push!(cumterms,nb)
        return(true,a1)
    else
        return (false,a1)
    end
end

function newcollectmmul!(mx::Mxpr)
    length(mx) < 2 && return mx
    a = margs(mx)
    n = 1
    count = 0
    coeffcount = 0
    cumterms = Array{Any}(undef, 0)
    while n < length(a)
        isa(a[n],Number) && (n += 1; continue)
        (success,fac) = _matchfacs2a(a[n],a[n+1],cumterms)
        if success
            count = 1
            @inbounds for i in (n+1):(length(a)-1)
                (success1,fac1) = _matchfacs2b(fac,a[i+1],cumterms)
                if success1
                    count += 1
                else
                    break
                end
            end
            newex =
                length(cumterms) == 0 ? 0 : length(cumterms) == 1 ? (cumterms[1] == 1 ? fac : mxpr(:Power,fac,cumterms[1])) :
                    mxpr(:Power,fac,mxpr(:Plus, cumterms...))
            splice!(a,n:n+count,[newex])
            empty!(cumterms)
        end
        n += 1
    end
    if  length(a) == 1
        return a[1]
    end
    if length(a) == 0
        return 1
    end
    return mx
end


text:
"""
    capm(asset_returns, benchmark_returns; risk_free=0.0)

Calculates the CAPM α and β coefficients based on sample covariance statistics and a simple linear regression model.

The linear regression model looks as follows:

        r_a - r_f = α + β(r_b - r_f) + ϵ

The α coefficient in this model is also known as Jensen's alpha. β is the slope coefficient, and ϵ is an error term.

# Arguments
- `asset_returns`:      Vector of asset returns.
- `benchmark_returns`:  Vector of benchmark returns (e.g. market portfolio returns for CAPM beta).
- `risk_free`:          Optional vector or scalar value denoting the risk-free return (must have same frequency as the provided returns, e.g. daily).

# Returns
Tuple (α, β) with the estimated α and β coefficients of the CAPM model.
"""
function capm(asset_returns, benchmark_returns; risk_free=0.0)
    asset_returns_ex = asset_returns .- risk_free
    benchmark_returns_ex = benchmark_returns .- risk_free
    μ1 = mean(asset_returns_ex)
    μ2 = mean(benchmark_returns_ex)
    β = sum((asset_returns_ex .- μ1).*(benchmark_returns_ex .- μ2)) / sum((benchmark_returns_ex .- μ2).^2)
    α = μ1 - β*μ2
    (α, β)
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

module Variography

using GeoStatsBase

using Printf
using Optim
using Distances
using Statistics
using SpecialFunctions: besselk, gamma
using StaticArrays
using RecipesBase
using Parameters

# extend methods for theoretical variograms
import Base: +

# variogram models
include("empirical.jl")
include("theoretical.jl")
include("pairwise.jl")
include("fitting.jl")

# plot recipes
include("plotrecipes/empirical.jl")
include("plotrecipes/theoretical.jl")
include("plotrecipes/hscatter.jl")
include("plotrecipes/varplane.jl")

export
  # empirical variograms
  EmpiricalVariogram,
  DirectionalVariogram,

  # theoretical variograms
  Variogram,
  GaussianVariogram,
  ExponentialVariogram,
  MaternVariogram,
  SphericalVariogram,
  CubicVariogram,
  PentasphericalVariogram,
  PowerVariogram,
  SineHoleVariogram,
  CompositeVariogram,
  isstationary,
  sill, nugget,
  pairwise,
  pairwise!,

  # fitting methods
  WeightedLeastSquares

end


text:
"""
    truncresponse(response, L, noisepower; warn=true)
Truncate the given array to the given length, printing a warning if the
truncated tail has more average power than the given noise floor.

`noisepower` can be a scalar or a row vector with one entry per channel

Warning can be disabled with the kwarg `warn=false`
"""
function truncresponse(response::AbstractArray, L, noisepower=eps(); warn=true)
    RL = size(response, 1)
    if RL > L
        truncenergy = sum(response[L+1:end, :].^2; dims=1)
        headroom = 1.05 # 5% over noise
        atol = sqrt(eps()) # absolute tolerance (for when noise floor = 0.0)
        if warn && any(truncenergy .> (atol + noisepower * headroom) * (RL-L))
            @warn string("nonsilent samples past end of analysis window. Check ",
                         "your test signal is long enough for the response")
            # Base.show_backtrace(STDERR, backtrace())
        end
        response = timeslice(response, 1:L)
    end

    response
end

# internal functions to do the actual truncation while keeping the
# dimensionality
timeslice(response::AbstractVector, I) = response[I]
timeslice(response::AbstractMatrix, I) = response[I, :]

# used to help parse arguments
striphz(f) = f
striphz(f::Frequency{T}) where T = T(f/Hz)
stripsec(t) = t
stripsec(t::Time{T}) where T = T(t/s)

energy(v) = sum(x->x^2, v)
power(v) = energy(v)/length(v)

"""
Normalize a frequency with a given samplerate to cycles/sample.
"""
freqnorm(f, sr) = uconvert(NoUnits, f/sr)

"""
Normalize a duration with a given samplerate to integer samples (rounded)
"""
durnorm(dur, sr) = round(Int, uconvert(NoUnits, dur*sr))

"""
    snr(mixture, signal)

Returns the signal-to-noise power ratio assuming that the mixture contains a
scaled version of the signal.
"""
function snr(mixture, signal)
    mixcoef = dot(mixture, signal) / dot(signal, signal)
    proj = mixcoef * signal

    energy(proj)/energy(mixture-proj)
end


text:
# IO routines for the Abaqus .inp file format
const ABAQUS_CPS3  = VTK_TRIANGLE
const ABAQUS_CPS4  = VTK_QUAD
const ABAQUS_CPS6  = VTK_QUADRATIC_TRIANGLE
const ABAQUS_CPS8  = VTK_QUADRATIC_QUAD
const ABAQUS_C3D4  = VTK_TETRA
const ABAQUS_C3D8  = VTK_HEXAHEDRON
const ABAQUS_C3D10 = VTK_QUADRATIC_TETRA
const ABAQUS_C3D20 = VTK_QUADRATIC_HEXAHEDRON

function _abaqus_element_string_to_int(x::String)
    if x === "CPS3"
        return ABAQUS_CPS3
    elseif x === "CPS4"
        return ABAQUS_CPS4
    elseif x === "CPS6"
        return ABAQUS_CPS6
    elseif x === "CPS8"
        return ABAQUS_CPS8
    elseif x === "C3D4"
        return ABAQUS_C3D4
    elseif x === "C3D8"
        return ABAQUS_C3D8
    elseif x === "C3D10"
        return ABAQUS_C3D10
    elseif x === "C3D20"
        return ABAQUS_C3D20
    else
        error("Invalid Abaqus element type.")
        return nothing
    end
end

function read_abaqus(path::String, ::Type{T}) where {T<:AbstractFloat}
    file = open(path, "r")
    try
        nodes = Point{3,T}[]
        element_types = UInt8[]
        offsets = UInt64[]
        elements = UInt64[]
        name = ""
        elsets = Dict{String, BitSet}()
        for line in eachline(file)
            if length(line) > 0
                if startswith(line, "**") # Comment
                    continue
                elseif "*Heading" === line
                    # Of the form " name.inp"
                    name = readline(file)[2:end-4]
                elseif "*NODE" === line
                    _read_abaqus_nodes!(file, nodes)
                elseif startswith(line, "*ELEMENT")
                    m = match(r"type=(.*?), ELSET", line)
                    if isnothing(m)
                        error("Incorrectly formatted Abaqus file?")
                    else
                        element_type = _abaqus_element_string_to_int(String(m.captures[1]))
                    end
                    _read_abaqus_elements!(file, element_type, element_types, offsets, elements)
                elseif startswith(line, "*ELSET")
                    m = match(r"(?<=ELSET=).*", line)
                    if isnothing(m)
                        error("Incorrectly formatted Abaqus file?")
                    else
                        set_name = m.match
                    end
                    elsets[set_name] = _read_abaqus_elset(file)
                end
            end
        end
        abaqus_2d = (ABAQUS_CPS3, ABAQUS_CPS4, ABAQUS_CPS6, ABAQUS_CPS8)
        abaqus_3d = (ABAQUS_C3D4, ABAQUS_C3D8, ABAQUS_C3D10, ABAQUS_C3D20)
        is2d = any(x->x ∈ abaqus_2d, element_types)
        is3d = any(x->x ∈ abaqus_3d, element_types)
        if is2d && is3d
            error("File contains both surface (CPS) and volume (C3D) elements."*
                  "Limit element types to be CPS or C3D, so mesh dimension may be determined.")
        end
        U = _select_uint_type(max(length(elements),
                                  length(nodes)))
        # Add the final offset
        push!(offsets, length(elements)+1)
        # Set materials if there are any
        materials = zeros(UInt8, length(element_types))
        material_names = String[]
        for key in keys(elsets)
            if startswith(key, "Material")
                push!(material_names, key)
            end
        end
        nmaterials = length(material_names)
        if 0 < nmaterials
            sort!(material_names)
            for (i, mat) in enumerate(material_names)
                for element_id in elsets[mat]
                    if materials[element_id] === UInt8(0)
                        materials[element_id] = i
                    else
                        error("Attempting to assign two materials to the same element.")
                    end
                end
                pop!(elsets, mat)
            end
        end
        trimmed_names = [mat[11:end] for mat in material_names]
        if is2d
            return VolumeMesh{2,T,U}(nodes, offsets, elements, element_types,
                                     materials, trimmed_names, name, elsets)
        else
            return VolumeMesh{3,T,U}(nodes, offsets, elements, element_types,
                                     materials, trimmed_names, name, elsets)
        end
    finally
        close(file)
    end
    return nothing
end

function _read_abaqus_nodes!(file::IOStream, nodes::Vector{Point{3,T}}) where {T<:AbstractFloat}
    # Count the number of nodes
    mark(file)
    nnodes = 0
    line = readline(file)
    while '*' !== line[1]
        nnodes += 1
        line = readline(file)
    end
    reset(file)
    # Allocate and populate a vector of nodes
    new_nodes = Vector{Point{3,T}}(undef, nnodes)
    for i in 1:nnodes
        line = readline(file)
        xyz = view(split(line, ','), 2:4)
        new_nodes[i] = Point{3,T}(ntuple(i->parse(Float64, xyz[i]), Val(3)))
    end
    append!(nodes, new_nodes)
    return nothing
end

function _read_abaqus_elements!(file::IOStream,
                                element_type::UInt8,
                                element_types::Vector{UInt8},
                                offsets::Vector{UInt64},
                                elements::Vector{UInt64})
    mark(file)
    line = readline(file)
    npts = points_in_vtk_type(element_type)
    nelements = 0
    while !('*' === line[1] || eof(file))
        nelements += 1
        line = readline(file)
    end
    reset(file)
    append!(element_types, fill(element_type, nelements))
    new_offsets = Vector{UInt64}(undef, nelements)
    new_elements = Vector{UInt64}(undef, npts*nelements)
    ectr = 1
    len_elements = length(elements)
    for i in 1:nelements
        line = readline(file)
        vids = view(split(line, ','), 2:npts+1)
        @. new_elements[ectr:ectr + npts - 1] = parse(UInt64, vids)
        new_offsets[i] = ectr + len_elements
        ectr += npts
    end
    append!(offsets, new_offsets)
    append!(elements, new_elements)
    return elements
end

function _read_abaqus_elset(file::IOStream)
    splitline = strip.(split(readuntil(file, "*")), [','])
    seek(file, position(file)-1)
    return BitSet(parse.(Int64, splitline))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Base: ImmutableDict

function binarize(e::T) where {T}
  !istree(e) && return e
  head = exprhead(e)
  if head == :call
    op = operation(e)
    args = arguments(e)
    meta = metadata(e)
    if op ∈ binarize_ops && arity(e) > 2
      return foldl((x, y) -> similarterm(e, op, [x, y], symtype(e); metadata = meta, exprhead = head), args)
    end
  end
  return e
end

"""
Recursive version of binarize
"""
function binarize_rec(e::T) where {T}
  !istree(e) && return e
  head = exprhead(e)
  op = operation(e)
  args = map(binarize_rec, arguments(e))
  meta = metadata(e)
  if head == :call
    if op ∈ binarize_ops && arity(e) > 2
      return foldl((x, y) -> similarterm(e, op, [x, y], symtype(e); metadata = meta, exprhead = head), args)
    end
  end
  return similarterm(e, op, args, symtype(e); metadata = meta, exprhead = head)
end

const binarize_ops = [:(+), :(*), (+), (*)]

function cleanast(e::Expr)
  # TODO better line removal
  if isexpr(e, :block)
    return Expr(e.head, filter(x -> !(x isa LineNumberNode), e.args)...)
  end

  # Binarize
  if isexpr(e, :call)
    op = e.args[1]
    if op ∈ binarize_ops && length(e.args) > 3
      return foldl((x, y) -> Expr(:call, op, x, y), @view e.args[2:end])
    end
  end
  return e
end

# Linked List interface
@inline assoc(d::ImmutableDict, k, v) = ImmutableDict(d, k => v)

struct LL{V}
  v::V
  i::Int
end

islist(x) = istree(x) || !isempty(x)

Base.empty(l::LL) = empty(l.v)
Base.isempty(l::LL) = l.i > length(l.v)

Base.length(l::LL) = length(l.v) - l.i + 1
@inline car(l::LL) = l.v[l.i]
@inline cdr(l::LL) = isempty(l) ? empty(l) : LL(l.v, l.i + 1)

# Base.length(t::Term) = length(arguments(t)) + 1 # PIRACY
# Base.isempty(t::Term) = false
# @inline car(t::Term) = operation(t)
# @inline cdr(t::Term) = arguments(t)

@inline car(v) = istree(v) ? operation(v) : first(v)
@inline function cdr(v)
  if istree(v)
    arguments(v)
  else
    islist(v) ? LL(v, 2) : error("asked cdr of empty")
  end
end

@inline take_n(ll::LL, n) = isempty(ll) || n == 0 ? empty(ll) : @views ll.v[(ll.i):(n + ll.i - 1)] # @views handles Tuple
@inline take_n(ll, n) = @views ll[1:n]

@inline function drop_n(ll, n)
  if n === 0
    return ll
  else
    istree(ll) ? drop_n(arguments(ll), n - 1) : drop_n(cdr(ll), n - 1)
  end
end
@inline drop_n(ll::Union{Tuple,AbstractArray}, n) = drop_n(LL(ll, 1), n)
@inline drop_n(ll::LL, n) = LL(ll.v, ll.i + n)

isliteral(::Type{T}) where {T} = x -> x isa T
is_literal_number(x) = isliteral(Number)(x)

# checking the type directly is faster than dynamic dispatch in type unstable code
_iszero(x) = x isa Number && iszero(x)
_isone(x) = x isa Number && isone(x)
_isinteger(x) = (x isa Number && isinteger(x)) || (x isa Symbolic && symtype(x) <: Integer)
_isreal(x) = (x isa Number && isreal(x)) || (x isa Symbolic && symtype(x) <: Real)

issortedₑ(args) = issorted(args, lt = <ₑ)
needs_sorting(f) = x -> is_operation(f)(x) && !issortedₑ(arguments(x))

# are there nested ⋆ terms?
function isnotflat(⋆)
  function (x)
    args = arguments(x)
    for t in args
      if istree(t) && operation(t) === (⋆)
        return true
      end
    end
    return false
  end
end

function hasrepeats(x)
  length(x) <= 1 && return false
  for i in 1:(length(x) - 1)
    if isequal(x[i], x[i + 1])
      return true
    end
  end
  return false
end

function merge_repeats(merge, xs)
  length(xs) <= 1 && return false
  merged = Any[]
  i = 1

  while i <= length(xs)
    l = 1
    for j in (i + 1):length(xs)
      if isequal(xs[i], xs[j])
        l += 1
      else
        break
      end
    end
    if l > 1
      push!(merged, merge(xs[i], l))
    else
      push!(merged, xs[i])
    end
    i += l
  end
  return merged
end

# Take a struct definition and make it be able to match in `@rule`
macro matchable(expr)
  @assert expr.head == :struct
  name = expr.args[2]
  if name isa Expr && name.head === :curly
    name = name.args[1]
  end
  fields = filter(x -> !(x isa LineNumberNode), expr.args[3].args)
  get_name(s::Symbol) = s
  get_name(e::Expr) = (@assert(e.head == :(::)); e.args[1])
  fields = map(get_name, fields)
  quote
    $expr
    TermInterface.istree(::$name) = true
    TermInterface.operation(::$name) = $name
    TermInterface.arguments(x::$name) = getfield.((x,), ($(QuoteNode.(fields)...),))
    TermInterface.arity(x::$name) = $(length(fields))
    Base.length(x::$name) = $(length(fields) + 1)
  end |> esc
end

using TimerOutputs

const being_timed = Ref{Bool}(false)

macro timer(name, expr)
  :(
    if being_timed[]
      @timeit $(esc(name)) $(esc(expr))
    else
      $(esc(expr))
    end
  )
end

macro iftimer(expr)
  esc(expr)
end

function timerewrite(f)
  reset_timer!()
  being_timed[] = true
  x = f()
  being_timed[] = false
  print_timer()
  println()
  x
end

"""
    @timerewrite expr

If `expr` calls `simplify` or a `RuleSet` object, track the amount of time
it spent on applying each rule and pretty print the timing.

This uses [TimerOutputs.jl](https://github.com/KristofferC/TimerOutputs.jl).

## Example:

```julia

julia> expr = foldr(*, rand([a,b,c,d], 100))
(a ^ 26) * (b ^ 30) * (c ^ 16) * (d ^ 28)

julia> @timerewrite simplify(expr)
 ────────────────────────────────────────────────────────────────────────────────────────────────
                                                         Time                   Allocations
                                                 ──────────────────────   ───────────────────────
                Tot / % measured:                     340ms / 15.3%           92.2MiB / 10.8%

 Section                                 ncalls     time   %tot     avg     alloc   %tot      avg
 ────────────────────────────────────────────────────────────────────────────────────────────────
 Rule((~y) ^ ~n * ~y => (~y) ^ (~n ...    667   11.1ms  21.3%  16.7μs   2.66MiB  26.8%  4.08KiB
   RHS                                       92    277μs  0.53%  3.01μs   14.4KiB  0.14%     160B
 Rule((~x) ^ ~n * (~x) ^ ~m => (~x)...    575   7.63ms  14.6%  13.3μs   1.83MiB  18.4%  3.26KiB
 (*)(~(~(x::!issortedₑ))) => sort_arg...    831   6.31ms  12.1%  7.59μs    738KiB  7.26%     910B
   RHS                                      164   3.03ms  5.81%  18.5μs    250KiB  2.46%  1.52KiB
   ...
   ...
 ────────────────────────────────────────────────────────────────────────────────────────────────
(a ^ 26) * (b ^ 30) * (c ^ 16) * (d ^ 28)
```
"""
macro timerewrite(expr)
  :(timerewrite(() -> $(esc(expr))))
end


text:
"""
    Hypergraph{T} <: AbstractMatrix{Union{T, Nothing}}

A hypergraph storing information about vertices and hyperedges.

**Constructors**

    Hypergraph{T}(n,k) where {T<:Real}
    Hypergraph{T,V}(n, k;
        v_meta=Vector{Union{V,Nothing}}(nothing, n)
        ) where {T<:Real, V}
    Hypergraph{T,V,E}(n, k;
        v_meta=Vector{Union{V,Nothing}}(nothing, n),
        he_meta=Vector{Union{E,Nothing}}(nothing, k)
        ) where {T<:Real, V, E}

Construct a hypergraph with a given number of vertices and hyperedges.
Optionally, values of type `V` can be stored at vertices and values of type `E`
can be stored at hyperedges.

    Hypergraph(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real}
    Hypergraph{V}(m::AbstractMatrix{Union{T, Nothing}};
        v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1))
        ) where {T<:Real, V}
    Hypergraph{V, E}(m::AbstractMatrix{Union{T, Nothing}};
        v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1)),
        he_meta=Vector{Union{E,Nothing}}(nothing, size(m,2))
        ) where {T<:Real, V, E}

Construct a hypergraph using its matrix representation.
In the matrix representation rows are vertices and columns are hyperedges.
Optionally, values of type `V` can be stored at vertices and values of type `E`
can be stored at hyperedges.

    Hypergraph(g::LightGraphs.Graph)

Constructs a hypergraph of degree 2 by making a deep copy of LightGraphs.Graph

**Arguments**

* `T` : type of weight values stored in the hypergraph
* `V` : type of values stored in the vertices of the hypergraph
* `E` : type of values stored in the edges of the hypergraph
* `n` : number of vertices
* `k` : number of hyperedges
* `m` : a matrix representation rows are vertices and columns are hyperedges
"""
struct Hypergraph{T,V,E} <: AbstractMatrix{Union{T, Nothing}}
    v2he::Vector{Dict{Int,T}}
    he2v::Vector{Dict{Int,T}}
    v_meta::Vector{Union{V,Nothing}}
    he_meta::Vector{Union{E,Nothing}}
    Hypergraph{T,V,E}(n, k,
            v_meta=Vector{Union{V,Nothing}}(nothing, n),
            he_meta=Vector{Union{E,Nothing}}(nothing, k)
            ) where {T<:Real, V, E} =
        new{T,V,E}(
            [Dict{Int,T}() for i in 1:n],[Dict{Int,T}() for i in 1:k],
            v_meta, he_meta)
end

Hypergraph{T}(n, k) where {T<:Real} =  Hypergraph{T,Nothing,Nothing}(n, k)

Hypergraph{T,V}(n, k;
              v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1))
             ) where {T<:Real, V} = Hypergraph{T,V,Nothing}(n, k; v_meta=v_meta)

function Hypergraph{V, E}(m::AbstractMatrix{Union{T, Nothing}};
                        v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1)),
                        he_meta=Vector{Union{E,Nothing}}(nothing, size(m,2))
                        ) where {T<:Real, V, E}
    n, k = size(m)
    h = Hypergraph{T, V, E}(n, k, v_meta, he_meta)
    h .= m
    h
end

Hypergraph{V}(m::AbstractMatrix{Union{T, Nothing}};
            v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1))) where {T<:Real, V} =
    Hypergraph{V, Nothing}(m; v_meta=v_meta)
Hypergraph(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real} =
    Hypergraph{Nothing, Nothing}(m)

function Hypergraph(g::LightGraphs.Graph)
    h = Hypergraph{Bool}(maximum(vertices(g)), ne(g))
    e = 0
    for edge in edges(g)
        e+=1
        h[edge.src,e] = true
        h[edge.dst,e] = true
    end
    h
end

"""
    Base.size(h::Hypergraph)

Returns the size of Hypergraph `h`.
The result is a tuple of the number of vertices and the number of hyperedges

"""
Base.size(h::Hypergraph) = (nhv(h), nhe(h))

"""
    Base.getindex(h::Hypergraph, idx::Vararg{Int,2})

Returns a value for a given vertex-hyperedge pair `idx` for a hypergraph `h`.
If a vertex does not belong to a hyperedge `nothing` is returned.

"""
@inline function Base.getindex(h::Hypergraph, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    get(h.v2he[idx[1]], idx[2], nothing)
end

"""
    Base.setindex!(h::Hypergraph, ::Nothing, idx::Vararg{Int,2})

Removes a vertex from a given hyperedge for a hypergraph `h` and a given vertex-hyperedge pair `idx`.
Note that trying to remove a vertex from a hyperedge when it is not present will not throw an error.

"""
@inline function Base.setindex!(h::Hypergraph, ::Nothing, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    pop!(h.v2he[idx[1]], idx[2], nothing)
    pop!(h.he2v[idx[2]], idx[1], nothing)
    h
end

"""
    Base.setindex!(h::Hypergraph, v::Real, idx::Vararg{Int,2})

Adds a vertex to a hyperedge (represented by indices `idx`) and assigns value
`v` to be stored with that assignment.

"""
@inline function Base.setindex!(h::Hypergraph, v::Real, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    h.v2he[idx[1]][idx[2]] = v
    h.he2v[idx[2]][idx[1]] = v
    h
end

"""
    getvertices(h::Hypergraph, he_id::Int)

Returns vertices from a hypergraph `a` for a given hyperedge `he_id`.

"""
@inline getvertices(h::Hypergraph, he_id::Int) = h.he2v[he_id]

"""
    gethyperedges(h::Hypergraph, v_id::Int)

Returns hyperedges for a given vertex `v_id` in a hypergraph `h`.

"""
@inline gethyperedges(h::Hypergraph, v_id::Int) = h.v2he[v_id]

"""
    add_vertex!(h::Hypergraph{T, V, E}; hyperedges::Dict{Int,T} = Dict{Int,T}(),
                v_meta::Union{V,Nothing} = nothing) where {T <: Real, V, E}

Adds a vertex to a given hypergraph `h`. Optionally, the vertex can be added
to existing hyperedges. The `hyperedges` parameter presents a dictionary
of hyperedge identifiers and values stored at the hyperedges.
Additionally, a value can be stored with the vertex using the `v_meta` keyword parameter.

"""
function add_vertex!(h::Hypergraph{T, V, E}; hyperedges::Dict{Int,T} = Dict{Int,T}(),
                     v_meta::Union{V, Nothing} = nothing) where {T <: Real, V, E}
    @boundscheck (checkbounds(h,1,k) for k in keys(hyperedges))
    push!(h.v2he,hyperedges)
    ix = length(h.v2he)
    for k in keys(hyperedges)
        h[ix,k]=hyperedges[k]
    end
    push!(h.v_meta, v_meta)
    ix
end

"""
    remove_vertex!(h::Hypergraph, v::Int)

Removes the vertex `v` from a given hypergraph `h`.
Note that running this function will cause reordering of vertices in the hypergraph:
the vertex `v` will replaced by the last vertex
of the hypergraph and the list of vertices will be shrunk.
"""
function remove_vertex!(h::Hypergraph, v::Int)
    n = nhv(h)
    if v < n
        h.v2he[v] = h.v2he[n]
        h.v_meta[v] = h.v_meta[n]
    end

    for hv in h.he2v
        if v < n && haskey(hv, n)
            hv[v] = hv[n]
            delete!(hv, n)
        else
            delete!(hv, v)
        end
    end
    resize!(h.v2he, length(h.v2he) - 1)
    h
end

"""
    add_hyperedge!(h::Hypergraph{T, V, E}; vertices::Dict{Int,T} = Dict{Int,T}(),
                   he_meta::Union{E,Nothing}=nothing) where {T <: Real, V, E}

Adds a hyperedge to a given hypergraph `h`. Optionally, existing vertices can be added
to the created hyperedge. The paramater `vertices` represents a dictionary
of vertex identifiers and values stored at the hyperedges.
Additionally, a value can be stored with the hyperedge using the `he_meta` keyword parameter.

"""
function add_hyperedge!(h::Hypergraph{T, V, E}; vertices::Dict{Int,T} = Dict{Int,T}(),
                        he_meta::Union{E,Nothing}=nothing) where {T <: Real, V, E}
    @boundscheck (checkbounds(h,k,1) for k in keys(vertices))
    push!(h.he2v,vertices)
    ix = length(h.he2v)
    for k in keys(vertices)
        h[k,ix]=vertices[k]
    end
    push!(h.he_meta, he_meta)
    ix
end

"""
    set_vertex_meta!(h::Hypergraph{T, V, E}, new_value::Union{V,Nothing}, id::Int)
        where {T <: Real, V, E}

Sets a new meta value `new_value` for the vertex `id` in the hypegraph `h`.

"""
function set_vertex_meta!(h::Hypergraph{T, V, E}, new_value::Union{V,Nothing}, id::Int) where {T <: Real, V, E}
    checkbounds(h.v_meta, id)
    h.v_meta[id] = new_value
    h.v_meta
end

"""
    get_vertex_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}

Returns a meta value stored at the vertex `id` in the hypergraph `h`.

"""
function get_vertex_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}
    checkbounds(h.v_meta, id)
    h.v_meta[id]
end

"""
    set_hyperedge_meta!(h::Hypergraph{T, V, E}, new_value::Union{E,Nothing}, id::Int)
        where {T <: Real, V, E}

Sets a new meta value `new_value` for the hyperedge `id` in the hypegraph `h`.

"""
function set_hyperedge_meta!(h::Hypergraph{T, V, E}, new_value::Union{E,Nothing}, id::Int) where {T <: Real, V, E}
    checkbounds(h.he_meta, id)
    h.he_meta[id] = new_value
    h.he_meta
end

"""
    get_hyperedge_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}

Returns a meta value stored at the hyperedge `id` in the hypergraph `h`.

"""
function get_hyperedge_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}
    checkbounds(h.he_meta, id)
    h.he_meta[id]
end

"""
    nhe(h::Hypergraph{T, V, E}) where {T <: Real, V, E}

Return the number of hyperedges in the hypergraph `h`.
"""
function nhe(h::Hypergraph{T, V, E}) where {T <: Real, V, E}
    length(h.he2v)
end

"""
    nhv(h::Hypergraph{T, V, E}) where {T <: Real, V, E}

Return the number of vertices in the hypergraph `h`.
"""
function nhv(h::Hypergraph{T, V, E}) where {T <: Real, V, E}
    length(h.v2he)
end

# TODO needs validate_hypergraph!(h::Hypergraph{T})


text:
using StructuresKit

#Test bimoment capacity equation against Bob Glauz's ballot example.
#see /testfiles/InternalForcesTest5/495A_Examples_v2.pdf
#check nominal and ASD calculations

Fy = 55
Cw = 11.9
Wn = 8.82

ASDorLRFD = 2   #nominal capacity
Bn = AISIS10024.h411(Cw, Fy, Wn, ASDorLRFD)

ASDorLRFD = 0   #ASD
BnASD = AISIS10024.h411(Cw, Fy, Wn, ASDorLRFD)

BnGlauz = 74.2
BnASDGlauz = 44.4

Error = abs((Bn - BnGlauz)/ BnGlauz)
ErrorASD = abs((BnASD - BnASDGlauz)/ BnASDGlauz)

#accept 1% error from numerical solution
@test Error <= 0.01
@test ErrorASD <= 0.01


text:
#= This file is a part of URIValidator. License is MIT:
https://github.com/rafpyprog/URIValidator.jl/blob/master/LICENSE.md =#

type URIComponents
    scheme
    authority
    path
    query
    fragment
    components

    function URIComponents(splited_URI)
        splited_URI_size = length(splited_URI)
        if splited_URI_size != 5
            message = "Splited URI length($splited_URI_size) is not valid. Should be 5"
            throw(ArgumentError(message))
        end

        this = new()
        this.scheme = splited_URI["scheme"]
        this.authority = splited_URI["authority"]
        this.path = splited_URI["path"]
        this.query = splited_URI["query"]
        this.fragment = splited_URI["fragment"]
        this.components = (this.scheme, this.authority, this.path, this.query,
                           this.fragment)
        return this
    end
end


text:
"""
	TreeStructure

This defines the tree structure.
A tree has a name, a list of parents, a list of states, list of probabilities, and a list of children.
Name - a string representing the name of the tree.
List of parents - shows the parent of each node in the tree.
List of states - contains the state of each node in the tree.
List of probabilities - shows the conditional probability of moving from one node to another.
List of children - contains the children of every parent node.
"""

using Distributions
using Statistics
using PyPlot

mutable struct Tree                                                                               #Tree structure
    name::String                                                                                #name of the tree
    parent::Vector{Int64}                                                                    # parents of nodes in the tree
    children::Vector{Vector{Int64}}                                             # successor nodes of each parent
    state::Matrix{Float64}                                                                 # states of nodes in the tree
    probability::Matrix{Float64}                                                     #probability to go from one node to another

"""
	Children(parent::Vector{Int64})

This function will create an array for the children of each of the parent nodes.
We need this array fo children for performance purposes in the stochastic approximation step.
So that instead of having it as a fucntion, we rather have it as an array so that we can only access the array
"""
    function Children(parent::Vector{Int64})
        allchildren = Vector{Vector{Int64}}([])
        for node in unique(parent)
            push!(allchildren,[i for i = 1:length(parent) if parent[i] == node])
        end
        return allchildren
    end
    Tree(name::String,parent::Vector{Int64},
    children::Vector{Vector{Int64}},
    state::Matrix{Float64},probability::Matrix{Float64})=new(name,parent,Children(parent),state,probability)

"""
	Tree(identifier::Int64)

This returns some examples of predefirned trees.
These are (0,302,303,304,305,306,307,402,404,405).
You can call any of the above tree and plot to see the properties of the tree.
"""
    function Tree(identifier::Int64)
        self= new()
        if identifier==0
            self.name = "Empty Tree"
            self.parent = zeros(Int64,0)
            self.children = Children(self.parent)
            self.state = zeros(Float64,0)'
            self.probability = ones(Float64,0)'
        elseif identifier == 302
            self.name = "Tree 1x2x2"
            self.parent =[0,1,1,2,2,3,3]
            self.children = Children(self.parent)
            self.state = [2.0 2.1 1.9 4.0 1.0 3.0]'
            self.probability = [1.0 0.5 0.5 0.5 0.5 0.5 0.5]'
        elseif identifier == 303
            self.name = "Tree 1x1x4"
            self.parent = [0,1,2,2,2,2]
            self.children = Children(self.parent)
            self.state = [3.0 3.0 6.0 4.0 2.0 0.0]'
            self.probability = [1.0 1.0 0.25 0.25 0.25 0.25]'
        elseif identifier == 304
            self.name = "Tree 1x4x1x1"
            self.parent = [0,1,2,0,4,5,0,7,8,0,10,11]
            self.children = Children(self.parent)
            self.state = [0.1 2.1 3.0 0.1 1.9 1.0 0.0 -2.9 -1.0 -0.1 -3.1 -4.0]'
            self.probability = [0.14 1.0 1.0 0.06 1.0 1.0 0.48 1.0 1.0 0.32 1.0 1.0]'
        elseif identifier == 305
            self.name = "Tree 1x1x4"
            self.parent = [0,1,2,2,2,2]
            self.children = Children(self.parent)
            self.state = [0.0 10.0 28.0 22.0 21.0 20.0]'
            self.probability = [1.0 1.0 0.25 0.25 0.25 0.25]'
        elseif identifier == 306
            self.name = "Tree 1x2x2"
            self.parent = [0,1,1,2,2,3,3]
            self.children = Children(self.parent)
            self.state = [0.0 10.0 10.0 28.0 22.0 21.0 20.0]'
            self.probability = [1.0 0.5 0.5 0.5 0.5 0.5 0.5]'
        elseif identifier == 307
            self.name = "Tree 1x4x1"
            self.parent = [0,1,1,1,1,2,3,4,5]
            self.children = Children(self.parent)
            self.state = [0.0 10.0 10.0 10.0 10.0 28.0 22.0 21.0 20.0]'
            self.probability = [1.0 0.25 0.25 0.25 0.25 1.0 1.0 1.0 1.0]'
        elseif identifier == 401
            self.name = "Tree 1x1x2x2"
            self.parent = [0,1,2,2,3,3,4,4]
            self.children = Children(self.parent)
            self.state = [10.0 10.0 8.0 12.0 9.0 6.0 10.0 13.0]'
            self.probability = [1.0 1.0 0.66 0.34 0.24 0.76 0.46 0.54]'
        elseif identifier == 402
            self.name = "Tree 1x2x2x2"
            self.parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            self.children = Children(self.parent)
            #self.state= [10.0 11.0 8.0 12.0 9.0 6.0 10.0 13.0 10.01 8.01 12.01 9.01 4.0 11.0 12.99]'
            self.state = [10.0 12.0 8.0 15.0 11.0 9.0 5.0 18.0 16.0 13.0 11.0 10.0 7.0 6.0 3.0]'
            self.probability = [1.0 0.8 0.7 0.3 0.2 0.8 0.4 0.6 0.2 0.5 0.5 0.4 0.6 0.7 0.3]'
        elseif identifier == 4022
            self.name = "Tree 1x2x2x2"
            self.parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            self.children = Children(self.parent)
            #self.state= [10.0 11.0 8.0 12.0 9.0 6.0 10.0 13.0 10.01 8.01 12.01 9.01 4.0 11.0 12.99]'
            self.state = [10.0 0.0; 12.0 1.0; 8.0 -1.0; 15.0 2.0; 11.0 1.0; 9.0 -0.5; 5.0 -2.0; 18.0 3.0; 16.0 1.8; 13.0 0.9; 11.0 0.2; 10.0 0.0; 7.0 -1.2; 6.0 -2.0; 3.0 -3.2]
            self.probability = [1.0 0.8 0.7 0.3 0.2 0.8 0.4 0.6 0.2 0.5 0.5 0.4 0.6 0.7 0.3]'
        elseif identifier ==404
            self.name = "Tree 1x2x2x2"
            self.parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            self.children = Children(self.parent)
            self.state = (0.2+0.6744).*[0.0 1.0 -1.0 2.0 0.1 0.0 -2.0 3.0 1.1 0.9 -1.1 1.2 -1.2 -0.8 -3.2]'
            self.probability = [1.0 0.3 0.7 0.2 0.8 0.1 0.9 0.5 0.5 0.6 0.4 0.4 0.6 0.3 0.7]'
        elseif identifier == 405
            self.name = "Tree 5"
            self.parent = [0,1,2,3,3,2,6,6,6,1,10,10,11,12,12,12]
            self.children = Children(self.parent)
            self.state = [7.0 12.0 14.0 15.0 13.0 9.0 11.0 10.0 8.0 4.0 6.0 2.0 5.0 0.0 1.0 3.0]'
            self.probability = [1.0 0.7 0.4 0.7 0.3 0.6 0.2 0.3 0.5 0.3 0.2 0.8 1.0 0.3 0.5 0.2]'
        end
        return self
    end
"""
	Tree(spec::Vector{Int64}, dimension=Int64[])

Instead of having a predefined tree, we can just provide the topology of the tree and then using it, we produce the tree.
Given the bushiness of the tree, the following function produces the tree.
For example, bushiness = 1x2x3x3 means that we have 18 leaves at the end and at stage 1, we have 1 node,
stage 2 has 2x1 nodes, stage 3 has 3x2x1 nodes and stage 4 has 3x3x2x1 nodes.
We generate normal random  variables for the states of the nodes of the tree.
"""

    function Tree(spec::Vector{Int64}, dimension=Int64[])
        self = new()
        if isempty(dimension)
            dimension = 1
        end
        leaves = spec
        for stage = 1 : length(spec)
            if stage == 1
                leaves = 1
                self.name = "Tree $(spec[1])"
                self.parent = zeros(Int64,spec[1])
                self.children = Children(self.parent)
                self.state = randn(spec[1],dimension)
                self.probability = ones(spec[1],1)
            else
                leaves = leaves * spec[stage-1]
                newleaves = vec(ones(Int64,spec[stage]) .* transpose(length(self.parent) .+ (1 .- leaves:0)))
                self.parent =  vcat(self.parent, newleaves)
                self.children = Children(self.parent)
                self.state = vcat(self.state,randn(length(newleaves),dimension))
                self.name = "$(self.name)x$(spec[stage])"
                tmp = rand(Uniform(0.3,1.0), spec[stage], leaves)
                tmp = tmp ./ (transpose(ones(1,spec[stage])) .* sum(tmp,dims=1))
                self.probability = vcat(self.probability,vec(tmp))
            end
        end
        return self
    end
end

"""
	stage(trr::Tree, node=Int64[])

The stage function returns the stage of each node in the tree
"""

function stage(trr::Tree,node = Int64[])
    if isempty(node)
        node = 1 : length(trr.parent)
    elseif isa(node, Int64)
        node = Int64[node]
    end
    #stage = zeros(length(collect(node))) #you can also use this
    stage = zero(node)
    for i = 1:length(node)
        pred = node[i]
        while pred > 0 && trr.parent[pred] > 0
            pred = trr.parent[pred]
            stage[i] += 1
        end
    end
    return stage
end

"""
	height(trr::Tree)

This function returns the height of the tree which is just the maximum number of the stages of each node
"""

function height(trr::Tree)
    return maximum(stage(trr))
end

"""
	leaves(trr::Tree,node=Int64[])

This function produces the terminal nodes of the tree.
The terminal nodes or the leaves are all those nodes which doesn't have children nodes.
They are the nodes which are not parents.
The function also returns the indexes (omegas) of these nodes as well as the conditional probabilities (prob) of reaching of the leaves from the root node
"""

function leaves(trr::Tree, node = Int64[])
    nodes = 1:length(trr.parent)
    leaves = setdiff(nodes,trr.parent)
    #leaves are all nodes which are not parents
    #setdiff(A,B) finds all the elements that belongs to A and are not in B
    omegas = 1:length(leaves)
    if !isempty(node) && isa(node,Int64)
        node = Int64[node]
    end
    if !isempty(node) && (0 ∉ node)
        omegas = Set{Int64}()
        nodes = leaves
        while any(j !=0 for j ∈ nodes)
            omegas = union(omegas,(ind for (ind,j) ∈ enumerate(nodes) if j ∈ node))
            nodes = Int64[trr.parent[max(0,j)] for j ∈ nodes]
        end
        omegas = collect(omegas)
    end
    leaves = Int64[leaves[j] for j ∈ omegas]
    prob = ones(Float64,length(leaves))
    nodes = leaves
    while any(j !=0 for j ∈ nodes)
        prob = prob .* trr.probability[[j for j in nodes]]
        nodes = Int64[trr.parent[max(0,j)] for j ∈ nodes]
    end
    return leaves,omegas, prob
end

"""
	nodes(trr::Tree,t=Int64[])

This function gives the nodes in the tree, generally the range of the nodes in the tree
"""

function nodes(trr::Tree, t = Int64[])
    nodes = 1 : length(trr.parent)
    if isempty(t) #if stage t is not given, return all nodes of the tree
        return nodes
    else # else return nodes at the given stage t
        stg = stage(trr)
        return Int64[i for i in nodes if stg[i] == t]
    end
end

"""
	root(trr::Tree,nodes=Int64[])
This function returns the starting node of the tree. It also provides the path from the starting node of the tree to the node indicated.
If the node is not indicated, it will return the starting node (which we call the root node).
Else, it will give a sequence of nodes from the root to the indicated node.
"""

function root(trr::Tree,nodes = Int64[])
    if isempty(nodes)
        nodes = trr.children[1]
    elseif isa(nodes,Int64)
        nodes = Int64[nodes]
    end
    root = 1 : length(trr.parent)
    for i in nodes
        iroot = Vector{Int64}([])
        tmp = i
        while tmp > 0
            push!(iroot,tmp)
            tmp = trr.parent[tmp]
        end
        root = Int64[i for i in reverse(iroot) if i in root]
    end
    return root
end

"""
	partTree(trr::Tree)

Sometimes, we are dealing with 2D-state trees and at the end of the stochastic approximation step, we get an aproximated tree with
states in higher dimension. So we can partition the tree into different trees with the same parents, children and probabilities but
different states. The following function helps us to do so.
"""

function partTree(trr::Tree)
    trees = Tree[]
    for col = 1:size(trr.state,2)
        subtree = Tree("Tree of state $col",trr.parent,trr.children,hcat(trr.state[:,col]),trr.probability)
        push!(trees,subtree)
    end
    return trees
end

"""
	buildProb!(trr::Tree,probabilities::Array{Float64,2})

In the stochastic approximation step, we only iterate over the probabilities of the leaves. The length of the array of probabilities we get at the
end is smaller than the length of the parents. So we have to build the probabilites for the remaining nodes.
"""
function buildProb!(trr::Tree,probabilities::Array{Float64,2})
    leaf,omegas,probaLeaf = leaves(trr)
    if length(probabilities) == length(nodes(trr))
        trr.probability .= probabilities
    else
        probabilities .= max.(0.0,probabilities)
        j = leaf;i = trr.parent[j];j .= j[i.> 0];i .= i[i .> 0]
        trr.probability .= zeros(length(trr.state[:,1]),1)
        trr.probability[j] = probabilities
        while !isempty(i)
            for k = 1:length(i)
                trr.probability[i[k]] = trr.probability[i[k]] + trr.probability[j[k]]
            end
            trr.probability[j] .= trr.probability[j] ./ trr.probability[i]
            trr.probability[isnan.(trr.probability)] .= 0.0
            j = unique(i); i = trr.parent[j]; j = j[i .> 0] ; i = i[i .> 0]
        end
    end
    return trr.probability
end

"""
	treeplot(trr::Tree,fig=1)

This function plots a tree in 1D. It also annotates the density of probabilities of the nodes in the tree.
"""

function treeplot(trr::Tree, fig= 1)
    if !isempty(fig)
        figure(fig)
    end
    #suptitle(trr.name, fontsize = 14)
    trs = subplot2grid((1,4), (0,0), colspan=3)
    title("states")
    stg = stage(trr)
    xticks(range(1, stop = height(trr) + 1))                                              #Set the ticks on the x-axis
    xlabel("stage, time")
    trs.spines["top"].set_visible(false)                                                       #remove the line of the box at the top
    trs.spines["right"].set_visible(false)		                                       # remove the line of the box at the right
    for i in range(1,stop = length(trr.parent))
        if stg[i] > 0
            trs.plot([stg[i]-1,stg[i]],[trr.state[trr.parent[i]],trr.state[i]])
        end
    end
    prs = subplot2grid((1,4), (0,3))
    title("probabilities")
    prs.spines["top"].set_visible(false)
    prs.spines["left"].set_visible(false)
    prs.spines["right"].set_visible(false)

    (Yi,_,probYi) = leaves(trr)
    Yi = [trr.state[i] for i in Yi]
    nY = length(Yi)
    h = 1.05*std(Yi)/ (nY^0.2) + 1e-3                                                           #Silverman rule of thumb
    trs.set_ylim(minimum(Yi)-h, maximum(Yi)+h)
    prs.set_ylim(minimum(Yi)-h, maximum(Yi)+h)
    yticks(())                                                                                #remove the ticks on probability plot
    t = LinRange(minimum(Yi)-h, maximum(Yi)+h, 100)                           #100 points on probability plot
    #density = zeros(length(collect(t)))
    density = zero(t)
    for (i, ti) in enumerate(t)
        for (j, xj) in enumerate(Yi)
            tmp = (xj - ti) / h
            density[i] += probYi[j]* 35/32 * max(1.0 -tmp^2, 0.0)^3 /h #triweight kernel
        end
    end
    plot(density, t)
    prs.fill_betweenx(t, 0 , density)
end

"""
	plotD(newtree::Tree)

This functions is used to plot a tree when we are dealing with a tree with states in 2D.
It returns a plot with two trees with differet states but the same parents, children and probabilities.
"""

function plotD(newtree::Tree)
    fig = figure(figsize = (10,6))
    stg = stage(newtree)
    for rw = 1:size(newtree.state,2)
      ax = subplot(1,size(newtree.state,2),rw)
      ax.spines["top"].set_visible(false)                                                  #remove the box top
      ax.spines["right"].set_visible(false)                                               #remove the box right
      for i in range(1,stop = length(newtree.parent))
          if stg[i] > 0
              ax.plot([stg[i]-1,stg[i]],[newtree.state[:,rw][newtree.parent[i]],newtree.state[:,rw][i]])
          end
          xlabel("stage, time")
          ylabel("states")
      end
      xticks(unique(stg))
    end
end


text:
module NFFT

using Base.Cartesian
using FFTW
using Distributed
using SparseArrays
using LinearAlgebra

export NFFTPlan, nfft, nfft_adjoint, ndft, ndft_adjoint

include("windowFunctions.jl")
include("precomputation.jl")

#=
Some internal documentation (especially for people familiar with the nfft)

- The window is precomputed during construction of the NFFT plan
  When performing the nfft convolution, the LUT of the window is used to
  perform linear interpolation. This approach is reasonable fast and does not
  require too much memory. There are, however alternatives known that are either
  faster or require no extra memory at all.

The non-exported functions apodization and convolve are implemented
using Cartesian macros, that may not be very readable.
This is a conscious decision where performance has outweighed readability.
More readable versions can be (and have been) written using the CartesianRange approach,
but at the time of writing this approach require *a lot* of memory.
=#

#=
D is the number of dimensions of the array to be transformed.
DIM is the dimension along which the array is transformed.
DIM == 0 is the ordinary NFFT, i.e., where all dimensions are transformed.
DIM is a type parameter since it allows the @generated macro to
compile more efficient methods.
=#
mutable struct NFFTPlan{D,DIM,T}
    N::NTuple{D,Int64}
    M::Int64
    x::Matrix{T}
    m::Int64
    sigma::T
    n::NTuple{D,Int64}
    K::Int64
    windowLUT::Vector{Vector{T}}
    windowHatInvLUT::Vector{Vector{T}}
    forwardFFT::FFTW.cFFTWPlan{Complex{T},-1,true,D}
    backwardFFT::FFTW.cFFTWPlan{Complex{T},1,true,D}
    tmpVec::Array{Complex{T},D}
    B::SparseMatrixCSC{T,Int64}
end

@inline dim(::NFFTPlan{D,DIM}) where {D, DIM} = DIM

"""
        NFFTPlan(x, N, ...) -> plan

Compute `D` dimensional NFFT plan for sampling locations `x` (a vector or a `D`-by-`M` matrix) that can be applied on arrays of size `N` (a tuple of length `D`).

The optional arguments control the accuracy.

It takes as optional keywords all the keywords supported by `plan_fft` function (like
`flags` and `timelimit`).  See documentation of `plan_fft` for reference.
"""
NFFTPlan

@enum PrecomputeFlags begin
  LUT = 1
  FULL = 2
end

function NFFTPlan(x::Matrix{T}, N::NTuple{D,Int}, m=4, sigma=2.0,
                       window=:kaiser_bessel, K=2000;
                       precompute::PrecomputeFlags=LUT, kwargs...) where {D,T}
    if D != size(x,1)
        throw(ArgumentError())
    end

    n = ntuple(d->round(Int,sigma*N[d]), D)

    tmpVec = zeros(Complex{T}, n)

    M = size(x,2)

    FP = plan_fft!(tmpVec; kwargs...)
    BP = plan_bfft!(tmpVec; kwargs...)

    # Create lookup table
    win, win_hat = getWindow(window)

    windowLUT = Vector{Vector{T}}(undef,D)
    windowHatInvLUT = Vector{Vector{T}}(undef,D)
    for d=1:D
        windowHatInvLUT[d] = zeros(T, N[d])
        for k=1:N[d]
            windowHatInvLUT[d][k] = 1. / win_hat(k-1-N[d]/2, n[d], m, sigma)
        end
    end

    if precompute == LUT
      Z = round(Int,3*K/2)
      for d=1:D
          windowLUT[d] = zeros(T, Z)
          for l=1:Z
              y = ((l-1) / (K-1)) * m/n[d]
              windowLUT[d][l] = win(y, n[d], m, sigma)
          end
      end

      B = sparse([],[],Float64[])
    else
      U1 = ntuple(d-> (d==1) ? 1 : (2*m+1)^(d-1), D)
      U2 = ntuple(d-> (d==1) ? 1 : prod(n[1:(d-1)]), D)
      B = precomputeB(win, x, n, m, M, sigma, T, U1, U2)
    end

    NFFTPlan{D,0,T}(N, M, x, m, sigma, n, K, windowLUT, windowHatInvLUT,
                    FP, BP, tmpVec, B )
end

NFFTPlan(x::AbstractMatrix{T}, N::NTuple{D,Int}, rest...; kwargs...) where {D,T} =
    NFFTPlan(collect(x), N, rest...; kwargs...)

NFFTPlan(x::AbstractVector, N::Integer, rest...; kwargs...) =
    NFFTPlan(reshape(x,1,length(x)), (N,), rest...; kwargs...)

# Directional NFFT
"""
        NFFTPlan(x, d, N, ...) -> plan

Compute *directional* NFFT plan:
A 1D plan that is applied along dimension `d` of a `D` dimensional array of size `N` with sampling locations `x` (a vector).

It takes as optional keywords all the keywords supported by `plan_fft` function (like
`flags` and `timelimit`).  See documentation of `plan_fft` for reference.
"""
function NFFTPlan(x::AbstractVector{T}, dim::Integer, N::NTuple{D,Int64}, m=4,
                       sigma=2.0, window=:kaiser_bessel, K=2000; kwargs...) where {D,T}
    n = ntuple(d->round(Int, sigma*N[d]), D)

    sz = [N...]
    sz[dim] = n[dim]
    tmpVec = Array{Complex{T}}(undef,sz...)

    M = length(x)

    FP = plan_fft!(tmpVec, dim; kwargs...)
    BP = plan_bfft!(tmpVec, dim; kwargs...)

    # Create lookup table
    win, win_hat = getWindow(window)

    windowLUT = Vector{Vector{T}}(undef,1)
    Z = round(Int, 3*K/2)
    windowLUT[1] = zeros(T, Z)
    for l = 1:Z
        y = ((l-1) / (K-1)) * m/n[dim]
        windowLUT[1][l] = win(y, n[dim], m, sigma)
    end

    windowHatInvLUT = Vector{Vector{T}}(undef,1)
    windowHatInvLUT[1] = zeros(T, N[dim])
    for k = 1:N[dim]
        windowHatInvLUT[1][k] = 1. / win_hat(k-1-N[dim]/2, n[dim], m, sigma)
    end

    B = sparse([],[],Float64[])

    NFFTPlan{D,dim,T}(N, M, reshape(x,1,M), m, sigma, n, K, windowLUT,
                      windowHatInvLUT, FP, BP, tmpVec, B)
end

function NFFTPlan(x::Matrix{T}, dim::Integer, N::NTuple{D,Int}, m=4, sigma=2.0,
                       window=:kaiser_bessel, K=2000; kwargs...) where {D,T}
    if size(x,1) != 1 && size(x,2) != 1
        throw(DimensionMismatch())
    end

    NFFTPlan(vec(x), dim, N, m, sigma, window, K; kwargs...)
end

function Base.show(io::IO, p::NFFTPlan{D,0}) where D
    print(io, "NFFTPlan with ", p.M, " sampling points for ", p.N, " array")
end

function Base.show(io::IO, p::NFFTPlan{D,DIM}) where {D,DIM}
    print(io, "NFFTPlan with ", p.M, " sampling points for ", p.N, " array along dimension ", DIM)
end

@generated function consistencyCheck(p::NFFTPlan{D,DIM,T}, f::AbstractArray{U,D},
                                     fHat::AbstractArray{Y}) where {D,DIM,T,U,Y}
    quote
        if $DIM == 0
            fHat_test = (p.M == length(fHat))
        elseif $DIM > 0
            fHat_test = @nall $D d -> ( d == $DIM ? size(fHat,d) == p.M : size(fHat,d) == p.N[d] )
        end

        if p.N != size(f) || !fHat_test
            throw(DimensionMismatch("Data is not consistent with NFFTPlan"))
        end
    end
end

### nfft functions ###

"""
        nfft!(p, f, fHat) -> fHat

Calculate the NFFT of `f` with plan `p` and store the result in `fHat`.

Both `f` and `fHat` must be complex arrays.
"""
function nfft!(p::NFFTPlan{D,DIM,T}, f::AbstractArray, fHat::StridedArray, verbose=false) where {D,DIM,T}
    consistencyCheck(p, f, fHat)

    fill!(p.tmpVec, zero(Complex{T}))
    t1 = @elapsed @inbounds apodization!(p, f, p.tmpVec)
    if nprocs() == 1
        t2 = @elapsed p.forwardFFT * p.tmpVec # fft!(p.tmpVec) or fft!(p.tmpVec, dim)
    else
        t2 = @elapsed dim(p) == 0 ? fft!(p.tmpVec) : fft!(p.tmpVec, dim(p))
    end
    t3 = @elapsed @inbounds convolve!(p, p.tmpVec, fHat)
    if verbose
      @info "Timing: apod=$t1 fft=$t2 conv=$t3"
    end
    return fHat
end

"""
        nfft(p, f) -> fHat

For a **non**-directional `D` dimensional plan `p` this calculates the NFFT of a `D` dimensional array `f` of size `N`.
`fHat` is a vector of length `M`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function nfft(p::NFFTPlan{D,0,T}, f::AbstractArray{U,D}, args...) where {D,T,U}
    fHat = zeros(Complex{T}, p.M)
    nfft!(p, f, fHat, args...)
    return fHat
end

function nfft(x, f::AbstractArray{T,D}, rest...; kwargs...) where {D,T}
    p = NFFTPlan(x, size(f), rest...; kwargs...)
    return nfft(p, f)
end

function nfft(p::NFFTPlan{D,DIM,T}, f::AbstractArray{U,D}, args...) where {D,DIM,T,U}
    sz = [p.N...]
    sz[DIM] = p.M
    fHat = Array{Complex{T}}(undef,sz...)
    nfft!(p, f, fHat, args...)
    return fHat
end

"""
        nfft_adjoint!(p, fHat, f) -> f

Calculate the adjoint NFFT of `fHat` and store the result in `f`.

Both `f` and `fHat` must be complex arrays.
"""
function nfft_adjoint!(p::NFFTPlan, fHat::AbstractArray, f::StridedArray, verbose=false)
    consistencyCheck(p, f, fHat)

    t1 = @elapsed @inbounds convolve_adjoint!(p, fHat, p.tmpVec)
    if nprocs() == 1
        t2 = @elapsed p.backwardFFT * p.tmpVec # bfft!(p.tmpVec) or bfft!(p.tmpVec, dim)
    else
        t2 = @elapsed dim(p) == 0 ? bfft!(p.tmpVec) : bfft!(p.tmpVec, dim(p))
    end
    t3 = @elapsed @inbounds apodization_adjoint!(p, p.tmpVec, f)
    if verbose
      @info "Timing: conv=$t1 fft=$t2 apod=$t3"
    end
    return f
end

"""
        nfft_adjoint(p, f) -> fHat

For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFFT of a length `M` vector `fHat`
`f` is a `D` dimensional array of size `N`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function nfft_adjoint(p::NFFTPlan{D,DIM,T}, fHat::AbstractArray{U}, args...) where {D,DIM,T,U}
    f = Array{Complex{T}}(undef,p.N)
    nfft_adjoint!(p, fHat, f, args...)
    return f
end

function nfft_adjoint(x, N, fHat::AbstractVector{T}, rest...; kwargs...) where T
    p = NFFTPlan(x, N, rest...; kwargs...)
    return nfft_adjoint(p, fHat)
end

include("directional.jl")
include("multidimensional.jl")
include("samplingDensity.jl")
include("NDFT.jl")

end


text:
module Voronoi

using MiniQhull
using LinearAlgebra

import ChainRulesCore: rrule, NoTangent

∧(x,y) = x[1,:].*y[2,:] .- x[2,:].*y[1,:]

"""
    boundary(d)

Returns the boundary vertices of a unit cube in `d` dimensions.
"""
function boundary(d)
    d == 2 && return Float32.([[-1;-1] [-1;+1] [+1;+1] [+1;-1]])

    b = zeros(Float32, d, 2^d)
    for i in 0:(2^d-1)
        for n in 0:(d-1)
            b[n+1,i+1] = (((i >> n) & 1) == 1) ? +1 : -1
        end
    end
    return b
end

const NB = 4

"""
    tessellation(q)

Construct a voronoi tessellation from generating points `q`.
Return _just_ vertices of construction.
"""
function tessellation(q)
    # q = hcat(boundary, x)
    triangulation = delaunay(q)

    # order triangulation
    areas = (q[:,triangulation[2,:]] .- q[:,triangulation[1,:]]) ∧ (q[:,triangulation[3,:]] .- q[:,triangulation[2,:]])
    orientation = sign.(areas)
    for i in findall(orientation .== -1)
        triangulation[2,i], triangulation[3,i] = triangulation[3,i], triangulation[2,i]
    end

    # compute voronoi vertices
    r = zeros(eltype(q), size(q,1), size(triangulation,2))
    for t in 1:size(triangulation,2)
        for i in 1:3
            j = (i % 3) + 1
            k = (j % 3) + 1

            α,β,γ = triangulation[i,t], triangulation[j,t], triangulation[k,t]

            qα² = sum((q[:,α]).^2)
            r[1,t] += qα²*(q[2,γ] - q[2,β])
            r[2,t] -= qα²*(q[1,γ] - q[1,β])
        end
    end

    r[1,:] = r[1,:] ./ abs.(areas)
    r[2,:] = r[2,:] ./ abs.(areas)

    triangulation, r
end

"""
    areas(x)

Compute the areas of all 2-dimensional simplices generated by the delaunay construction for pointcloud `x`.
Assumes `x` is sized ``N \\times 2``
"""
function areas(x)
    q = hcat(boundary(2), x)
    triangulation = delaunay(q)

    a = 0.5*[
        let
            q[1,t[1]]*(q[2,t[2]]-q[2,t[3]]) +
            q[1,t[2]]*(q[2,t[3]]-q[2,t[1]]) +
            q[1,t[3]]*(q[2,t[1]]-q[2,t[2]])
        end for t in eachcol(triangulation)
    ]
    s = sign.(a)

    return s.*a
end

function rrule(::typeof(areas), x)
    q = hcat(boundary(2), x)
    triangulation = delaunay(q)

    a = 0.5*[
        let
            q[1,t[1]]*(q[2,t[2]]-q[2,t[3]]) +
            q[1,t[2]]*(q[2,t[3]]-q[2,t[1]]) +
            q[1,t[3]]*(q[2,t[1]]-q[2,t[2]])
        end for t in eachcol(triangulation)
    ]
    s = sign.(a)

    return s.*a, (∂a) -> let
        ∂x = zeros(size(x))

        for (i,t) in enumerate(eachcol(triangulation))
            if t[1] > NB
                ∂x[1,t[1]-NB] += (q[2,t[2]]-q[2,t[3]])*∂a[i]*s[i]
                ∂x[2,t[1]-NB] -= (q[1,t[2]]-q[1,t[3]])*∂a[i]*s[i]
            end
            if t[2] > 4
                ∂x[1,t[2]-NB] += (q[2,t[3]]-q[2,t[1]])*∂a[i]*s[i]
                ∂x[2,t[2]-NB] -= (q[1,t[3]]-q[1,t[1]])*∂a[i]*s[i]
            end
            if t[3] > 4
                ∂x[1,t[3]-NB] += (q[2,t[1]]-q[2,t[2]])*∂a[i]*s[i]
                ∂x[2,t[3]-NB] -= (q[1,t[1]]-q[1,t[2]])*∂a[i]*s[i]
            end
        end

        (NoTangent(), 0.5*∂x)
    end

end

# derivative of det given by adjugate (which is just rescaled inverse)
# see https://en.wikipedia.org/wiki/Adjugate_matrix

affine(simplex) = vcat(simplex, ones(Float32, 1, size(simplex,2)))
volume(simplex) = det(affine(simplex))

"""
    volumes(x)

Compute the volume of all d-dimensional simplices generated by the delaunay construction for pointcloud `x`.
Assumes `x` is sized ``N \\times d``
"""
function volumes(x)
    d = size(x,1)
    q = hcat(boundary(d), x)

    simplices = delaunay(q)

    Ω = [ volume(hcat((q[:,i] for i in simplex)...)) for simplex in eachcol(simplices) ]
    s = sign.(Ω)

    return s.*Ω
end

function rrule(::typeof(volumes), x)
    d = size(x,1)
    b = boundary(d)
    q = hcat(b, x)

    v₀ = size(b,2)
    simplices = delaunay(q)

    Ω = [ volume(hcat((q[:,i] for i in simplex)...)) for simplex in eachcol(simplices) ]
    s = sign.(Ω)

    Z = factorial(d)
    return (s.*Ω) ./ Z, function(∂Ω)
        ∂x = zeros(Float32,size(x))

        for (n,simplex) in enumerate(eachcol(simplices))
            ω = s[n]*Ω[n]*∂Ω[n]
            try
                Q = inv(affine(hcat((q[:,i] for i in simplex)...)))

                for (i,v) in enumerate(simplex)
                    if v > v₀
                        for j in 1:d
                            ∂x[j,v-v₀] += ω*Q[i,j]
                        end
                    end
                end
            catch error
                # XXX: do we need to do something more intelligent here?
                if isa(error, LinearAlgebra.SingularException)
                    continue
                else
                    throw(error)
                end
            end
        end

        return (NoTangent(), ∂x ./ Z)
    end
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
type Person
    name:: String
    height::Float64
end

names(Person) #> 2-element Array{Symbol,1}: :name  :height
Person.types  #> (String,Float64)

code_lowered(+, (Int, Int))
#1-element Array{Any,1}:
# :($(Expr(:lambda, {:x,:y}, {{},{{:x,:Any,0},{:y,:Any,0}},{}}, :(begin  # int.jl
#, line 33:
#        return box(Int64,add_int(unbox(Int64,x),unbox(Int64,y)))
#    end))))

code_typed(+, (Int, Int))
# 1-element Array{Any,1}:
#  :($(Expr(:lambda, {:x,:y}, {{},{{:x,Int64,0},{:y,Int64,0}},{}}, :(begin  # int.
# jl, line 33:
#         return box(Int64,add_int(x::Int64,y::Int64))::Int64
#     end::Int64))))

max_stars_repo_path: Module 1/Chapter07/reflection.jl

max_issues_repo_path: Module 1/Chapter07/reflection.jl

max_forks_repo_path: Module 1/Chapter07/reflection.jl


text:
import Dolo
import YAML
using AxisArrays

path = Dolo.pkg_path
@testset "testing model algos" begin

    @testset "testing mc models" begin
        fn = joinpath(path, "examples", "models", "rbc_mc.yaml")
        # fn = joinpath(path, "examples", "models", "LAMP.yaml")
        model_mc = Dolo.yaml_import(fn)

        drc = Dolo.ConstantDecisionRule(model_mc.calibration[:controls])
        # @time tid_res = Dolo.time_iteration_direct(model_mc, drc; maxit=20, verbose=true)
        # @time ti_res = Dolo.time_iteration(model_mc, tid_res.dr; maxit=20, verbose=false)
        @time ti_res = Dolo.time_iteration(model_mc; maxit=20, verbose=false)
        @time iti_res = Dolo.improved_time_iteration(model_mc; verbose=false, maxit=10000)
        @test Dolo.converged(iti_res)
        # @time drv = Dolo.evaluate_policy(model_mc, tid_res.dr; maxit=20, verbose=false)

        sim = Dolo.simulate(model_mc, ti_res.dr, model_mc.exogenous) #; N=100, T=20)

        sim = Dolo.simulate(model_mc, ti_res.dr, i0= 2)
        sim = Dolo.simulate(model_mc, ti_res.dr)
        #
        k = sim[Axis{:V}(:k)]
        n = sim[Axis{:V}(:n)]
        z = sim[Axis{:V}(:z)]
        y = sim[Axis{:V}(:y)]
        alpha= model_mc.calibration.flat[:alpha]
        exp.(z).*k.^alpha.*n.^(1-alpha)== y

        @test true
    end

    # # compare with prerecorded values
    # kvec = linspace(dr.grid.min[1], dr.grid.max[1], 10)
    # nvec = [dr(1, [k])[1] for k in kvec]
    # ivec = [dr(1, [k])[2] for k in kvec]
    # # compare  time_iteration_direct
    # nvec_d = [drd(1, [k])[1] for k in kvec]
    # ivec_d = [drd(1, [k])[2] for k in kvec]
    # @assert maximum(abs, nvec_d-nvec)<1e-4
    #
    # # compare  vfi
    # nvec_0 = [dr0(1, [k])[1] for k in kvec]
    # ivec_0 = [dr0(1, [k])[2] for k in kvec]
    # @assert maximum(abs, nvec_0-nvec)<1e-4

    # let's redo when model is stable !
    # ivec_test = [0.295977,  0.257538,  0.21566,  0.173564,  0.132103,  0.0915598,  0.0520067,  0.0134661,  7.01983e-6, 3.40994e-17]
    # nvec_test = [ 0.391997,  0.348033,  0.318369,  0.296276,  0.278821,  0.264487,  0.25239 ,  0.241974,  0.236604,  0.233779 ]
    # @assert maximum(abs(ivec-ivec_test))<1e-5
    # @assert maximum(abs(nvec-nvec_test))<1e-5

    @testset "testing iid models" begin

        fn = joinpath(path, "examples", "models", "rbc_iid.yaml")
        model = Dolo.yaml_import(fn)

        @time dr = Dolo.perturb(model)

        drc = Dolo.ConstantDecisionRule(model.calibration[:controls])

        # @time tid_res = Dolo.time_iteration_direct(model, drc; maxit=20, verbose=true)
        # @time ti_res = Dolo.time_iteration(model, tid_res.dr; maxit=20, verbose=false)
        @time ti_res = Dolo.time_iteration(model; maxit=20, verbose=false)
        @time iti_res = Dolo.improved_time_iteration(model; maxit=20, verbose=false)
        @test Dolo.converged(iti_res)
        #
        # @time sol_v = Dolo.value_iteration(model, tid_res.dr; maxit=20, verbose=true)

        # Dolo.simulate(model, tid_res.dr)

        s0 = model.calibration[:states] .+ 0.1
        sim = Dolo.simulate(model, ti_res.dr; s0=s0)
        Dolo.simulate(model, ti_res.dr; N=10)

        res = Dolo.response(model.exogenous, [0.01])
        irf = Dolo.response(model, ti_res.dr, :e_z)
        irf = Dolo.response(model, ti_res.dr, :e_z, -0.01)

        # check with random and Smolyak grid
        # @time Dolo.time_iteration(model, maxit=20, grid=Dict(:tag => :Random, :N => 200))
        # @time Dolo.time_iteration(model, maxit=20, grid=Dict(:tag => :Smolyak, :mu => 3))

        @test true
    end

    # kvec = linspace(dr.grid.min[1], dr.grid.max[1], 10)
    # nvec = [dr(1, [k])[1] for k in kvec]
    # ivec = [dr(1, [k])[2] for k in kvec]
    # nvec_d = [drd(1, [k])[1] for k in kvec]
    # ivec_d = [drd(1, [k])[2] for k in kvec]
    # nvec_0 = [dr0(1, [k])[1] for k in kvec]
    # ivec_0 = [dr0(1, [k])[2] for k in kvec]
    #
    # @assert maximum(abs, nvec_d-nvec)<1e-5
    # @assert maximum(abs, nvec_0-nvec)<1e-5 # not satisfied right now (see tol. of optimizer)

    # #TODO : reactivate
    # @testset "testing ar1 models" begin

    #     # AR1 model: this one should be exactly equivalent to rbc_dtcc_ar1
    #     import Dolo
    #     fn = joinpath(Dolo.pkg_path, "examples", "models", "rbc.yaml")
    #     model = Dolo.yaml_import(fn)
    #     dp = Dolo.discretize(model.exogenous)

    #     @time dr = Dolo.perturb(model)
    #     # @time tid_res = Dolo.time_iteration_direct(model; maxit=20, verbose=true)
    #     @time ti_res = Dolo.time_iteration(model; maxit=20, verbose=false)
    #     @time iti_res = Dolo.improved_time_iteration(model; maxit=20, verbose=false)
    #     @test Dolo.converged(iti_res)
    #     # @time sol_v = Dolo.value_iteration(model, ti_res.dr; maxit=3, verbose=true)

    #     model.symbols[:exogenous]
    #     #
    #     Dolo.simulate(model, ti_res.dr; m0=[0.015], N=10)
    #     Dolo.response(model.exogenous, [0.01])
    #     #
    #     sim = Dolo.response(model, ti_res.dr, :z)
    #     #
    #     irf = Dolo.response(model, ti_res.dr, :z, -0.01)
    #     @test true
    # end

end


text:
module SubPackage
using Pkg.Artifacts

# All this module will do is reference its `arty` Artifact.
arty = artifact"arty"
end


text:
export Embed
using Knet: atype
using AutoGrad: Param

"""
    Embed(vocabsize, embedsize; atype, init)
    Embed(weights)

References:
* torch.nn.Embed
* tf.keras.layers.Embed
"""
struct Embed; w;
    Embed(w)=new(w isa Param ? w : Param(w))
end

function Embed(vocabsize::Integer, embedsize::Integer; atype=atype(), init=𝑵(1))
    Embed(convert(atype, init(embedsize,vocabsize)))
end

function (l::Embed)(x)
    l.w[:,x]
end

# function (l::Embed)(x::MaskedArray)
#     a = l(x.array)
#     m = (x.mask === nothing ? nothing : reshape(x.mask, 1, size(x.mask)...))
#     MaskedArray(a, m)
# end


text:
module TestTimestepArrays

using Mimi
using Test

import Mimi:
    FixedTimestep, VariableTimestep, TimestepVector, TimestepMatrix, next_timestep, hasvalue,
    isuniform, first_period, last_period, first_and_step, reset_compdefs

reset_compdefs()

a = collect(reshape(1:16,4,4))

## quick check of isuniform
@test isuniform([]) == false
@test isuniform([1]) == true
@test isuniform([1,2,3]) == true
@test isuniform([1,2,3,5]) == false
@test isuniform(1) == true

@test first_and_step([2,3,4]) == (2,1)
@test first_and_step([1:2:10...]) == (1,2)

#------------------------------------------------------------------------------
# 1. Test TimestepVector - Fixed Timestep
#------------------------------------------------------------------------------
years = collect(2000:1:2003)

#1a.  test constructor, lastindex, and length (with both
# matching years and mismatched years)

x = TimestepVector{FixedTimestep{2000, 1}, Int}(a[:,3])
@test length(x) == 4
@test lastindex(x) == 4

#1b.  test hasvalue, getindex, and setindex! (with both matching years and
# mismatched years)

t = FixedTimestep{2001, 1, 3000}(1)

@test hasvalue(x, t)
@test !hasvalue(x, FixedTimestep{2000, 1, 2012}(10))
@test x[t] == 10

t2 = next_timestep(t)

@test x[t2] == 11

x[t2] = 99
@test x[t2] == 99

t3 = FixedTimestep{2000, 1, 2003}(1)
@test x[t3] == 9

x[t3] = 100
@test x[t3] == 100

x[1] = 101
@test x[1] == 101

#------------------------------------------------------------------------------
# 2. Test TimestepVector - Variable Timestep
#------------------------------------------------------------------------------

years = (2000, 2005, 2015, 2025)
x = TimestepVector{VariableTimestep{years}, Int}(a[:,3])

#2a.  test hasvalue, getindex, and setindex! (with both matching years and
# mismatched years)

y2 = Tuple([2005:5:2010; 2015:10:3000])
t = VariableTimestep{y2}()

@test hasvalue(x, t)
@test !hasvalue(x, VariableTimestep{years}(10))
@test x[t] == 10

t2 =  next_timestep(t)

@test x[t2] == 11
#@test indices(x) == (2000:2003,) #may remove this function

x[t2] = 99
@test x[t2] == 99

t3 = VariableTimestep{years}()
@test x[t3] == 9

x[t3] = 100
@test x[t3] == 100

x[1] = 101
@test x[1] == 101

#------------------------------------------------------------------------------
# 3. Test TimestepMatrix - Fixed Timestep
#------------------------------------------------------------------------------
years = Tuple(2000:1:2003)

#3a.  test constructor (with both matching years
# and mismatched years)

y = TimestepMatrix{FixedTimestep{2000, 1}, Int}(a[:,1:2])

#3b.  test hasvalue, getindex, and setindex! (with both matching years and
# mismatched years)

t = FixedTimestep{2001, 1, 3000}(1)

@test hasvalue(y, t, 1)
@test !hasvalue(y, FixedTimestep{2000, 1, 3000}(10), 1)
@test y[t,1] == 2
@test y[t,2] == 6

t2 = next_timestep(t)

@test y[t2,1] == 3
@test y[t2,2] == 7

y[t2, 1] = 5
@test y[t2, 1] == 5

t3 = FixedTimestep{2000, 1, 2005}(1)

@test y[t3, 1] == 1
@test y[t3, 2] == 5

y[t3, 1] = 10
@test y[t3,1] == 10

y[:,:] = 11
@test all([y[i,1] == 11 for i in 1:4])
@test all([y[1,j] == 11 for j in 1:2])

#3c.  interval wider than 1
z = TimestepMatrix{FixedTimestep{2000, 2}, Int}(a[:,3:4])
t = FixedTimestep{1980, 2, 3000}(11)

@test z[t,1] == 9
@test z[t,2] == 13

t2 = next_timestep(t)
@test z[t2,1] == 10
@test z[t2,2] == 14

#------------------------------------------------------------------------------
# 4. Test TimestepMatrix - Variable Timestep
#------------------------------------------------------------------------------

years = Tuple([2000:5:2005; 2015:10:2025])
y = TimestepMatrix{VariableTimestep{years}, Int}(a[:,1:2])

#4a.  test hasvalue, getindex, setindex!, and lastindex (with both matching years and
# mismatched years)

t = VariableTimestep{Tuple([2005:5:2010; 2015:10:3000])}()

@test hasvalue(y, t, 1)
@test !hasvalue(y, VariableTimestep{years}(10))
@test y[t,1] == 2
@test y[t,2] == 6

t2 = next_timestep(t)

@test y[t2,1] == 3
@test y[t2,2] == 7

y[t2, 1] = 5
@test y[t2, 1] == 5

t3 = VariableTimestep{years}()

@test y[t3, 1] == 1
@test y[t3, 2] == 5

y[t3, 1] = 10
@test y[t3,1] == 10

y[:,:] = 11
@test all([y[i,1] == 11 for i in 1:4])
@test all([y[1,j] == 11 for j in 1:2])

#------------------------------------------------------------------------------
# 5. Test TimestepArray methods
#------------------------------------------------------------------------------

x_years = Tuple(2000:5:2015) #fixed
y_years = Tuple([2000:5:2005; 2015:10:2025]) #variable

x_vec = TimestepVector{FixedTimestep{2000, 5}, Int}(a[:,3])
x_mat = TimestepMatrix{FixedTimestep{2000, 5}, Int}(a[:,1:2])
y_vec = TimestepVector{VariableTimestep{y_years}, Int}(a[:,3])
y_mat = TimestepMatrix{VariableTimestep{y_years}, Int}(a[:,1:2])

@test first_period(x_vec) == first_period(x_mat) == x_years[1]
@test first_period(y_vec) == first_period(y_mat) == y_years[1]
@test last_period(x_vec) == last_period(x_mat) == x_years[end]
@test last_period(y_vec) == last_period(y_mat) == y_years[end]

@test size(x) == size(a[:,3])
@test size(y) == size(a[:,1:2])
@test size(y,2) == size(a[:,1:2],2)

@test ndims(x) == 1
@test ndims(y) == 2

@test eltype(x) == eltype(a)
@test eltype(y) == eltype(a)

fill!(x, 2)
fill!(y, 2)
@test x.data == fill(2, (4))
@test y.data == fill(2, (4, 2))

#------------------------------------------------------------------------------
# 6. Test that getindex for TimestepArrays doesn't allow access to `missing`
#       values during `run` that haven't been computed yet.
#------------------------------------------------------------------------------

@defcomp foo begin
    par = Parameter(index=[time])
    var = Variable(index=[time])
    function run_timestep(p, v, d, t)
        if is_last(t)
            v.var[t] = 0
        else
            v.var[t] = p.par[t+1]   # This is where the error will be thrown, if connected to an internal variable that has not yet been computed.
        end
    end
end

years = 2000:2010

m = Model()
set_dimension!(m, :time, years)
add_comp!(m, foo, :first)
add_comp!(m, foo, :second)
connect_param!(m, :second=>:par, :first=>:var)
set_param!(m, :first, :par, 1:length(years))

@test_throws ErrorException run(m)

end #module


text:
using Forge.DSL

window = create_window(1024, 768, "3d Surface Demo")
show_window(window)

const XMIN = -32
const XMAX =  32
const YMIN = -32
const YMAX =  32
const DX = 0.25
const XSIZE = (XMAX-XMIN)/DX
const YSIZE = (YMAX-YMIN)/DX

chart = create_chart_3d()
set_axes_limits(chart, XMIN-2, XMAX+2, YMIN-2, YMAX+2, -0.5, 1)
set_axes_titles(chart, "x-axis", "y-axis", "z-axis")

surface = create_surface(chart, XSIZE, YSIZE)

data = []
for x in XMIN:DX:XMAX+DX, y in YMIN:DX:YMAX+DX
    z = sqrt(x*x + y*y) + 2.2204e-16
    push!(data, x, y, sin(z)/z)
end

update_vertex_buffer(surface, data)

draw_chart(!isclosed, window, chart)
save_framebuffer(window, "surface.png")