text:
#####
##### Building blocks for ordering columns and tables.
#####

export TrustOrdering, TryOrdering, VerifyOrdering

####
#### ColumnOrder
####

"""
$(TYPEDEF)

Ordering specification for a column. `K::Symbol` is a key for ordering by `isless`,
`R::Bool == true` reverses that for this key. *Internal.*
"""
struct ColumnOrdering{K, R}
    function ColumnOrdering{K, R}() where {K, R}
        @argcheck K isa Symbol
        @argcheck R isa Bool
        new{K, R}()
    end
end

"""
$(SIGNATURES)

Accessor for key. *Internal.*
"""
orderkey(co::ColumnOrdering{K}) where {K} = K

orderrev(co::ColumnOrdering{K, R}) where {K, R} = R

"""
$(SIGNATURES)

Process a column ordering specification, called by [`table_ordering`](@ref).
"""
@inline ColumnOrdering(key::Symbol, rev::Bool = false) = ColumnOrdering{key, rev}()
@inline ColumnOrdering(keyrev::Pair{Symbol, typeof(reverse)}) = ColumnOrdering(first(keyrev), true)
@inline ColumnOrdering(cs::ColumnOrdering) = cs
ColumnOrdering(x) = throw(ArgumentError("Unrecognized ordering specification $(x)."))

"Types accepted as valid column ordering specifications in the user API."
const ColumnOrderingSpecification =
    Union{Symbol, Pair{Symbol, typeof(reverse)}, ColumnOrdering}

ordering_repr(cs::ColumnOrdering{K, R}) where {K, R} = (R ? "↓" : "↑") * string(K)

"""
Shorthand for table ordering. *Internal.*
"""
const TableOrdering = Tuple{Vararg{ColumnOrdering}}

"""
$(SIGNATURES)

A string representation of an ordering, eg for use in `show`.
"""
function ordering_repr(ordering::TableOrdering)
    isempty(ordering) ? "no ordering" : "ordering " * join(ordering_repr.(ordering), " ")
end

"""
$(SIGNATURES)

Process ordering specifications for columns (an iterable or possibly a TableOrdering),
return a `TableOrdering`. Check for uniqueness (but not validity) of keys.

Accepted syntax:

- `:key`, for ordering a column in ascending order,

- `:key => reverse`, for ordering a column in descending order.

All functions which accept ordering specifications should use this, but the function itself
is not part of the API.
"""
table_ordering(column_ordering_specifications) =
    map(ColumnOrdering, tuple(column_ordering_specifications...))

"""
$(SIGNATURES)

When `invert == false`, keep the initial part of ordering that has keys in `keys`. Not
having a key in `keys` invalidates the tail ordering from that point. This is useful for
selecting subsets of orderings.

When `invert == true`, *having a key* in `keys` invalidates the ordering. This is useful for
orderings of merged and dropped columns.
"""
Base.@pure function mask_ordering(ordering::TableOrdering, keys::Keys, invert::Bool = false)
    kept = ColumnOrdering[]
    for o in ordering
        if key_in(orderkey(o), keys) ⊻ invert
            push!(kept, o)
        else
            # skipping a ColumnOrdering invalidates the rest
            break
        end
    end
    (kept..., )
end

"""
$(SIGNATURES)

Extend `ordering` so that it can be used to split a table with `splitkeys`.
"""
function split_compatible_ordering(ordering::TableOrdering, splitkeys::Keys)
    @argcheck allunique(splitkeys)
    o_first, o_rest = first(ordering), Base.tail(ordering)
    s_first = first(splitkeys)
    if orderkey(o_first) ≡ s_first
        (o_first, split_compatible_ordering(o_rest, Base.tail(splitkeys))...)
    else
        o_matched = findfirst(o -> orderkey(o) ≡ s_first, o_rest)
        (o_matched ≡ nothing ? ColumnOrdering{s_first, false}() : o_rest[o_matched],
         split_compatible_ordering(ordering, Base.tail(splitkeys))...)
    end
end

split_compatible_ordering(ordering::TableOrdering, splitkeys::Tuple{}) = ()

split_compatible_ordering(ordering::Tuple{}, splitkeys::Tuple{}) = ()

split_compatible_ordering(ordering::Tuple{}, splitkeys::Keys) = table_ordering(splitkeys)

####
#### Comparisons
####

cmp_ordering(::ColumnOrdering{K, R}, a, b) where {K, R} =
    cmp(getproperty(a, K), getproperty(b, K)) * (R ? -1 : 1)

"""
$(SIGNATURES)

Compare `a` and `b` given the `ordering`.

*Internal*.
"""
@inline cmp_ordering(ordering::TableOrdering, a, b) = _cmp_table_ordering(a, b, ordering...)

_cmp_table_ordering(a, b) = 0

function _cmp_table_ordering(a, b, column_ordering::ColumnOrdering, rest...)
    r = cmp_ordering(column_ordering, a, b)
    r ≠ 0 && return r
    _cmp_table_ordering(a, b, rest...)
end

@inline isless_ordering(ordering, a, b) = cmp_ordering(ordering, a, b) == -1

"""
$(SIGNATURES)

Return the (initial) part of `ordering` under which `a ≤ b`.
"""
retained_ordering(ordering::TableOrdering, a, b) = _retained_table_ordering(a, b, ordering...)

_retained_table_ordering(a, b) = ()

function _retained_table_ordering(a, b, column_ordering::ColumnOrdering, rest...)
    if cmp_ordering(column_ordering, a, b) ≤ 0
        (column_ordering, _retained_table_ordering(a, b, rest...)...)
    else
        ()
    end
end

####
#### Ordering rules
####

"""
$(TYPEDEF)

Rule for dealing with specified orderings. Verifies key uniqueness.

See [`VerifyOrdering`](@ref), [`TrustOrdering`](@ref), and [`TryOrdering`](@ref).

This type and its methods are *internal*.
"""
struct OrderingRule{R, O <: TableOrdering}
    ordering::O
    function OrderingRule{R}(ordering::O) where {R, O <: TableOrdering}
        @argcheck R ∈ (:trust, :verify, :try)
        @argcheck allunique(orderkey.(ordering)) "Duplicate order keys."
        new{R, O}(ordering)
    end
end

OrderingRule{R}(column_ordering_specifications::Tuple{Vararg{ColumnOrderingSpecification}}
                ) where {R} = OrderingRule{R}(table_ordering(column_ordering_specifications))

OrderingRule{R}(column_ordering_specifications::ColumnOrderingSpecification...) where {R} =
    OrderingRule{R}(column_ordering_specifications)

"Verify that the specified ordering holds. This is the default ordering rule."
const VerifyOrdering = OrderingRule{:verify}

"""
Accept the specified ordering to hold without any checks (except for verifying that column
names are valid).

!!! note
    This can lead to incorrect results, use cautiously. [`VerifyOrdering`](@ref) is
    recommended instead as it has little overhead.
"""
const TrustOrdering = OrderingRule{:trust}

"Try the specified ordering, when necesary degrade to a subset of it that holds."
const TryOrdering = OrderingRule{:try}

# forced conversion
OrderingRule{R}(ordering_rule::OrderingRule) where {R} =
    OrderingRule{R}(ordering_rule.ordering)

"""
$(SIGNATURES)

For `TryOrdering`, return a masked ordering rule that is is defined on `keys` so that
comparisons make sense, otherwise return the original `ordering_rule` (that will just error
for undefined keys).
"""
Base.@pure function mask_try_ordering(ordering_rule::OrderingRule{R}, keys::Keys) where {R}
    R ≡ :try ? OrderingRule{R}(mask_ordering(ordering_rule.ordering, keys)) : ordering_rule
end


text:
# tests MNL loss

srand(1);
m,n,k = 200,50,2;
kfit = k+1
K = 4; # number of categories
d = n*K;
# matrix to encode
X_real, Y_real = randn(m,k), randn(k,d);
XY = X_real*Y_real;
# subtract the mean so we can compare the truth with the fit;
# the loss function is invariant under shifts
losses = fill(MultinomialLoss(K),n)
yidxs = get_yidxs(losses)
for i=1:m
	for j=1:n
		mef = mean(XY[i,yidxs[j]])
		XY[i,yidxs[j]] = XY[i,yidxs[j]] - mef
	end
end

A = zeros(Int, (m, n))
for i=1:m
	for j=1:n
		wv = WeightVec(Float64[exp(-XY[i, K*(j-1) + l]) for l in 1:K])
		l = sample(wv)
		A[i,j] = l
	end
end

# and the model
losses = Loss[fill(MultinomialLoss(K),n)...]
rx, ry = QuadReg(), QuadReg();
glrm = GLRM(A,losses,rx,ry,kfit, scale=false, offset=false, X=randn(kfit,m), Y=randn(kfit,d));
gg = GGLRM(A, losses, rx, ry, kfit, X=randn(kfit,m), Y=randn(kfit,d))
fit!(glrm);
fit!(gg);
@test_approx_eq_eps(whole_objective(gg, gg.X'gg.Y), whole_objective(gg, glrm.X'glrm.Y), 0.1*length(A))

ext: jl


text:
using QuickTypes
using QuickTypes: construct, roottypeof, fieldsof, type_parameters, roottype,
      tuple_parameters, @d
using Test
using ConstructionBase: setproperties

abstract type Vehicle end

@qstruct Car{T<:Number, U}(size::T, nwheels::Int=4; manufacturer::U=nothing,
                           brand::String="off-brand") <: Vehicle

c = Car(10; manufacturer=("Danone", "Hershey"))

@test c.size==10
@test c.nwheels==4
@test c.manufacturer==("Danone", "Hershey")
@test c.brand=="off-brand"
c2 = @inferred setproperties(c, (size=42, nwheels=8))
@test c2.nwheels == 8
@test c2.size == 42
@test c2.brand == c.brand
# Check that the fields are in the right order
@test collect(fieldnames(Car)) == [:size, :nwheels, :manufacturer, :brand]
# This is essentially the definition of these functions.
@test construct(roottypeof(c), fieldsof(c)...) == c
@test type_parameters(Vector{Int}) == Base.Core.svec(Int64, 1)
@test tuple_parameters(Tuple{Int, Float64}) == Base.Core.svec(Int64, Float64)
@inferred roottypeof(1=>2) == Pair

################################################################################

@qstruct Empty()
Empty()
@test setproperties(Empty(), NamedTuple()) === Empty()

# Used to yield:
#     WARNING: static parameter T does not occur in signature for Type.
#     The method will not be callable.
@qstruct Blah{T}()

################################################################################

@qstruct Boring(x::Int)
@inferred Boring(10)
@test Boring(10).x == 10
@test Boring(10.0).x == 10   # check that convert is called correctly
@qstruct ParametricBoring{X}(x::X; _concise_show=true)
@inferred ParametricBoring(10)
@test ParametricBoring(10).x === 10
o = ParametricBoring(1)
@test setproperties(o, x=:one).x === :one

@qstruct Kwaroo(x; y=10)
@test Kwaroo(5) == Kwaroo(5; y=10)
o = Kwaroo(5, y=10)
o2 = @inferred setproperties(o, (x=:five, y=100.0))
@test o2 isa Kwaroo
@test o2.x === :five
@test o2.y === 100.0

################################################################################
# Slurping

@qstruct Slurp(x, y=1, args...; kwargs...)
s = Slurp(1,2,3,4,5,6,7; x=1, y=10+2)
@test s.args == (3,4,5,6,7)
@test s.kwargs == pairs((x=1, y=12))
s2 = @inferred setproperties(s, x=:hello)
@test s2 isa Slurp
@test s2.x == :hello
@test s2.y == s.y

let
    @unpack_Slurp Slurp(10)
    @test x == 10
    @test y == 1
end

# Slurping (parametric case)

@qstruct SlurpParam{T}(x::AbstractVector{T}, y=1, args...; kwargs...)
s = SlurpParam([1,2,3,4,5,6,7], 8, 9,10; x=1, y=10+2)
@test s.args == (9,10)
@test s.kwargs == pairs((x=1, y=12))

################################################################################

@qmutable Foo{T}(x::T; y=2) do
    @assert x < 10
end

@test_throws AssertionError Foo(11; y=10.0)
@test_throws AssertionError construct(Foo, 11, 10.0)

################################################################################
# Fully-parametric

@qstruct_fp Plane(nwheels::Number; brand=:zoomba) do
    @assert nwheels < 100
end <: Vehicle

@test_throws MethodError Plane{Int, Symbol}(2; brand=12)
@test Plane{Int, Symbol}(2; brand=:zoomba).brand == :zoomba
@test supertype(Plane) == Vehicle
# This used to be a MethodError, but since we moved the outer constructor inside
# the type, it has become a TypeError. Not sure why!
@test_throws TypeError Plane("happy")

@qstruct_fp NoFields()   # was an error before it was special-cased

o = Plane(4)
o2 = @inferred setproperties(o, brand=10, nwheels=o.nwheels)
@test o2 isa Plane
@test o2.brand === 10
@test o2.nwheels === o.nwheels

################################################################################
# Narrowly-parametric

@qstruct_fp Foo_fp(a, b)
@qstruct_np Foo_np(a, b)
convert_f(foo) = convert(foo.a, 10)
@test_throws(Exception, @inferred convert_f(Foo_fp(Int, 2)))
@inferred convert_f(Foo_np(Int, 2))
@test fieldtype(typeof(Foo_np(Int, 2)), :a) == Type{Int64}

@qstruct Issue11(;no_default_value)
@test_throws UndefKeywordError Issue11()

################################################################################
# @qfunctor

@qfunctor function Action(a; kw=100)(x)
    return a + x + kw
end

@test Action(2)(10) == 112

@qfunctor ParamAction{X}(a::X)(b::T) where T = (a, b, X, T)
@test ParamAction(1)(2.0) == (1, 2.0, Int, Float64)

################################################################################
# @destruct

@destruct foo(Ref(x)) = x+2
@destruct foo(Ref{Float64}(x)) = x+10
@test foo(Ref(10)) == 12
@test foo(Ref(10.0)) == 20
@destruct foo(a, (Ref{T} where T)(x)) = a + x

struct LongerStruct{X}
    a
    b
    c::X
end

@destruct function kwfun(LongerStruct{X}(u,v; c, bof=b)) where X
    return u, v, c, bof
end
@test kwfun(LongerStruct(4,5,6)) == (4,5,6,5)

@destruct nested(LongerStruct(Ref(Ref(a)))) = a
@test nested(LongerStruct(Ref(Ref(44)), 3, 4)) == 44

@destruct tup_destruct(Ref((a,Ref(b)))) = (a, b)
@test tup_destruct(Ref((1,Ref(2)))) == (1,2)

@d Ref(x) := Ref(111)
@test x == 111

@destruct for (LongerStruct(Ref(xx)), Ref(yy)) in [(LongerStruct(Ref(55), 10, 20), Ref(66))]
    @test (xx, yy) == (55, 66)
end

@d LongerStruct(x)(y) = (x, y)
@test LongerStruct(10, 20, 30)(5) == (10, 5)

@test @d((Ref(x), Ref(y)) -> x+2)(Ref(10), Ref(20)) == 12

@d with_type(Ref(a::Int)) = a
@test with_type(Ref(1)) === 1
@test with_type(Ref(2.0)) === 2   # This is not necessarily _desirable_, but it's reasonable...


text:
@testset "155.min-stack.jl" begin
    @testset "155-test case 1" begin
        stack = MinStack()
        @test push!(stack, -2) === nothing
        @test push!(stack, 0) === nothing
        @test push!(stack, -3) === nothing
        @test get_min(stack) == -3
        @test pop!(stack) === nothing
        @test top(stack) == 0
        @test get_min(stack) == -2
    end
    @testset "155-test case 2" begin
        stack = MinStack()
        push!(stack, 2) === nothing
        @test push!(stack, 0) === nothing
        @test push!(stack, 3) === nothing
        @test push!(stack, 0) === nothing
        @test get_min(stack) == 0
        @test pop!(stack) === nothing
        @test get_min(stack) == 0
        @test pop!(stack) === nothing
        @test get_min(stack) == 0
        @test pop!(stack) === nothing
        @test get_min(stack) == 2
    end
end


text:
"""
    NodeType

Each node in the directed graph can be classified into the following types
- VARIABLE: Denotes a decision variable.
- PARAMETER: An adjustable parameter value (not a decision variable).
- CONSTANT: A constant value
- EXPRESSION: Any other expression that isn't a subexpression
- SUBEXPRESSION: Any expression referencing a different graph representation.
"""
@enum(NodeClass, CONSTANT, PARAMETER, VARIABLE, EXPRESSION, SUBEXPRESSION)
#=
- SELECT: These nodes do not store values rather they reference the output
          stored in the buffer assocatied with an EXPRESSION. These nodes
          only occur after a multiple output function. In the mimo JuMP
          extension they correspond to intermediate variables introduced
          by the user.
=#

struct Variable end
struct Parameter end
struct Constant end
struct Expression end
struct Subexpression end
struct User end

function user end
function usern end

const ATOM_TYPE_INSTANCES = instances(AtomType)

abstract type AbstractNode end

"""
$(TYPEDEF)

Describes connectivity and expression represented by node.
"""
struct Node <: AbstractNode
    node_class::NodeClass
    ex_type::AtomType
    first_index::Int
    second_index::Int
    arity::Int
    children::Vector{Int}
end

for (t, s, a) in ((Variable, VARIABLE, VAR_ATOM),
                  (Parameter, PARAMETER, PARAM_ATOM),
                  (Constant, CONSTANT, CONST_ATOM),
                  (Subexpression, SUBEXPRESSION, SUBEXPR),)
    @eval Node(::$t, i::Int) = Node($s, $a, i, 0, 0, Int[])
end

for v in (PLUS, MINUS, MULT, POW, DIV, MAX, MIN)
    @eval Node(::Val{true}, ::Val{$v}, c::Vector{Int}) = Node(EXPRESSION, $v, 0, 0, length(c), c)
end
for d in ALL_ATOM_TYPES
    @eval Node(::Val{false}, ::Val{$d}, c::Vector{Int}) = Node(EXPRESSION, $d, 0, 0, 1, c)
end
Node(::Val{true}, ::Val{USER}, i::Int, c::Vector{Int}) = Node(EXPRESSION, USER, i, 0, 1, c)
Node(::Val{true}, ::Val{USERN}, i::Int, c::Vector{Int}) = Node(EXPRESSION, USERN, i, 0, length(c), c)

Node(::Val{:first_index}, n::Node, i::Int) = Node(node_class(n), ex_type(n), i, second_index(n), arity(n), children(n))

node_class(n::Node)   = n.node_class
ex_type(n::Node)      = n.ex_type
first_index(n::Node)  = n.first_index
second_index(n::Node) = n.node_second_index
arity(n::Node)        = n.arity
children(n::Node)     = n.children
child(n::Node, i)     = @inbounds getindex(n.children, i)

node_is_class(::Variable, n::Node) = node_class(n) == VARIABLE
node_is_class(::Parameter, n::Node) = node_class(n) == PARAMETER
node_is_class(::Constant, n::Node) = node_class(n) == CONSTANT

mv_eago_not_jump = setdiff(JuMP._Derivatives.operators,
                           union(Symbol[k for k in keys(REV_BIVARIATE_ATOM_DICT)],
                                 Symbol[k for k in keys(REV_NARITY_ATOM_DICT)]))
eago_mv_switch = quote end
for s in mv_eago_not_jump
    global eago_mv_switch = quote
        $eago_mv_switch
        (d == $s) && (return Node(Val(true), Val($s), v[c]))
    end
end
@eval function _create_call_node(i, v, c::UnitRange{Int}, op::OperatorRegistry)
    if i == 1
        return Node(Val(true), Val(PLUS), v[c])
    elseif i == 2
        return Node(Val(true), Val(MINUS), v[c])
    elseif i == 3
        return Node(Val(true), Val(MULT), v[c])
    elseif i == 4
        return Node(Val(true), Val(POW), v[c])
    elseif i == 5
        return Node(Val(true), Val(DIV), v[c])
    elseif i == 6
        error("If-else currently unsupported...")
    elseif i >= JuMP._Derivatives.USER_OPERATOR_ID_START
        i_mv = i - JuMP._Derivatives.USER_OPERATOR_ID_START + 1
        d = op.multivariate_id[i_mv]
        $eago_mv_switch
        return Node(Val(true), Val(USERN), i_mv, v[c])
    end
end

function binary_switch_typ(ids, exprs)
    if length(exprs) <= 3
        out = Expr(:if, Expr(:call, :(==), :i, ids[1]),
                   :(Node(Val(false), Val($(exprs[1])), v[c])))
        if length(exprs) > 1
            push!(out.args, binary_switch_typ(ids[2:end], exprs[2:end]))
        end
        return out
    else
        mid = length(exprs) >>> 1
        return Expr(:if, Expr(:call, :(<=), :i, ids[mid]),
                         binary_switch_typ(ids[1:mid], exprs[1:mid]),
                         binary_switch_typ(ids[mid+1:end], exprs[mid+1:end]))
    end
end

indx_JuMP = Int[]
indx_EAGO = AtomType[]
for k in univariate_operators
    if haskey(REV_UNIVARIATE_ATOM_DICT, k)
        k_EAGO = REV_UNIVARIATE_ATOM_DICT[k]
        push!(indx_JuMP, univariate_operator_to_id[k])
        push!(indx_EAGO, k_EAGO)
    end
end

uni_eago_not_jump = setdiff(univariate_operators, Symbol[k for k in keys(REV_UNIVARIATE_ATOM_DICT)])
uni_eago_not_jump = push!(uni_eago_not_jump, :-)
eago_uni_switch = quote end
for s in uni_eago_not_jump
    global eago_uni_switch = quote
        $eago_uni_switch
        (d == $s) && (return Node(Val(false), Val($s), v[c]))
    end
end
atom_switch = binary_switch_typ(indx_JuMP, indx_EAGO)
@eval function _create_call_node_uni(i::Int, v, c::UnitRange{Int}, op::OperatorRegistry)

    if i >= JuMP._Derivatives.USER_UNIVAR_OPERATOR_ID_START
        j = i - JuMP._Derivatives.USER_UNIVAR_OPERATOR_ID_START + 1
        dop = op.univariate_operator_id[j]
        d = op.univariate_operator_to_id[dop]
        $eago_uni_switch
        return Node(Val(true), Val(USER), j, v[c])
    end
    $atom_switch
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using SmartAsserts: SmartAsserts, @smart_assert, @show_assertion, group_function_args
using Test

@testset "SmartAsserts.jl" begin

    @test group_function_args(:(f(x1; kw1=1, kw2, a.kw3, kw4=2)).args) ==
        (:f, Any[:x1], Pair[:kw1 => 1, :kw2 => :kw2, :kw3 => :(a.kw3), :kw4 => 2])

    @test begin
        @smart_assert 1 < 2
        @smart_assert 1 == 1 "This should pass"
        @smart_assert ≈(1, 1.2, atol=0.4) "This should pass"
        true
    end

    @test_throws AssertionError begin
        a = 5
        @smart_assert a < 1
    end

    @test_throws AssertionError begin
        a = 5
        @smart_assert a < 1 "This should fail. 2a = $(2 * a)"
    end

    @test_throws AssertionError begin
        @smart_assert ≈(1, 1.2, atol=0.1; rtol=0.1)
    end

    @test begin
        a = 5
        occursin("`a` evaluates to 5", @show_assertion(a < 1)[2])
    end

    @test begin
        occursin(
            "`reverse([1, 1, 2])` evaluates to [2, 1, 1]",
            @show_assertion(allunique(reverse([1, 1, 2])))[2],
        )
    end

    @test begin
        f_with_kw(a; b, c) = a == b + c
        c = 2
        reason = @show_assertion(f_with_kw(4, b=1; c))[2]
        !occursin("`4` evaluates to 4", reason) &&
            occursin("`c` evaluates to 2", reason) &&
            !occursin("`1` evaluates to 1", reason)
    end

    @test begin
        a = 1
        b = 3
        reason = @show_assertion(a <= b <= (b - 1))[2]
        occursin("`a` evaluates to 1", reason) &&
            occursin("`b` evaluates to 3", reason) &&
            occursin("`b - 1` evaluates to 2", reason)
    end

    @test_throws AssertionError begin
        @smart_assert 1 + 1 == 0 "Should fail"
    end
end

SmartAsserts.set_enabled(false)
@testset "Turning off" begin
    @test begin
        @smart_assert 1 + 1 == 0 "Should be turned off"
        true
    end
end

SmartAsserts.set_enabled(true)
@testset "Turning back on" begin
    @test_throws AssertionError begin
        @smart_assert 1 + 1 == 0 "Should fail again"
    end
end


text:
@testset "syntax error reports" begin
    let
        s = """
        begin
            c = rand(Bool)
        end

        end
        """

        res, interp = profile_text′(s)

        @test !isempty(res.toplevel_error_reports)
        @test first(res.toplevel_error_reports) isa SyntaxErrorReport
    end
end

@testset "fix self-reference of virtual module" begin
    let
        res, interp = @profile_toplevel begin
            const foo = sum
            Main.foo("julia") # `Main.sum` should be resolved as constant
        end
        test_sum_over_string(res)
    end

    let
        res, interp = @profile_toplevel begin
            let
                Main = "julia" # local `Main` should not be resolved to virtual module
                sum(Main)
            end
        end
        test_sum_over_string(res)
    end
end

@testset "fix toplevel global `Symbol`" begin
    # this case otherwise will throw an error in `CC.typ_for_val` in optimization
    let
        res, interp = @profile_toplevel begin
            v = '1'
            v = if rand(Bool)
                rand(Int)
            else
                v
            end
            sin(v)
        end
        @test true
    end

    # `c` wrapped in `GotoIfNot` node should be transformed into `GlobalRef(vmod, :c)`,
    # otherwise `NonBooleanCondErrorReport` (and `UncaughtExceptionReport`) will be reported
    let
        res, interp = @profile_toplevel begin
            c = false
            if c
                throw("should be ignored")
            end
        end
        @test isempty(res.inference_error_reports)
    end
end

@testset "\"toplevel definitions\"" begin
    let
        vmod = gen_virtual_module()
        @profile_toplevel vmod begin
            # function
            foo() = nothing

            # local function, shouldn't be evaluated
            let
                c = rand(Bool)
                foo′(a) = a || c
            end

            # macro
            macro foo(ex) ex end

            # abstract type
            abstract type Foo end

            # struct
            struct Foo1 <: Foo
                val::Int
            end
            mutable struct Foo2 <: Foo
                val::Int
            end

            # primitive type
            primitive type Foo3 32 end

            # import, using
            using Base: Fix1
            import Base: getproperty

            function getproperty(foo::Foo, sym::Symbol)
                return if sym === :val
                    getfield(foo, sym)::Int
                else
                    getfield(foo, sym)
                end
            end
        end

        @test is_concrete(vmod, :foo)
        @test !is_concrete(vmod, :foo′)
        @test is_concrete(vmod, Symbol("@foo"))
        @test is_concrete(vmod, :Foo)
        @test is_concrete(vmod, :Foo1)
        @test is_concrete(vmod, :Foo2)
        @test is_concrete(vmod, :Foo3)
        @test is_concrete(vmod, :Fix1)
        @test !isempty(methodswith(getfield(vmod, :Foo), getproperty))
    end

    # basic profiling with user-defined types
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            gb = rand(Bool)

            struct Foo
                b::Bool
                Foo(b = gb) = new(b)
            end

            foo = Foo(gb)
        end

        # global variables aren't evaluated but kept in `interp` instead
        @test is_abstract(vmod, :gb)
        @test isa_abstract(vmod.gb, Bool)
        @test is_concrete(vmod, :Foo)
        @test is_abstract(vmod, :foo)
        @test isa_abstract(vmod.foo, vmod.Foo)
        @test isempty(res.toplevel_error_reports)
        @test isempty(res.inference_error_reports)
    end

    # definitions using type aliases
    let
        vmod = gen_virtual_module()
        @profile_toplevel vmod begin
            const BT = Bool

            struct Foo
                b::BT
            end

            const c = Foo(rand(Bool))
            uc = Foo(rand(Bool))
        end
        @test is_concrete(vmod, :BT)
        @test is_concrete(vmod, :Foo)
        @test is_abstract(vmod, :c)
        @test is_abstract(vmod, :uc)
    end

    # a toplevel definition within a block
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            begin
                struct Foo
                    bar
                end

                foo = Foo(:bar)
                println(foo)
            end
        end
        @test is_concrete(vmod, :Foo)
        @test isa_abstract(vmod.foo, vmod.Foo)
    end

    # toplevel definitions within a block
    # somewhat related upstream issue: https://github.com/JuliaDebug/LoweredCodeUtils.jl/issues/47
    # well, the actual error here is world age error ...
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            begin
                abstract type Foo end
                struct Foo1 <: Foo
                    foo
                end

                foo = Foo1(:foo)
                println(foo)
            end
        end
        @test isempty(res.toplevel_error_reports)
        @test is_abstract(vmod, :foo)
    end

    @testset "toplevel definitions by `eval` calls" begin
        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                # these definitions shouldn't be abstracted away
                for fname in (:foo, :bar, :baz)
                    @eval begin
                        @inline ($(Symbol("is", fname)))(a) = a === $(QuoteNode(fname))
                    end
                end

                # these should be abstracted away (i.e. shouldn't throw)
                isfoo(:foo) && throw("foo")                   # should be reported
                isbar(:foo) && isbaz(:baz) && throw("foobaz") # shouldn't be reported
                isbar(:bar) && isbaz(:baz) && throw("barbaz") # should be reported
            end

            @test is_concrete(vmod, :isfoo)
            @test is_concrete(vmod, :isbar)
            @test is_concrete(vmod, :isbaz)
            @test length(res.inference_error_reports) == 2
            @test all(er->isa(er, UncaughtExceptionReport), res.inference_error_reports)
        end
    end
end

@testset "macro expansions" begin
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            @inline foo(a) = identity(a)
            foo(10)
        end

        @test is_concrete(vmod, :foo)
        @test isempty(res.toplevel_error_reports)
        @test isempty(res.inference_error_reports)
    end

    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            macro foo(ex)
                @assert Meta.isexpr(ex, :call)
                push!(ex.args, 1)
                return ex
            end

            @foo sin() # otherwise NoMethodError
        end

        @test is_concrete(vmod, Symbol("@foo"))
        @test isempty(res.toplevel_error_reports)
        @test isempty(res.inference_error_reports)
    end

    # macro expansions with access to global variables will fail
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            const arg = rand(Bool)

            macro foo(ex)
                @assert Meta.isexpr(ex, :call)
                arg && push!(ex.args, 1)
                return ex
            end

            @foo sin()
        end

        @test is_concrete(vmod, Symbol("@foo"))
        @test_broken isempty(res.toplevel_error_reports)
        @test isempty(res.inference_error_reports)
    end
end

@testset "remove `const`" begin
    let
        res, interp = @profile_toplevel begin
            const s = "julia"
            sum(s)
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end

    let
        res, interp = @profile_toplevel begin
            let
                const s = "julia"
                sum(s)
            end
        end

        @test !isempty(res.toplevel_error_reports)
        @test first(res.toplevel_error_reports) isa SyntaxErrorReport
    end
end

@testset "expression flattening" begin
    let
        @test (begin
            # internal error shouldn't occur
            res, interp = @profile_toplevel begin
                r = rand(); s = sin(a); c = cos(b); tan(c)
            end
        end; true)
    end
end

@testset "handle `include`" begin
    let
        f1 = normpath(FIXTURE_DIR, "include1.jl")
        f2 = normpath(FIXTURE_DIR, "include1.jl")

        vmod = gen_virtual_module()
        res, interp = profile_file′(f1, vmod)

        @test f1 in res.included_files
        @test f2 in res.included_files
        @test is_concrete(vmod, :foo)
        @test isempty(res.toplevel_error_reports)
        @test isempty(res.inference_error_reports)
    end

    let
        f = normpath(FIXTURE_DIR, "nonexistinclude.jl")
        res, interp = profile_file′(f)

        @test f in res.included_files
        @test !isempty(res.toplevel_error_reports)
        @test first(res.toplevel_error_reports) isa ActualErrorWrapped
        @test !isempty(res.inference_error_reports)
        @test first(res.inference_error_reports) isa GlobalUndefVarErrorReport
    end

    let
        f = normpath(FIXTURE_DIR, "selfrecursiveinclude.jl")
        res, interp = profile_file′(f)

        @test f in res.included_files
        @test !isempty(res.toplevel_error_reports)
        @test first(res.toplevel_error_reports) isa RecursiveIncludeErrorReport
    end

    let
        f1 = normpath(FIXTURE_DIR, "chainrecursiveinclude1.jl")
        f2 = normpath(FIXTURE_DIR, "chainrecursiveinclude2.jl")
        res, interp = profile_file′(f1)

        @test f1 in res.included_files
        @test f2 in res.included_files
        @test !isempty(res.toplevel_error_reports)
        let
            report = first(res.toplevel_error_reports)
            @test report isa RecursiveIncludeErrorReport
            @test report.duplicated_file == f1
            @test f1 in report.files
            @test f2 in report.files
            @test report.file == f2
            @test report.line == 1
        end
    end
end

@testset "module definition" for s in ("""
                                       module foo end
                                       """,
                                       """
                                       baremodule foo end
                                       """,
                                       """
                                       module foo
                                       module bar end
                                       end
                                       """
                                       )
    res, interp = profile_text′(s)
    @test isempty(res.toplevel_error_reports)
end

@testset "module usage" begin
    # using
    let
        res, interp = @profile_toplevel begin
            module foo

            using Base.Meta: isexpr

            isexpr(:(foo(bar)), :call)
            isexpr2(:(foo(bar)), :call)

            end
        end

        @test isempty(res.toplevel_error_reports)
        @test !isempty(res.inference_error_reports)
        @test first(res.inference_error_reports) isa GlobalUndefVarErrorReport
        @test occursin("isexpr2", first(res.inference_error_reports).msg)
    end

    # error handling for module usages
    let
        res, interp = @profile_toplevel begin
            using Base: foo
        end

        @test !isempty(res.toplevel_error_reports)
        er = first(res.toplevel_error_reports)
        @test er isa ActualErrorWrapped
        @test er.err == UndefVarError(:foo)
        @test er.file == (@__FILE__) && er.line == (@__LINE__) - 7
    end

    # sequential usage
    let
        res, interp = @profile_toplevel begin
            module foo

            bar(s) = sum(s)

            module baz

            using ..foo

            bar("julia") # -> GlobalUndefVarErrorReport

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        @test !isempty(res.inference_error_reports)
        @test first(res.inference_error_reports) isa GlobalUndefVarErrorReport
    end

    # usage of global objects
    let
        res, interp = @profile_toplevel begin
            module foo

            bar(s) = sum(s)

            module baz

            using ..foo

            foo.bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end

    let
        res, interp = @profile_toplevel begin
            module foo

            bar(s) = sum(s)

            module baz

            using ..foo: bar

            bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end

    # module usage within a block
    let
        res, interp = @profile_toplevel begin
            module foo

            bar(s) = sum(s)

            module baz

            begin
                using ..foo: bar
                bar("julia") # -> NoMethodErrorReports
            end

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end

    @testset "module usage of virtual global variable" begin
        let
            res, interp = @profile_toplevel begin
                module foo

                const bar = sum

                module baz

                using ..foo: bar

                bar("julia")

                end # module bar

                end # module foo
            end

            @test isempty(res.toplevel_error_reports)
            test_sum_over_string(res)
        end

        let
            res, interp = @profile_toplevel begin
                module foo

                const bar = "julia"

                module baz

                using ..foo: bar

                sum(bar)

                end # module bar

                end # module foo
            end

            @test isempty(res.toplevel_error_reports)
            test_sum_over_string(res)
        end

        let
            res, interp = @profile_toplevel begin
                module foo

                const bar = "julia"

                export bar

                end # module foo

                using .foo
                sum(bar)
            end

            @test isempty(res.toplevel_error_reports)
            test_sum_over_string(res)
        end
    end

    # export
    let
        res, interp = @profile_toplevel begin
            module foo

            bar(s) = sum(s)

            export bar

            end

            using .foo

            bar("julia") # -> NoMethodErrorReports
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end
end

@testset "sequential" begin
    let
        res, interp = @profile_toplevel begin
            foo(1:1000)

            foo(a) = length(a)
        end
        @test length(res.inference_error_reports) === 1
        @test first(res.inference_error_reports) isa GlobalUndefVarErrorReport
    end

    let
        res, interp = @profile_toplevel begin
            foo(a) = length(a)
            foo("julia") # should not error

            foo(a) = sum(a)
            foo(1:1000)  # should not error too
        end
        @test isempty(res.inference_error_reports)
    end
end

@testset "virtual global variables" begin
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            var = rand(Bool)
            const constvar = rand(Bool)
        end

        @test is_abstract(vmod, :var)
        @test isa_abstract(vmod.var, Bool)
        @test is_abstract(vmod, :constvar)
        @test isa_abstract(vmod.constvar, Bool)
    end

    @testset "scope" begin
        # function
        # --------

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                begin
                    local localvar = rand(Bool)
                    global globalvar = rand(Bool)
                end
            end

            @test !isdefined(vmod, :localvar)
            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Bool)
        end

        # blocks
        # ------

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                begin
                    globalvar = rand(Bool)
                end
            end

            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Bool)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                begin
                    local localvar = rand(Bool)
                    globalvar = localvar
                end
            end

            @test !isdefined(vmod, :localvar)
            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Bool)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                begin
                    local localvar
                    localvar = rand(Bool) # this shouldn't be annotated as `global`
                    globalvar = localvar
                end
            end

            @test !isdefined(vmod, :localvar)
            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Bool)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                globalvar2 = begin
                    local localvar = rand(Bool)
                    globalvar1 = localvar
                end
            end

            @test !isdefined(vmod, :locbalvar)
            @test is_abstract(vmod, :globalvar1)
            @test isa_abstract(vmod.globalvar1, Bool)
            @test is_abstract(vmod, :globalvar2)
            @test isa_abstract(vmod.globalvar2, Bool)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                let
                    localvar = rand(Bool)
                end
            end

            @test !isdefined(vmod, :localvar)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                globalvar = let
                    localvar = rand(Bool)
                    localvar
                end
            end

            @test !isdefined(vmod, :localvar)
            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Bool)
        end

        # loops
        # -----

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                for i = 1:100
                    localvar = i
                end
            end

            @test !isdefined(vmod, :localvar)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                for i in 1:10
                    localvar = rand(Bool)
                    global globalvar = localvar
                end
            end

            @test !isdefined(vmod, :localvar)
            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Bool)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                while true
                    localvar = rand(Bool)
                end
            end

            @test !isdefined(vmod, :localvar)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                while true
                    localvar = rand(Bool)
                    global globalvar = localvar
                end
            end

            @test !isdefined(vmod, :localvar)
            @test is_abstract(vmod, :globalvar)
            @test isa_abstract(vmod.globalvar, Bool)
        end
    end

    @testset "multiple declaration/assignment" begin
        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                s, c = sincos(1)
            end

            @test is_abstract(vmod, :s)
            @test isa_abstract(vmod.s, Float64)
            @test is_abstract(vmod, :c)
            @test isa_abstract(vmod.c, Float64)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                begin
                    local s, c
                    s, c = sincos(1)
                end
            end

            @test !isdefined(vmod, :s)
            @test !isdefined(vmod, :c)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                let
                    global s, c
                    s, c = sincos(1)
                end
            end

            @test is_abstract(vmod, :s)
            @test isa_abstract(vmod.s, Float64)
            @test is_abstract(vmod, :c)
            @test isa_abstract(vmod.c, Float64)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                so, co = let
                    si, ci = sincos(1)
                    si, ci
                end
            end

            @test !isdefined(vmod, :si)
            @test !isdefined(vmod, :ci)
            @test is_abstract(vmod, :so)
            @test isa_abstract(vmod.so, Float64)
            @test is_abstract(vmod, :co)
            @test isa_abstract(vmod.co, Float64)
        end

        let
            vmod = gen_virtual_module()
            res, interp = @profile_toplevel vmod begin
                begin
                    local l
                    l, g = sincos(1)
                end
            end
            @test !isdefined(vmod, :l)
            @test is_abstract(vmod, :g)
            @test isa_abstract(vmod.g, Float64)
        end
    end
end

@testset "toplevel throw" begin
    res, interp = @profile_toplevel begin
        throw("throw me")
    end
    @test length(res.inference_error_reports) == 1
    @test first(res.inference_error_reports) isa UncaughtExceptionReport
end

@testset "error handling within ConcreteInterpreter" begin
    let
        res, interp = @profile_toplevel begin
            struct A <: B end # UndefVarError(:B) should be handled into `res.toplevel_error_reports`
        end

        @test !isempty(res.toplevel_error_reports)
        er = first(res.toplevel_error_reports)
        @test er isa ActualErrorWrapped
        @test er.err == UndefVarError(:B)
        @test er.file == (@__FILE__) && er.line == (@__LINE__) - 7
    end

    # stacktrace cropping
    let
        res, interp = @profile_toplevel begin
            foo() = throw("don't call me, pal")
            struct A <: foo() end
        end

        @test !isempty(res.toplevel_error_reports)
        er = first(res.toplevel_error_reports)
        @test er isa ActualErrorWrapped
        @test er.err == "don't call me, pal"
        @test er.file == (@__FILE__) && er.line == (@__LINE__) - 7
        @test length(er.st) == 1
        sf = first(er.st)
        @test sf.file === Symbol(@__FILE__) && sf.line == (@__LINE__) - 11
    end
end

@testset "invalid constant redefinition" begin
    # for virtual global assignment
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            fib(n) = n≤2 ? n : fib(n-1)+fib(n-1)

            const foo = fib(1000000000000) # ::Int

            foo = fib(1000000000000.) # ::Float64
        end

        @test is_abstract(vmod, :foo)
        @test length(res.inference_error_reports) == 1
        er = first(res.inference_error_reports)
        @test er isa InvalidConstantRedefinition
        @test er.name === :foo
    end

    # for concretized constants
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            fib(n) = n≤2 ? n : fib(n-1)+fib(n-1)
            const T = typeof(fib(1000000000000)) # never terminates, yes

            T = Nothing
        end

        @test is_concrete(vmod, :T) # wao, this is concretized
        @test length(res.inference_error_reports) == 1
        er = first(res.inference_error_reports)
        @test er isa InvalidConstantRedefinition
        @test er.name === :T
    end
end

@testset "control flows in toplevel frames" begin
    let
        vmod = gen_virtual_module()
        res, interp = @profile_toplevel vmod begin
            v = rand(Int)

            if rand(Bool)
                v = rand(Char)
            end

            sin(v) # union-split no method error should be reported
        end

        @test is_abstract(vmod, :v)
        @test length(res.inference_error_reports) === 1
        er = first(res.inference_error_reports)
        @test er isa NoMethodErrorReport
        @test er.unionsplit
    end
end


text:
using Documenter, Literate, Individual

# Set Literate.jl config if not being compiled on recognized service.
config = Dict{String,String}()
if !(haskey(ENV, "GITHUB_ACTIONS") || haskey(ENV, "GITLAB_CI"))
  config["nbviewer_root_url"] = "https://nbviewer.jupyter.org/github/AlgebraicJulia/AlgebraicPetri.jl/blob/gh-pages/dev"
  config["repo_root_url"] = "https://github.com/AlgebraicJulia/AlgebraicPetri.jl/blob/master/docs"
end

const literate_dir = joinpath(@__DIR__, "..", "examples")
const generated_dir = joinpath(@__DIR__, "src", "examples")

for (root, dirs, files) in walkdir(literate_dir)
  out_dir = joinpath(generated_dir, relpath(root, literate_dir))
  for file in files
    f,l = splitext(file)
    if l == ".jl" && !startswith(f, "_")
      Literate.markdown(joinpath(root, file), out_dir;
        config=config, documenter=true, credit=false)
    end
  end
end

makedocs(
  sitename  = "Individual.jl",
  pages     = [
    "Individual.jl" => "index.md",
    "Examples" => Any[
        "examples/sir-basic.md",
        "examples/sir-scheduling.md",
        "examples/sir-age.md"
      ],
    "Library Reference" => "api.md",
    "Contributing" => "contributing.md"
  ]
)

deploydocs(
  repo = "github.com/slwu89/Individual.jl.git",
  branch = "gh-pages"
)


text:
const merge_algs = Dict(
    :molzahn => merge_molzahn,
    :sliwak => merge_sliwak
)

function merge(cliques::AbstractVector, cliquetree::AbstractVector, type::Symbol; kwargs...)
    return merge_algs[type](cliques, cliquetree; kwargs...)
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
import ACE.OrthPolys: OrthPolyBasis, discrete_jacobi

@doc raw"""
`Scal1pBasis`

One-particle basis of the form $P_n(x_i)$ for a general scalar, invariant
input `x`. This type basically just translates the `TransformedPolys` into a valid
one-particle basis.
"""
function Scal1pBasis(varsym::Symbol, varidx::Union{Integer, Nothing}, idxsym::Symbol,
                     P, label::String = "P$idxsym", trans=nothing)
   getval = isnothing(varidx) ? ACE.GetVal{varsym}() : ACE.GetVali{varsym, varidx}()
   spec = [ NamedTuple{(idxsym,)}((i,)) for i = 1:length(P) ]
   degrees = collect(0:(length(P)-1))
   P1 = (trans == nothing) ? P : chain(trans, P)
   return B1pComponent(P1, getval, spec, degrees, label)
end

function scal1pbasis(varsym::Symbol, idxsym::Symbol,
                     maxdeg::Integer, trans, rcut::Real, rin::Real=0.0;
                     varidx = nothing, label::String = "P$idxsym",
                     kwargs...)
   J = discrete_jacobi(maxdeg; trans=trans, xcut=rcut, xin=rin, kwargs...)
   return Scal1pBasis(varsym, varidx, idxsym, J, label, trans)
end


text:
using Plots

x=1:10
ys=rand(10,2)
scatter(x,ys)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module GMM

using PDMats
using ForwardDiff
using MathProgBase
using Ipopt
using StatsBase
using Reexport
using Distributions
@reexport using CovarianceMatrices
import MathProgBase: SolverInterface
import CovarianceMatrices.RobustVariance

export gmm, status, coef, objval, momentfunction, jacobian, mfvcov, vcov,
       optimal_W, TwoStepGMM, OneStepGMM, IterativeGMM

# ------------------ #
# Iteration managers #
# ------------------ #

abstract IterationManager

immutable OneStepGMM <: IterationManager
    k::RobustVariance
end

immutable TwoStepGMM <: IterationManager
    k::RobustVariance
end

immutable IterativeGMM <: IterationManager
    k::RobustVariance
    tol::Float64
    maxiter::Int
end

# kwarg constructors with default values
OneStepGMM(;k::RobustVariance=HC0()) = OneStepGMM(k)

TwoStepGMM(;k::RobustVariance=HC0()) = TwoStepGMM(k)

function IterativeGMM(;k::RobustVariance=HC0(), tol::Float64=1e-12,
                       maxiter::Int=500)
    IterativeGMM(k, tol, maxiter)
end

type IterationState
    n::Int
    change::Float64
    prev::Array  # previous value
end

finished(::OneStepGMM, ist::IterationState) = ist.n >= 1
finished(::TwoStepGMM, ist::IterationState) = ist.n >= 2
function finished(mgr::IterativeGMM, ist::IterationState)
    ist.n > mgr.maxiter || abs(ist.change) <= mgr.tol
end

# -------------- #
# Main GMM Types #
# -------------- #

type GMMNLPE <: SolverInterface.AbstractNLPEvaluator
    mf::Function
    smf::Function
    Dmf::Function
    mgr::IterationManager
    W::Array{Float64, 2}
end

abstract MomentEstimatorResult

type GMMResult <: MomentEstimatorResult
    status::Symbol
    objval::Real
    coef::Array{Float64, 1}
    nmom::Integer
    npar::Integer
    nobs::Integer
end

abstract MomentEstimator

type GMMEstimator <: MomentEstimator
    e::GMMNLPE
    r::GMMResult
end

# --------------- #
# Display methods #
# --------------- #

function show_extra(me::GMMEstimator)
    j, p = J_test(me, me.e.mgr.k)
    "\nJ-test: $(round(j, 3)) (P-value: $(round(p, 3)))\n"
end

# default to nothing
show_extra(me::MomentEstimator) = ""

function Base.writemime{T<:MomentEstimator}(io::IO, ::MIME"text/plain", me::T)
    # get coef table and j-test
    ct = coeftable(me, me.e.mgr.k)

    # show info for our model
    println(io, "$(T): $(npar(me)) parameter(s) with $(nmom(me)) moment(s)")

    # Show extra information for this type
    println(io, show_extra(me))

    # print coefficient table
    println(io, "Coefficients:\n")

    # then show coeftable
    show(io, ct)
end

include("util.jl")
include("gmm_mathprogbase.jl")
include("api.jl")
include("post_estimation.jl")

end


text:
"""
An attempt to understand Hirose-Ono paper entitled:

  Direct minimization to generate electronic states with proper
  occupation numbers
  Phys. Rev. B. 64 085105
"""

function test_main()
    #srand(1234)
    Nstates = 5
    Npoints = 15

    ϕ = rand(Npoints,Nstates) + im*rand(Npoints,Nstates)

    println("")
    for ist = 2:Nstates
        d = dot( ϕ[:,ist], ϕ[:,1] )
        @printf("States: %4d d = [%18.10f,%18.10f]\n", ist, real(d), imag(d))
    end

    U = ϕ'*ϕ

    # Calculate unitary matrix U which orthogonalizes ϕ
    Um = inv(sqrtm(U))  # Um = U^{-1/2}
    #println(abs.(eigvals(Um)))

    ϕ_ortho = ϕ * Um

    # Check orthogonalization wrt 1st state
    println("")
    println("Check orthogonality w.r.t 1st state")
    for ist = 2:Nstates
        d = dot( ϕ_ortho[:,ist], ϕ_ortho[:,1] )
        @printf("States: %4d d = [%18.10f,%18.10f]\n", ist, real(d), imag(d))
    end

    #c = copy(Um)
    c = sqrtm(U)   # <----- this is the c_{l,i}

    phi = zeros(Complex128,Npoints,Nstates)
    for i = 1:Nstates
        for l = 1:Nstates
            phi[:,i] = phi[:,i] + c[l,i] * ϕ_ortho[:,l]
        end
    end

    # Check orthogonalization wrt 1st states
    println("")
    for ist = 2:Nstates
        d = dot( phi[:,ist], phi[:,1] )
        @printf("States: %4d d = [%18.10f,%18.10f]\n", ist, real(d), imag(d))
    end

    T = zeros(Complex128,Nstates,Nstates)
    for k = 1:Nstates
        for l = 1:Nstates
            for i = 1:Nstates
                T[k,l] = T[k,l] + c[k,i]*c[l,i]
            end
        end
    end

    λ = eigvals(T)
    Focc = real((2 - λ).*λ)  # the occupation number, should not exceed 1.0

    println("\nOccupation numbers")
    for ist = 1:Nstates
        @printf("%4d %18.10f\n", ist, Focc[ist])
    end
    println("sum(Focc) = ", sum(Focc))

end

test_main()


text:
cd(@__DIR__)
using Pkg
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
Pkg.activate(@__DIR__)
CI && Pkg.instantiate()
using ClimateBase
using Documenter, DimensionalData, NCDatasets, StatsBase

# %%
using DocumenterTools: Themes
# download the themes
for file in ("juliaclimate-lightdefs.scss", "juliaclimate-darkdefs.scss", "juliaclimate-style.scss")
    download("https://raw.githubusercontent.com/JuliaClimate/doctheme/master/$file", joinpath(@__DIR__, file))
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliaclimate-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliaclimate-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliaclimate-$(w).scss"), header*"\n"*theme)
end
# compile the themes
Themes.compile(joinpath(@__DIR__, "juliaclimate-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "juliaclimate-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

# %% Build the actual docs
makedocs(modules = [ClimateBase, DimensionalData],
    sitename= "ClimateBase",
    authors = "George Datseris, Phillippe Roy, and contributors.",
    doctest = false,
    format = Documenter.HTML(
        prettyurls = CI,
        assets = [
            # "assets/logo.ico",
            asset("https://fonts.googleapis.com/css?family=Quicksand|Montserrat|Source+Code+Pro|Lora&display=swap", class=:css),
            ],
        ),
    pages = [
        "Introduction" => "index.md",
        ],
)

if CI
    deploydocs(
        repo = "github.com/JuliaClimate/ClimateBase.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
function GroundFieldOfViewDetector(arg0::Frame, arg1::FieldOfView)
    return GroundFieldOfViewDetector((Frame, FieldOfView), arg0, arg1)
end

function g(obj::GroundFieldOfViewDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_fov(obj::GroundFieldOfViewDetector)
    return jcall(obj, "getFOV", FieldOfView, ())
end

function get_frame(obj::GroundFieldOfViewDetector)
    return jcall(obj, "getFrame", Frame, ())
end


text:
using Base.Test
using AGRN

function get_diff(parent::GRN, child::GRN)
    diffs = 0
    for i in 1:100
        inputs = rand(parent.nin)
        parent_outs = step!(parent, inputs)
        child_outs = step!(parent, inputs)
        if any(parent_outs .!= child_outs)
            diffs += 1
        end
    end
    diffs
end

@testset "Crossover tests" begin
    config = AGRN.Config()
    nin = rand(1:10)
    nout = rand(1:10)
    nreg = rand(1:10)
    parent1 = GRN(nin, nout, rand(1:10), config)
    parent2 = GRN(nin, nout, rand(1:10), config)
    child = crossover(parent1, parent2, config)
    @test (child.nin == parent1.nin) && (child.nin == parent2.nin)
    @test (child.nout == parent1.nout) && (child.nout == parent2.nout)
    @test (child.nreg <= max(parent1.nreg, parent2.nreg))
    @test (child.beta == parent1.beta) || (child.beta == parent2.beta)
    @test (child.delta == parent1.delta) || (child.delta == parent2.delta)
    @test get_diff(parent1, child) > 0
    @test get_diff(parent2, child) > 0
    @test distance(parent1, child, config) <= distance(parent1, parent2, config)
    @test distance(parent2, child, config) <= distance(parent1, parent2, config)
end


text:
using GMC_NS, Distributions, ConjugatePriors

sample_dist=Normal(100.,5.)
samples=rand(sample_dist, 10000)

prior=NormalGamma(300.,2e-3,1.,10.) #a lousy prior

box=[0. 1000.
     1/1000. 1/eps()]

gmc=GMC_DEFAULTS
gmc[1]=5
gmc[2]=eps()

e=Normal_Ensemble("NE_test", 10, samples, prior, box, gmc...)

uds=Vector{Vector{Function}}([[convergence_display],[evidence_display],[info_display]])

lds=Vector{Vector{Function}}([[model_obs_display]])

converge_ensemble!(e,backup=(true,100),upper_displays=uds, lower_displays=lds, disp_rot_its=1000)


text:
m = AnSchorfheide()
system = compute_system(m)
Φ, Ψ, F_ϵ, F_u = DSGE.compute_system_function(system)
zero_sys = DSGE.zero_system_constants(system)

@testset "Check System, Transition, Measurement, and PseduoMeasurement access functions" begin
    @test typeof(system) == System{Float64}
    @test typeof(system[:transition]) == Transition{Float64}
    @test typeof(system[:measurement]) == Measurement{Float64}
    @test typeof(system[:pseudo_measurement]) == PseudoMeasurement{Float64}
    @test typeof(system[:transition][:TTT]) == Matrix{Float64}
    @test typeof(system[:transition][:RRR]) == Matrix{Float64}
    @test typeof(system[:transition][:CCC]) == Vector{Float64}
    @test typeof(system[:measurement][:ZZ]) == Matrix{Float64}
    @test typeof(system[:measurement][:DD]) == Vector{Float64}
    @test typeof(system[:measurement][:QQ]) == Matrix{Float64}
    @test typeof(system[:measurement][:EE]) == Matrix{Float64}
    @test typeof(system[:pseudo_measurement][:ZZ_pseudo]) == Matrix{Float64}
    @test typeof(system[:pseudo_measurement][:DD_pseudo]) == Vector{Float64}
end

@testset "Check miscellaneous functions acting on System types" begin
    @test sum(zero_sys[:CCC]) == 0.
    @test sum(zero_sys[:DD]) == 0.
    @test sum(zero_sys[:DD_pseudo]) == 0.
    @test Φ(ones(size(system[:transition][:TTT],2)), zeros(size(system[:transition][:RRR],2))) ==
        system[:transition][:TTT] * ones(size(system[:transition][:TTT],2))
    @test Ψ(ones(size(system[:transition][:TTT],2))) ==
        system[:measurement][:ZZ] * ones(size(system[:transition][:TTT],2)) + system[:measurement][:DD]
    @test F_ϵ.μ == zeros(3)
    @test F_ϵ.Σ.mat == system[:measurement][:QQ]
    @test F_u.μ == zeros(3)
    @test F_u.Σ.mat == system[:measurement][:EE]
end

@testset "Using compute_system to update an existing system" begin

    # Check updating
    system1 = compute_system(m, system)
    system2 = compute_system(m, system; observables = collect(keys(m.observables))[1:end - 1],
                                  pseudo_observables = collect(keys(m.pseudo_observables))[1:end - 1],
                                  shocks = collect(keys(m.exogenous_shocks))[1:end - 1],
                                  states = collect(keys(m.endogenous_states))[1:end - 1])
    system3 = compute_system(m, system; zero_DD = true)
    system4 = compute_system(m, system; zero_DD_pseudo = true)

    @test system1[:TTT] ≈ system[:TTT]
    @test system1[:RRR] ≈ system[:RRR]
    @test system1[:CCC] ≈ system[:CCC]
    @test system1[:ZZ] ≈ system[:ZZ]
    @test system1[:DD] ≈ system[:DD]
    @test system1[:QQ] ≈ system[:QQ]
    @test system1[:EE] ≈ zeros(size(system[:ZZ], 1), size(system[:ZZ], 1))
    @test system1[:ZZ_pseudo] ≈ system[:ZZ_pseudo]
    @test system1[:DD_pseudo] ≈ system[:DD_pseudo]

    @test system2[:TTT] ≈ system[:TTT][1:end - 1, 1:end - 1]
    @test system2[:RRR] ≈ system[:RRR][1:end - 1, 1:end - 1]
    @test sum(abs.(system2[:CCC])) ≈ sum(abs.(system[:CCC][1:end - 1])) ≈ 0.
    @test system2[:ZZ] ≈ system[:ZZ][1:end - 1, 1:end - 1]
    @test system2[:DD] ≈ system[:DD][1:end - 1]
    @test system2[:QQ] ≈ system[:QQ][1:end - 1, 1:end - 1]
    @test system2[:EE] ≈ zeros(size(system[:ZZ], 1) - 1, size(system[:ZZ], 1) - 1)
    @test system2[:ZZ_pseudo] ≈ system[:ZZ_pseudo][1:end - 1, 1:end - 1]
    @test system2[:DD_pseudo] ≈ system[:DD_pseudo][1:end - 1]

    @test system3[:TTT] ≈ system[:TTT]
    @test system3[:RRR] ≈ system[:RRR]
    @test system3[:CCC] ≈ system[:CCC]
    @test system3[:ZZ] ≈ system[:ZZ]
    @test system3[:DD] ≈ zeros(length(system[:DD]))
    @test system3[:QQ] ≈ system[:QQ]
    @test system3[:EE] ≈ zeros(size(system[:ZZ], 1), size(system[:ZZ], 1))
    @test system3[:ZZ_pseudo] ≈ system[:ZZ_pseudo]
    @test system3[:DD_pseudo] ≈ system[:DD_pseudo]

    @test system4[:TTT] ≈ system[:TTT]
    @test system4[:RRR] ≈ system[:RRR]
    @test system4[:CCC] ≈ system[:CCC]
    @test system4[:ZZ] ≈ system[:ZZ]
    @test system4[:DD] ≈ system[:DD]
    @test system4[:QQ] ≈ system[:QQ]
    @test system4[:EE] ≈ zeros(size(system[:ZZ], 1), size(system[:ZZ], 1))
    @test system4[:ZZ_pseudo] ≈ system[:ZZ_pseudo]
    @test system4[:DD_pseudo] ≈ zeros(length(system[:DD_pseudo]))

    # Check errors
    @test_throws ErrorException compute_system(m, system; observables = [:blah])
    @test_throws ErrorException compute_system(m, system; states = [:blah])
    @test_throws KeyError compute_system(m, system; shocks = [:blah])
    @test_throws ErrorException compute_system(m, system; pseudo_observables = [:blah])
end

@testset "VAR approximation of state space" begin
    m = Model1002("ss10"; custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true)))

    system = compute_system(m)
    system = compute_system(m, system; observables = [:obs_hours, :obs_gdpdeflator,
                                                      :laborshare_t, :NominalWageGrowth],
                            shocks = collect(keys(m.exogenous_shocks)))
    yyyyd, xxyyd, xxxxd = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                                 system[:DD], system[:ZZ], system[:EE],
                                                 zeros(size(system[:ZZ], 1),
                                                       DSGE.n_shocks_exogenous(m)),
                                                 4; get_population_moments = true)
    yyyydc, xxyydc, xxxxdc = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                                         system[:DD], system[:ZZ], system[:EE],
                                                         zeros(size(system[:ZZ], 1),
                                                               DSGE.n_shocks_exogenous(m)),
                                                         4; get_population_moments = true,
                                                         use_intercept = true)
    β, Σ = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ], system[:DD],
                                  system[:ZZ], system[:EE],
                                  zeros(size(system[:ZZ], 1), DSGE.n_shocks_exogenous(m)),
                                  4; get_population_moments = false)
    βc, Σc = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ], system[:DD],
                                  system[:ZZ], system[:EE],
                                  zeros(size(system[:ZZ], 1), DSGE.n_shocks_exogenous(m)),
                                  4; get_population_moments = false, use_intercept = true)

    expmat = load("reference/exp_var_approx_state_space.jld2")
    @test @test_matrix_approx_eq yyyyd expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyyd expmat["xxyyd"][2:end, :]
    @test @test_matrix_approx_eq xxxxd expmat["xxxxd"][2:end, 2:end]
    @test @test_matrix_approx_eq yyyydc expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyydc expmat["xxyyd"]
    @test @test_matrix_approx_eq xxxxdc expmat["xxxxd"]

    expβ = \(expmat["xxxxd"][2:end, 2:end], expmat["xxyyd"][2:end, :])
    expΣ = expmat["yyyyd"] - expmat["xxyyd"][2:end, :]' * expβ
    expΣ += expΣ'
    expΣ ./= 2
    @test @test_matrix_approx_eq β expβ
    @test @test_matrix_approx_eq Σ expΣ

    expβc = \(expmat["xxxxd"], expmat["xxyyd"])
    expΣc = expmat["yyyyd"] - expmat["xxyyd"]' * expβc
    expΣc += expΣc'
    expΣc ./= 2
    @test @test_matrix_approx_eq βc expβc
    @test @test_matrix_approx_eq Σc expΣc

    # Check DSGEVAR automates this properly
    dsgevar = DSGEVAR(m)
    DSGE.update!(dsgevar, shocks = collect(keys(m.exogenous_shocks)),
                 observables = [:obs_hours, :obs_gdpdeflator, :laborshare_t, :NominalWageGrowth],
                 lags = 4, λ = Inf)
    yyyyd, xxyyd, xxxxd = compute_system(dsgevar; get_population_moments = true)
    yyyydc, xxyydc, xxxxdc = compute_system(dsgevar; get_population_moments = true, use_intercept = true)
    β, Σ = compute_system(dsgevar)
    βc, Σc = compute_system(dsgevar; use_intercept = true)

    @test @test_matrix_approx_eq yyyyd expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyyd expmat["xxyyd"][2:end, :]
    @test @test_matrix_approx_eq xxxxd expmat["xxxxd"][2:end, 2:end]
    @test @test_matrix_approx_eq yyyydc expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyydc expmat["xxyyd"]
    @test @test_matrix_approx_eq xxxxdc expmat["xxxxd"]

    expβ = \(expmat["xxxxd"][2:end, 2:end], expmat["xxyyd"][2:end, :])
    expΣ = expmat["yyyyd"] - expmat["xxyyd"][2:end, :]' * expβ
    expΣ += expΣ'
    expΣ ./= 2
    @test @test_matrix_approx_eq β expβ
    @test @test_matrix_approx_eq Σ expΣ

    expβc = \(expmat["xxxxd"], expmat["xxyyd"])
    expΣc = expmat["yyyyd"] - expmat["xxyyd"]' * expβc
    expΣc += expΣc'
    expΣc ./= 2
    @test @test_matrix_approx_eq βc expβc
    @test @test_matrix_approx_eq Σc expΣc
end

@testset "VAR using DSGE as a prior" begin
    m = Model1002("ss10"; custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true)))
    dsgevar = DSGEVAR(m)
    jlddata = load(joinpath(dirname(@__FILE__), "reference/test_dsgevar_lambda_irfs.jld2"))
    DSGE.update!(dsgevar, shocks = collect(keys(m.exogenous_shocks)),
                 observables = [:obs_hours, :obs_gdpdeflator, :laborshare_t, :NominalWageGrowth],
                 lags = 4, λ = Inf)
    data = jlddata["data"]
    yyyydc1, xxyydc1, xxxxdc1 = compute_system(dsgevar; get_population_moments = true, use_intercept = true)
    βc1, Σc1 = compute_system(dsgevar; use_intercept = true)

    yyyyd2, xxyyd2, xxxxd2 = compute_system(dsgevar, data; get_population_moments = true)
    β2, Σ2 = compute_system(dsgevar, data)

    # Check when λ = Inf
    @test @test_matrix_approx_eq yyyydc1 yyyyd2
    @test @test_matrix_approx_eq xxyydc1 xxyyd2
    @test @test_matrix_approx_eq xxxxdc1 xxxxd2
    @test @test_matrix_approx_eq βc1 β2
    @test @test_matrix_approx_eq Σc1 Σ2

    # Check when λ is finite
    DSGE.update!(dsgevar, λ = 1.)
Random.seed!(1793) # need to seed for this
    yyyyd, xxyyd, xxxxd = compute_system(dsgevar, data; get_population_moments = true)
    β, Σ = compute_system(dsgevar, data)

    @test @test_matrix_approx_eq jlddata["exp_data_beta"] β
    @test @test_matrix_approx_eq jlddata["exp_data_sigma"] Σ
end

@testset "VECM approximation of state space" begin
    matdata = load("reference/vecm_approx_state_space.jld2")
    nobs = Int(matdata["nvar"])
    p = Int(matdata["nlags"])
    coint = Int(matdata["coint"])
    TTT = matdata["TTT"]
    RRR = matdata["RRR"]
    ZZ = matdata["ZZ"]
    DD = vec(matdata["DD"])
    QQ = matdata["QQ"]
    EE = matdata["EE"]
    MM = matdata["MM"]
    yyyyd, xxyyd, xxxxd = DSGE.vecm_approx_state_space(TTT, RRR, QQ,
                                                       DD, ZZ, EE, MM, nobs,
                                                       p, coint; get_population_moments = true,
                                                       test_GA0 = matdata["GA0"])
    yyyydc, xxyydc, xxxxdc = DSGE.vecm_approx_state_space(TTT, RRR, QQ,
                                                          DD, ZZ, EE, MM, nobs,
                                                          p, coint; get_population_moments = true,
                                                          use_intercept = true, test_GA0 = matdata["GA0"])
    β, Σ = DSGE.vecm_approx_state_space(TTT, RRR, QQ, DD,
                                        ZZ, EE, MM, nobs, p, coint;
                                        get_population_moments = false,
                                        test_GA0 = matdata["GA0"])
    βc, Σc = DSGE.vecm_approx_state_space(TTT, RRR, QQ, DD,
                                          ZZ, EE, MM, nobs, p, coint;
                                          get_population_moments = false, use_intercept = true,
                                          test_GA0 = matdata["GA0"])

    no_int_inds = vcat(1:coint, coint + 2:size(matdata["xxyyd"], 1))
    @test @test_matrix_approx_eq yyyyd matdata["yyyyd"]
    @test @test_matrix_approx_eq xxyyd matdata["xxyyd"][no_int_inds, :]
    @test @test_matrix_approx_eq xxxxd matdata["xxxxd"][no_int_inds, no_int_inds]
    @test @test_matrix_approx_eq yyyydc matdata["yyyyd"]
    @test @test_matrix_approx_eq xxyydc matdata["xxyyd"]
    @test @test_matrix_approx_eq xxxxdc matdata["xxxxd"]

    expβ = \(matdata["xxxxd"][no_int_inds, no_int_inds], matdata["xxyyd"][no_int_inds, :])
    expΣ = matdata["yyyyd"] - matdata["xxyyd"][no_int_inds, :]' * expβ
    @test @test_matrix_approx_eq β expβ
    @test @test_matrix_approx_eq Σ expΣ

    expβc = matdata["xxxxd"] \ matdata["xxyyd"]
    expΣc = matdata["yyyyd"] - matdata["xxyyd"]' * expβc
    @test @test_matrix_approx_eq βc expβc
    @test @test_matrix_approx_eq Σc expΣc

    # Check DSGEVECM automates this properly
    # m = Model1002()
    # dsgevecm = DSGEVECM(m)
    # DSGE.update!(dsgevecm, shocks = collect(keys(m.exogenous_shocks)),
    #              observables = [:obs_hours, :obs_gdpdeflator, :laborshare_t, :NominalWageGrowth],
    #              lags = 4, λ = Inf)
    # yyyyd, xxyyd, xxxxd = compute_system(dsgevecm; get_population_moments = true)
    # yyyydc, xxyydc, xxxxdc = compute_system(dsgevecm; get_population_moments = true, use_intercept = true)
    # β, Σ = compute_system(dsgevecm)
    # βc, Σc = compute_system(dsgevecm; use_intercept = true)

    # @test @test_matrix_approx_eq yyyyd matdata["yyyyd"]
    # @test @test_matrix_approx_eq xxyyd matdata["xxyyd"][2:end, :]
    # @test @test_matrix_approx_eq xxxxd matdata["xxxxd"][2:end, 2:end]
    # @test @test_matrix_approx_eq yyyydc matdata["yyyyd"]
    # @test @test_matrix_approx_eq xxyydc matdata["xxyyd"]
    # @test @test_matrix_approx_eq xxxxdc matdata["xxxxd"]

    # expβ = \(matdata["xxxxd"][2:end, 2:end], matdata["xxyyd"][2:end, :])
    # expΣ = matdata["yyyyd"] - matdata["xxyyd"][2:end, :]' * expβ
    # expΣ += expΣ'
    # expΣ ./= 2
    # @test @test_matrix_approx_eq β expβ
    # @test @test_matrix_approx_eq Σ expΣ

    # expβc = \(matdata["xxxxd"], matdata["xxyyd"])
    # expΣc = matdata["yyyyd"] - matdata["xxyyd"]' * expβc
    # expΣc += expΣc'
    # expΣc ./= 2
    # @test @test_matrix_approx_eq βc expβc
    # @test @test_matrix_approx_eq Σc expΣc
end

@testset "Updating a system for a DSGEVECM" begin
    dsge = AnSchorfheide()
    m = DSGEVECM(dsge)
    sys = compute_system(dsge)
    Dout1 = DSGE.compute_DD_coint_add(m, sys, [:obs_gdp, :obs_cpi])
    @info "The following 2 warnings about an empty vector are expected"
    sys_dsgevecm, Dout2 = compute_system(m, sys; get_DD_coint_add = true,
                                         cointegrating_add = [:obs_gdp, :obs_cpi])
    Dout3 = DSGE.compute_DD_coint_add(m, sys, Vector{Symbol}(undef, 0))
    _, Dout4 = compute_system(m, sys; get_DD_coint_add = true)

    # Check computing DD_coint_add
    @test Dout1 == sys[:DD][1:2]
    @test Dout2 == sys[:DD][1:2]
    @test isempty(Dout3)
    @test isempty(Dout4)

    # Check the system looks right
    @test @test_matrix_approx_eq sys[:TTT] sys_dsgevecm[:TTT]
    @test @test_matrix_approx_eq sys[:RRR] sys_dsgevecm[:RRR]
    @test @test_matrix_approx_eq sys[:CCC] sys_dsgevecm[:CCC]
    @test @test_matrix_approx_eq sys[:ZZ] sys_dsgevecm[:ZZ]
    @test @test_matrix_approx_eq sys[:DD] sys_dsgevecm[:DD]
    @test @test_matrix_approx_eq sys[:QQ] sys_dsgevecm[:QQ]
    @test @test_matrix_approx_eq sys_dsgevecm[:EE] zeros(size(sys_dsgevecm[:EE]))
end

nothing


text:
# Working around my inability to come up with algorithms
# A collection of functions returning sequences of flops for
# generating all binary functios at the given Hamming distance

# USAGE:
#
#   hamseq[depth](length)
#
# where depth is the Hamming distance and length is the number
# of bits in the string
# OUTPUT:
#   [depth, C] array of flip indices, where C is the number of strings at the given HD

function hamseq1(NB)
    out = [i for i in 1:NB]
    return Int.(out)
end

function hamseq2(NB)
	NF = 2
	conf = zeros(NF,1)
	out = zeros(NF, 1)

	for conf[1] in 1:(NB - 1)
		Nk = 0
		for conf[2] in (conf[1]+1):(NB - Nk)
			out = [out conf]
		end
	end
    return Int.(out[:,2:end])
end

function hamseq3(NB)
	NF = 3
	conf = zeros(NF,1)
	out = zeros(NF, 1)

	for conf[1] in 1:(NB - NF + 1)
		Nk = NF - 2

		for conf[2] in (conf[1]+1):(NB - Nk)
			Nk = NF - 3

			for conf[3] in (conf[2]+1):(NB - Nk)
				out = [out conf]
			end

		end
	end
    return Int.(out[:,2:end])
end

function hamseq4(NB)
	NF = 4
	conf = zeros(NF,1)
	out = zeros(NF, 1)

	for conf[1] in 1:(NB - NF + 1)
		Nk = NF - 2

		for conf[2] in (conf[1]+1):(NB - Nk)
			Nk = NF - 3

			for conf[3] in (conf[2]+1):(NB - Nk)
				Nk = NF - 4

				for conf[4] in (conf[3]+1):(NB - Nk)
					out = [out conf]
				end
			end
		end
	end
    return Int.(out[:,2:end])
end

function hamseq5(NB)
	NF = 5
	conf = zeros(NF,1)
	out = zeros(NF, 1)

	for conf[1] in 1:(NB - NF + 1)
		Nk = NF - 2

		for conf[2] in (conf[1]+1):(NB - Nk)
			Nk = NF - 3

			for conf[3] in (conf[2]+1):(NB - Nk)
				Nk = NF - 4

				for conf[4] in (conf[3]+1):(NB - Nk)
					Nk = NF - 5

					for conf[5] in (conf[4]+1):(NB - Nk)
						out = [out conf]

    				end
    			end
    		end
		end
	end
    return Int.(out[:,2:end])
end

hamseq = [hamseq1, hamseq2, hamseq3, hamseq4, hamseq5]

max_stars_repo_name: merement/Dice

max_issues_repo_name: merement/Dice

max_forks_repo_name: merement/Dice


text:
using Match

const PulseTable = map(x -> 95.52f0 / (8128.0f0/Float32(x) + 100.0f0), 0:30)

const TNDTable = map(x -> 163.67f0 / (24329.0f0/Float32(x) + 100.0f0), 0:202)

function apustep!(console::Console, apu::APU)
  apu.cycle += 1
  steptimer!(console, apu)
  #f1 = floor(Int32, Float64(cycle1) / frameCounterRate)
  #f2 = floor(Int32, Float64(cycle2) / frameCounterRate)
  apu.frameCounter += 240
  if apu.frameCounter >= cpuFrequency
    stepframecounter!(console)
    apu.frameCounter -= cpuFrequency
  end
  apu.sampleCounter += 1.0
  if apu.sampleCounter >= apu.sampleRate
    sendsample!(apu)
    apu.sampleCounter -= apu.sampleRate
  end
  #s1 = floor(Int32, Float64(cycle1) / apu.sampleRate)
  #s2 = floor(Int32, Float64(cycle2) / apu.sampleRate)
  #if s1 != s2
    #sendsample!(apu)
  #end
end

# mode 0:    mode 1:       function
# ---------  -----------  -----------------------------
#  - - - f    - - - - -    IRQ (if bit 6 is clear)
#  - l - l    l - l - -    Length counter and sweep
#  e e e e    e e e e -    Envelope and linear counter
function stepframecounter!(console::Console)
  apu = console.apu
  @match apu.framePeriod begin
    0x04 => begin
      apu.frameValue = (apu.frameValue + 0x01) & 0x03
      @match apu.frameValue begin
        0x00 => begin
          stepenvelope!(apu)
        end
        0x01 => begin
          stepenvelope!(apu)
          stepsweep!(apu)
          steplength!(apu)
        end
        0x02 => begin
          stepenvelope!(apu)
        end
        0x03 => begin
          stepenvelope!(apu)
          stepsweep!(apu)
          steplength!(apu)
          fireirq!(console)
        end
      end
    end
    0x05 => begin
      apu.frameValue = (apu.frameValue + 0x01) % 0x05
      @match apu.frameValue begin
        0x00 => begin
          stepenvelope!(apu)
        end
        0x01 => begin
          stepenvelope!(apu)
          stepsweep!(apu)
          steplength!(apu)
        end
        0x02 => begin
          stepenvelope!(apu)
        end
        0x04 => begin
          stepenvelope!(apu)
          stepsweep!(apu)
          steplength!(apu)
        end
      end
    end
  end
end

function sendsample!(apu::APU)
  out = step!(apu.filterChain, output!(apu))
  write!(apu.channel, out)
end

function output!(apu::APU)::Float32
  p1 = output!(apu.pulse1)
  p2 = output!(apu.pulse2)
  t = output!(apu.triangle)
  n = output!(apu.noise)
  d = output!(apu.dmc)
  pulseOut = PulseTable[p1 + p2 + 1]
  tndOut = TNDTable[3 * t + 2 * n + d + 1]
  return pulseOut + tndOut
end

function stepenvelope!(apu::APU)
  stepenvelope!(apu.pulse1)
  stepenvelope!(apu.pulse2)
  stepcounter!(apu.triangle)
  stepenvelope!(apu.noise)
end

function stepsweep!(apu::APU)
  stepsweep!(apu.pulse1)
  stepsweep!(apu.pulse2)
end

function steplength!(apu::APU)
  steplength!(apu.pulse1)
  steplength!(apu.pulse2)
  steplength!(apu.triangle)
  steplength!(apu.noise)
end

function readregister(apu::APU, address::UInt16)::UInt8
  @match address begin
    0x4015 => readstatus(apu)
    _ => begin
      # default:
      # TODO: add logging.
      # log.Fatalf("unhandled apu register read at address: 0x%04X", address)
      0x00
    end
  end
end

function apuwriteregister!(console::Console, address::UInt16, val::UInt8)
  apu = console.apu
  @match address begin
    0x4000 => begin
      writecontrol!(apu.pulse1, val)
    end
    0x4001 => begin
      writesweep!(apu.pulse1, val)
    end
    0x4002 => begin
      writetimerlow!(apu.pulse1, val)
    end
    0x4003 => begin
      writetimerhigh!(apu.pulse1, val)
    end
    0x4004 => begin
      writecontrol!(apu.pulse2, val)
    end
    0x4005 => begin
      writesweep!(apu.pulse2, val)
    end
    0x4006 => begin
      writetimerlow!(apu.pulse2, val)
    end
    0x4007 => begin
      writetimerhigh!(apu.pulse2, val)
    end
    0x4008 => begin
      writecontrol!(apu.triangle, val)
    end
    0x4009 => begin
    end
    0x4010 => begin
      writecontrol!(apu.dmc, val)
    end
    0x4011 => begin
      writevalue!(apu.dmc, val)
    end
    0x4012 => begin
      writeaddress!(apu.dmc, val)
    end
    0x4013 => begin
      writelength!(apu.dmc, val)
    end
    0x400A => begin
      writetimerlow!(apu.triangle, val)
    end
    0x400B => begin
      writetimerhigh!(apu.triangle, val)
    end
    0x400C => begin
      writecontrol!(apu.noise, val)
    end
    0x400D => begin
    end
    0x400E => begin
      writeperiod!(apu.noise, val)
    end
    0x400F => begin
      writelength!(apu.noise, val)
    end
    0x4015 => begin
      writecontrol!(apu, val)
    end
    0x4017 => begin
      writeframecounter!(apu, val)
    end
    _ => begin
      # default:
      #   log.Fatalf("unhandled apu register write at address: 0x%04X", address)
    end
  end
end

function readstatus(apu::APU)::UInt8
  result = 0x00
  if apu.pulse1.lengthValue > 0
    result |= 0x01
  end
  if apu.pulse2.lengthValue > 0
    result |= 0x02
  end
  if apu.triangle.lengthValue > 0
    result |= 0x04
  end
  if apu.noise.lengthValue > 0
    result |= 0x08
  end
  if apu.dmc.currentLength > 0
    result |= 0x10
  end
  return result
end

function writecontrol!(apu::APU, val::UInt8)
  pulse1 = apu.pulse1
  pulse1.enabled = (val & 1) == 1
  if !pulse1.enabled
    pulse1.lengthValue = 0
  end
  pulse2 = apu.pulse2
  pulse2.enabled = (val & 2) == 2
  if !pulse2.enabled
    pulse2.lengthValue = 0
  end
  triangle = apu.triangle
  triangle.enabled = (val & 4) == 4
  if !triangle.enabled
    triangle.lengthValue = 0
  end
  noise = apu.noise
  noise.enabled = (val & 8) == 8
  if !noise.enabled
    noise.lengthValue = 0
  end
  dmc = apu.dmc
  dmc.enabled = (val & 16) == 16
  if !dmc.enabled
    dmc.currentLength = 0
  else
    if dmc.currentLength == 0
      restart!(dmc)
    end
  end
end

function writeframecounter!(apu::APU, val::UInt8)
  apu.framePeriod = 0x04 + ((val >> 7) & 0x01)
  apu.frameIRQ = ((val >> 6) & 0x01) == 0x00
  # frameValue = 0
  if apu.framePeriod == 0x05
    stepenvelope!(apu)
    stepsweep!(apu)
    steplength!(apu)
  end
end


text:
"""

## fsparm!

Function fsparm assembles fin_el matrices into a Julia sparse
global stiffness matrix.

### Method
```julia
fsparm!(gsm, g, km)
```

### Arguments
```julia
* gsm::SparseArrays{Float64, Float64}   : Sparse stiffnes matrix (Updated)
* g::Vector{Int}                      : Global coordinate vector.
* km::Matrix{Float64}                   : Stiffness matrix.
```
"""
function fsparm!(gsm, g, km)
  #println("g = $g")
  ndof = size(g, 1)
  #println("ndof = $ndof")
  for i in 1:ndof
    k = g[i]
    #println("g[$i] = k = $k")
    if k !== 0
      for j in 1:ndof
        #println("i = $i, j = $j, g[$j] = $(g[j])")
        if g[j] !== 0
          iw = k - g[j]
          #println("iw = $iw")
          if iw >= 0
            if g[i] == g[j]
              gsm[g[i], g[j]] += km[i, j]
            else
              gsm[g[i], g[j]] += km[i, j]
              gsm[g[j], g[i]] += km[i, j]
            end
            #println("gsm[$(g[i]), $(g[j])] += $(km[i,j])")
            #full(gsm) |> display
            #println()
          end
        end
      end
    end
  end
end


text:
module TestComposite

using Test
using Mimi

import Mimi:
    ComponentId, ComponentPath, ComponentDef, AbstractComponentDef,
    CompositeComponentDef, ModelDef, time_labels, compdef, find_comp,
    import_params!, CompositeVariableDef, CompositeParameterDef, components,
    dim_names

@defcomp Comp1 begin
    par_1_1 = Parameter(index=[time])      # external input
    var_1_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_1_1[t] = p.par_1_1[t]
    end
end

@defcomp Comp2 begin
    par_2_1 = Parameter(index=[time])      # connected to Comp1.var_1_1
    par_2_2 = Parameter(index=[time])      # external input
    var_2_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_2_1[t] = p.par_2_1[t] + p.foo * p.par_2_2[t]
    end
end

@defcomp Comp3 begin
    par_3_1 = Parameter(index=[time])      # connected to Comp2.var_2_1
    var_3_1 = Variable(index=[time])       # external output
    foo = Parameter(default=30)

    function run_timestep(p, v, d, t)
        # @info "Comp3 run_timestep"
        v.var_3_1[t] = p.par_3_1[t] * 2
    end
end

@defcomp Comp4 begin
    par_4_1 = Parameter(index=[time])      # connected to Comp2.var_2_1
    var_4_1 = Variable(index=[time])        # external output
    foo = Parameter(default=300)

    function run_timestep(p, v, d, t)
        # @info "Comp4 run_timestep"
        v.var_4_1[t] = p.par_4_1[t] * 2
    end
end

m = Model()
set_dimension!(m, :time, 2005:2020)

@defcomposite A begin
    Component(Comp1)
    Component(Comp2)

    foo1 = Parameter(Comp1.foo)
    foo2 = Parameter(Comp2.foo)

    var_2_1 = Variable(Comp2.var_2_1)

    connect(Comp2.par_2_1, Comp1.var_1_1)
    connect(Comp2.par_2_2, Comp1.var_1_1)
end

@defcomposite B begin
    Component(Comp3)
    Component(Comp4)

    foo3 = Parameter(Comp3.foo)
    foo4 = Parameter(Comp4.foo)

    var_3_1 = Variable(Comp3.var_3_1)
end

@defcomposite top begin
    Component(A)

    fooA1 = Parameter(A.foo1)
    fooA2 = Parameter(A.foo2)

    # TBD: component B isn't getting added to mi
    Component(B)
    foo3 = Parameter(B.foo3)
    foo4 = Parameter(B.foo4)

    var_3_1 = Variable(B.var_3_1)

    connect(B.par_3_1, A.var_2_1)
    connect(B.par_4_1, B.var_3_1)
end

# We have created the following composite structure:
#
#          top
#        /    \
#       A       B
#     /  \     /  \
#    1    2   3    4

top_ref = add_comp!(m, top, nameof(top))
top_comp = compdef(top_ref)

md = m.md

@test find_comp(md, :top) == top_comp

#
# Test various ways to access sub-components
#
c1 = find_comp(md, ComponentPath(:top, :A, :Comp1))
@test c1.comp_id == Comp1.comp_id

c2 = md[:top][:A][:Comp2]
@test c2.comp_id == Comp2.comp_id

c3 = find_comp(md, "/top/B/Comp3")
@test c3.comp_id == Comp3.comp_id

add_shared_param!(m, :model_fooA1, 1)
connect_param!(m, :top, :fooA1, :model_fooA1)

add_shared_param!(m, :model_fooA2, 2)
connect_param!(m, :top, :fooA2, :model_fooA2)

@test_throws ErrorException add_shared_param!(m, :model_par_1_1, collect(1:length(time_labels(md)))) # need to give index
add_shared_param!(m, :model_par_1_1, collect(1:length(time_labels(md))), dims = [:time])
connect_param!(m, :top, :par_1_1, :model_par_1_1)

Mimi.build!(m)

run(m)

mi = m.mi

# test parameters and variables fields of CompositeComponentInstance
top_var_keys = keys(mi[:top].variables)
top_par_keys = keys(mi[:top].parameters)
for item in md[:top].namespace
    if isa(item.second, CompositeVariableDef)
        @test in(item.first, top_var_keys)
    elseif isa(item.second, CompositeParameterDef)
        @test in(item.first, top_par_keys)
    end
end

# test access methods
@test mi[:top][:A][:Comp2, :par_2_2] == collect(1.0:16.0)
@test mi["/top/A/Comp2", :par_2_2] == collect(1.0:16.0)

@test mi["/top/A/Comp2", :var_2_1] == collect(3.0:3:48.0)
@test mi["/top/A/Comp1", :var_1_1] == collect(1.0:16.0)
@test mi["/top/B/Comp4", :par_4_1] == collect(6.0:6:96.0)

@test m[:top, :fooA1] == 1
@test m[:top, :foo3] == 30.
@test m[:top, :var_3_1] == collect(6.0:6:96.0)

# test ways to drill down into composites to get information
@test collect(keys(components(top))) == [:A, :B]
@test collect(keys(components(m, :top))) == [:A, :B]
@test collect(keys(components(Comp1))) == []
@test collect(keys(components(Comp1))) == []

@test dim_names(m, :Comp1, :par_1_1) == [:time]
@test dim_names(m, :Comp1, :var_1_1) == [:time]
@test dim_names(m, :A, :var_2_1) == [:time]

# Test joining external params.
#
m2 = Model()
set_dimension!(m2, :time, 2005:2020)

@defcomposite top2 begin
    Component(Comp1)
    Component(Comp2)

    connect(Comp2.par_2_1, Comp1.var_1_1)
    connect(Comp2.par_2_2, Comp1.var_1_1)

    foo = Parameter(Comp1.foo, Comp2.foo)
end

top2_ref = add_comp!(m2, top2, nameof(top2))

#
# Test _fix_comp_path on internal connections 3 levels down
#

@defcomposite top3 begin
    Component(top)
end

@test top3[:top][:A].internal_param_conns[1].src_comp_path == Mimi.ComponentPath(:top3, :top, :A, :Comp1)

path1 = ComponentPath(:a, :b)
path2 = ComponentPath(:c, :d)
@test ComponentPath(path1, path2) == ComponentPath(:a, :b, :c, :d)

# renaming

@defcomp A begin
    p1 = Parameter()
    p2 = Parameter()

    v1 = Variable()
end

@defcomp B begin
    v2 = Variable()
end

@defcomposite C begin

    foo = Component(A)
    bar = Component(B)

    rename_p1 = Parameter(foo.p1)

    connect(foo.p2, bar.v2)

    rename_v1 = Variable(foo.v1)
end

for key in [:foo, :bar, :rename_p1, :rename_v1]
    @test key in keys(C.namespace)
end

end # module

nothing