text:
# the current julia logo and graphics

# not exported

const julia_blue    = (0.251, 0.388, 0.847)
const julia_purple  = (0.584, 0.345, 0.698)
const julia_green   = (0.22, 0.596, 0.149)
const julia_red     = (0.796, 0.235, 0.2)

# to be deprecated
const darker_blue    = (0.251, 0.388, 0.847)
const lighter_blue   = (0.4, 0.51, 0.878)
const darker_purple  = (0.584, 0.345, 0.698)
const lighter_purple = (0.667, 0.475, 0.757)
const darker_green   = (0.22, 0.596, 0.149)
const lighter_green  = (0.376, 0.678, 0.318)
const darker_red     = (0.796, 0.235, 0.2)
const lighter_red    = (0.835, 0.388, 0.361)
const purples        = (darker_purple, lighter_purple)
const greens         = (darker_green, lighter_green)
const reds           = (darker_red, lighter_red)

"""
    julialogo(;
        action=:fill,
        color=true,
        bodycolor=colorant"black",
        centered=false)

Draw the Julia logo. The default action is to fill the logo and use the colors:

    julialogo()

If `color` is `false`, the `bodycolor` color is used for the logo.

The function uses the current drawing state (position, scale, etc).

The `centered` keyword lets you center the logo at its
mathematical center, but the optical center might lie
somewhere else - it's difficult to position well due to its
asymmetric design.

To use the logo as a clipping mask:

```
julialogo(action=:clip)
```

(In this case the `color` setting is automatically ignored.)

To obtain a stroked (outlined) version:

```
julialogo(action=:path)
sethue("red")
strokepath()
```
"""
function julialogo(;
        action    = :fill,
        color     = true,
        bodycolor = colorant"black",
        centered  = false)

    if centered == true
        translate(-Point(330, 213)/2)
    end

    # save current color
    r, g, b, a = Luxor.get_current_redvalue(), Luxor.get_current_greenvalue(), Luxor.get_current_bluevalue(), Luxor.get_current_alpha()

    # "j" without dot
    _j  =  [
    PathMove(Point(72.87109375, 177.3125)),
    PathCurve(Point(72.87109375, 184.84765625), Point(72.0234375, 190.93359375), Point(70.328125, 195.56640625)),
    PathCurve(Point(68.6328125, 200.203125), Point(66.22265625, 203.80078125), Point(63.09375, 206.36328125)),
    PathCurve(Point(59.96875, 208.92578125), Point(56.21875, 210.640625), Point(51.84765625, 211.5078125)),
    PathCurve(Point(47.4765625, 212.37109375), Point(42.61328125, 212.8046875), Point(37.265625, 212.8046875)),
    PathCurve(Point(30.02734375, 212.8046875), Point(24.48828125, 211.67578125), Point(20.6484375, 209.4140625)),
    PathCurve(Point(16.8046875, 207.15234375), Point(14.8828125, 204.44140625), Point(14.8828125, 201.2734375)),
    PathCurve(Point(14.8828125, 198.63671875), Point(15.953125, 196.4140625), Point(18.1015625, 194.60546875)),
    PathCurve(Point(20.25, 192.796875), Point(23.1328125, 191.89453125), Point(26.75, 191.89453125)),
    PathCurve(Point(29.46484375, 191.89453125), Point(31.6328125, 192.62890625), Point(33.25, 194.09765625)),
    PathCurve(Point(34.87109375, 195.56640625), Point(36.2109375, 197.01953125), Point(37.265625, 198.44921875)),
    PathCurve(Point(38.46875, 200.03125), Point(39.48828125, 201.0859375), Point(40.31640625, 201.61328125)),
    PathCurve(Point(41.14453125, 202.140625), Point(41.8984375, 202.40625), Point(42.578125, 202.40625)),
    PathCurve(Point(44.0078125, 202.40625), Point(45.1015625, 201.55859375), Point(45.85546875, 199.86328125)),
    PathCurve(Point(46.609375, 198.16796875), Point(46.984375, 194.87109375), Point(46.984375, 189.97265625)),
    PathLine(Point(46.984375, 97.05078125)),
    PathLine(Point(72.87109375, 89.9296875)),
    PathClose(),]

    # "u"
    _u = [
    PathMove(Point(109.73828125, 92.4140625)),
    PathLine(Point(109.73828125, 152.21484375)),
    PathCurve(Point(109.73828125, 153.875), Point(110.05859375, 155.4375), Point(110.69921875, 156.90625)),
    PathCurve(Point(111.33984375, 158.375), Point(112.2265625, 159.640625), Point(113.35546875, 160.6953125)),
    PathCurve(Point(114.48828125, 161.75), Point(115.8046875, 162.59765625), Point(117.3125, 163.23828125)),
    PathCurve(Point(118.8203125, 163.87890625), Point(120.44140625, 164.19921875), Point(122.17578125, 164.19921875)),
    PathCurve(Point(124.1328125, 164.19921875), Point(126.359375, 163.1015625), Point(129.0703125, 161.203125)),
    PathCurve(Point(133.36328125, 158.1953125), Point(135.96484375, 156.12890625), Point(135.96484375, 153.68359375)),
    PathCurve(Point(135.96484375, 153.09765625), Point(135.96484375, 92.4140625), Point(135.96484375, 92.4140625)),
    PathLine(Point(161.73828125, 92.4140625)),
    PathLine(Point(161.73828125, 177.3125)),
    PathLine(Point(135.96484375, 177.3125)),
    PathLine(Point(135.96484375, 169.3984375)),
    PathCurve(Point(132.57421875, 172.26171875), Point(128.95703125, 174.55859375), Point(125.11328125, 176.29296875)),
    PathCurve(Point(121.26953125, 178.02734375), Point(117.5390625, 178.89453125), Point(113.921875, 178.89453125)),
    PathCurve(Point(109.703125, 178.89453125), Point(105.78125, 178.1953125), Point(102.1640625, 176.80078125)),
    PathCurve(Point(98.546875, 175.40625), Point(95.3828125, 173.50390625), Point(92.671875, 171.09375)),
    PathCurve(Point(89.95703125, 168.68359375), Point(87.828125, 165.85546875), Point(86.28125, 162.61328125)),
    PathCurve(Point(84.73828125, 159.375), Point(83.96484375, 155.90625), Point(83.96484375, 152.21484375)),
    PathLine(Point(83.96484375, 92.4140625)),
    PathClose()]

    # "l"
    _l = [
    PathMove(Point(197.8828125, 177.3125)),
    PathLine(Point(172.22265625, 177.3125)),
    PathLine(Point(172.22265625, 58.27734375)),
    PathLine(Point(197.8828125, 51.15625)),
    PathClose()]

    #  "u"
    _i = [
    PathMove(Point(208.6015625, 97.05078125)),
    PathLine(Point(234.375, 89.9296875)),
    PathLine(Point(234.375, 177.3125)),
    PathLine(Point(208.6015625, 177.3125)),
    PathClose()]

    _a = [
    PathMove(Point(288.2265625, 133.44921875)),
    PathCurve(Point(285.73828125, 134.5078125), Point(283.23046875, 135.73046875), Point(280.70703125, 137.125)),
    PathCurve(Point(278.18359375, 138.51953125), Point(275.8828125, 140.046875), Point(273.8125, 141.703125)),
    PathCurve(Point(271.73828125, 143.359375), Point(270.0625, 145.1328125), Point(268.78125, 147.015625)),
    PathCurve(Point(267.5, 148.8984375), Point(266.859375, 150.859375), Point(266.859375, 152.89453125)),
    PathCurve(Point(266.859375, 154.4765625), Point(267.06640625, 156.00390625), Point(267.48046875, 157.47265625)),
    PathCurve(Point(267.89453125, 158.94140625), Point(268.48046875, 160.203125), Point(269.234375, 161.2578125)),
    PathCurve(Point(269.98828125, 162.3125), Point(270.81640625, 163.16015625), Point(271.72265625, 163.80078125)),
    PathCurve(Point(272.625, 164.44140625), Point(273.60546875, 164.76171875), Point(274.66015625, 164.76171875)),
    PathCurve(Point(276.76953125, 164.76171875), Point(278.8984375, 164.12109375), Point(281.046875, 162.83984375)),
    PathCurve(Point(283.1953125, 161.55859375), Point(285.5859375, 159.94140625), Point(288.2265625, 157.98046875)),
    PathClose(),
    PathMove(Point(314.109375, 177.3125)),
    PathLine(Point(288.2265625, 177.3125)),
    PathLine(Point(288.2265625, 170.52734375)),
    PathCurve(Point(286.79296875, 171.734375), Point(285.3984375, 172.84765625), Point(284.04296875, 173.86328125)),
    PathCurve(Point(282.6875, 174.87890625), Point(281.16015625, 175.765625), Point(279.46484375, 176.51953125)),
    PathCurve(Point(277.76953125, 177.2734375), Point(275.8671875, 177.85546875), Point(273.75390625, 178.2734375)),
    PathCurve(Point(271.64453125, 178.6875), Point(269.15625, 178.89453125), Point(266.296875, 178.89453125)),
    PathCurve(Point(262.375, 178.89453125), Point(258.8515625, 178.328125), Point(255.7265625, 177.19921875)),
    PathCurve(Point(252.59765625, 176.06640625), Point(249.94140625, 174.5234375), Point(247.7578125, 172.5625)),
    PathCurve(Point(245.5703125, 170.60546875), Point(243.89453125, 168.28515625), Point(242.7265625, 165.609375)),
    PathCurve(Point(241.55859375, 162.9375), Point(240.97265625, 160.015625), Point(240.97265625, 156.8515625)),
    PathCurve(Point(240.97265625, 153.609375), Point(241.59375, 150.671875), Point(242.83984375, 148.03125)),
    PathCurve(Point(244.08203125, 145.39453125), Point(245.77734375, 143.0234375), Point(247.92578125, 140.91015625)),
    PathCurve(Point(250.07421875, 138.80078125), Point(252.578125, 136.91796875), Point(255.44140625, 135.2578125)),
    PathCurve(Point(258.3046875, 133.6015625), Point(261.37890625, 132.07421875), Point(264.65625, 130.6796875)),
    PathCurve(Point(267.93359375, 129.28515625), Point(271.34375, 128.0078125), Point(274.88671875, 126.83984375)),
    PathCurve(Point(278.42578125, 125.671875), Point(281.93359375, 124.55859375), Point(285.3984375, 123.50390625)),
    PathLine(Point(288.2265625, 122.82421875)),
    PathLine(Point(288.2265625, 114.4609375)),
    PathCurve(Point(288.2265625, 109.03515625), Point(287.1875, 105.19140625), Point(285.1171875, 102.9296875)),
    PathCurve(Point(283.04296875, 100.66796875), Point(280.2734375, 99.5390625), Point(276.80859375, 99.5390625)),
    PathCurve(Point(272.73828125, 99.5390625), Point(269.91015625, 100.51953125), Point(268.328125, 102.4765625)),
    PathCurve(Point(266.74609375, 104.4375), Point(265.953125, 106.80859375), Point(265.953125, 109.59765625)),
    PathCurve(Point(265.953125, 111.1796875), Point(265.78515625, 112.7265625), Point(265.4453125, 114.234375)),
    PathCurve(Point(265.109375, 115.7421875), Point(264.5234375, 117.05859375), Point(263.6953125, 118.19140625)),
    PathCurve(Point(262.8671875, 119.3203125), Point(261.6796875, 120.2265625), Point(260.1328125, 120.90234375)),
    PathCurve(Point(258.58984375, 121.58203125), Point(256.6484375, 121.921875), Point(254.3125, 121.921875)),
    PathCurve(Point(250.6953125, 121.921875), Point(247.7578125, 120.8828125), Point(245.49609375, 118.8125)),
    PathCurve(Point(243.234375, 116.73828125), Point(242.10546875, 114.12109375), Point(242.10546875, 110.953125)),
    PathCurve(Point(242.10546875, 108.015625), Point(243.1015625, 105.28515625), Point(245.09765625, 102.76171875)),
    PathCurve(Point(247.09765625, 100.234375), Point(249.7890625, 98.06640625), Point(253.18359375, 96.26171875)),
    PathCurve(Point(256.57421875, 94.44921875), Point(260.4921875, 93.01953125), Point(264.9375, 91.96484375)),
    PathCurve(Point(269.3828125, 90.91015625), Point(274.09375, 90.3828125), Point(279.06640625, 90.3828125)),
    PathCurve(Point(285.171875, 90.3828125), Point(290.4296875, 90.9296875), Point(294.83984375, 92.01953125)),
    PathCurve(Point(299.24609375, 93.11328125), Point(302.8828125, 94.67578125), Point(305.74609375, 96.7109375)),
    PathCurve(Point(308.609375, 98.74609375), Point(310.71875, 101.1953125), Point(312.07421875, 104.05859375)),
    PathCurve(Point(313.43359375, 106.921875), Point(314.109375, 110.12890625), Point(314.109375, 113.66796875)),
    PathClose()]

    _red_dot = [
    PathMove(Point(240.2734375, 68.08984375)),
    PathCurve(Point(240.2734375, 77.7578125), Point(232.4375, 85.58984375), Point(222.7734375, 85.58984375)),
    PathCurve(Point(213.10546875, 85.58984375), Point(205.2734375, 77.7578125), Point(205.2734375, 68.08984375)),
    PathCurve(Point(205.2734375, 58.42578125), Point(213.10546875, 50.58984375), Point(222.7734375, 50.58984375)),
    PathCurve(Point(232.4375, 50.58984375), Point(240.2734375, 58.42578125), Point(240.2734375, 68.08984375)),
    ]

    _blue_dot = [
    PathMove(Point(77.953125, 68.08984375)),
    PathCurve(Point(77.953125, 77.7578125), Point(70.1171875, 85.58984375), Point(60.453125, 85.58984375)),
    PathCurve(Point(50.7890625, 85.58984375), Point(42.953125, 77.7578125), Point(42.953125, 68.08984375)),
    PathCurve(Point(42.953125, 58.42578125), Point(50.7890625, 50.58984375), Point(60.453125, 50.58984375)),
    PathCurve(Point(70.1171875, 50.58984375), Point(77.953125, 58.42578125), Point(77.953125, 68.08984375)),
    ]

    _purple_dot = [
    PathMove(Point(282.3203125, 68.08984375)),
    PathCurve(Point(282.3203125, 77.7578125), Point(274.484375, 85.58984375), Point(264.8203125, 85.58984375)),
    PathCurve(Point(255.15625, 85.58984375), Point(247.3203125, 77.7578125), Point(247.3203125, 68.08984375)),
    PathCurve(Point(247.3203125, 58.42578125), Point(255.15625, 50.58984375), Point(264.8203125, 50.58984375)),
    PathCurve(Point(274.484375, 50.58984375), Point(282.3203125, 58.42578125), Point(282.3203125, 68.08984375)),
    ]

    _green_dot = [
    PathMove(Point(261.30078125, 31.671875)),
    PathCurve(Point(261.30078125, 41.3359375), Point(253.46484375, 49.171875), Point(243.80078125, 49.171875)),
    PathCurve(Point(234.1328125, 49.171875), Point(226.30078125, 41.3359375), Point(226.30078125, 31.671875)),
    PathCurve(Point(226.30078125, 22.0078125), Point(234.1328125, 14.171875), Point(243.80078125, 14.171875)),
    PathCurve(Point(253.46484375, 14.171875), Point(261.30078125, 22.0078125), Point(261.30078125, 31.671875))
    ]

    if action == :clip
        map.(drawpath, (_j, _u, _l, _i, _a))
        map(drawpath, _blue_dot)
        map(drawpath, _green_dot)
        map(drawpath, _purple_dot)
        map(drawpath, _red_dot)
        do_action(:clip)
    elseif action == :path
        map.(drawpath, (_j, _u, _l, _i, _a))
        map(drawpath, _blue_dot)
        map(drawpath, _green_dot)
        map(drawpath, _purple_dot)
        map(drawpath, _red_dot)
    elseif color == true
        setcolor(bodycolor)
        map.(drawpath, (_j, _u, _l, _i, _a))
        do_action(action)

        sethue(Luxor.julia_blue)
        map(drawpath, _blue_dot)
        do_action(action)

        sethue(Luxor.julia_green)
        map(drawpath, _green_dot)
        do_action(action)

        sethue(Luxor.julia_purple)
        map(drawpath, _purple_dot)
        do_action(action)

        sethue(Luxor.julia_red)
        map(drawpath, _red_dot)
        do_action(action)
    else
        sethue(bodycolor)
        map.(drawpath, (_j, _u, _l, _i, _a))
        map(drawpath, _blue_dot)
        map(drawpath, _green_dot)
        map(drawpath, _purple_dot)
        map(drawpath, _red_dot)
        do_action(action)
    end
    # restore saved color
    setcolor(r, g, b, a)

    # and position
    if centered == true
        translate(Point(330, 213)/2)
    end
end

"""
    juliacircles(radius=100;
        outercircleratio=0.75,
        innercircleratio=0.65,
        action=:fill)

Draw the three Julia circles ("dots") in color centered at the origin.

The distance of the centers of each circle from the origin is `radius`.

The optional keyword argument `outercircleratio` (default 0.75) determines the radius of each circle relative to the main radius. So the default is to draw circles of radius 75 points around a larger circle of radius 100.

Return the three centerpoints.

The `innercircleratio` (default 0.65) no longer does anything useful (it used to draw the smaller circles) and will be deprecated.
"""
function juliacircles(radius=100;
        outercircleratio=0.75,
        innercircleratio=0.65,
        action=:fill)
    # clockwise, from bottom left
    color_sequence = [julia_red, julia_green, julia_purple]
    points = ngon(O, radius, 3, pi/6, vertices=true)
    @layer begin
        for (n, p) in enumerate(points)
            setcolor(color_sequence[n]...)
            circle(p, outercircleratio * radius, :path)
            if action != :fill
                newsubpath()
            else
                do_action(action)
            end
        end
        if action != :fill
            do_action(action)
        end
    end
    return points
end


text:
"""
    struct HomSpace{S<:ElementarySpace, P1<:CompositeSpace{S}, P2<:CompositeSpace{S}}
        codomain::P1
        domain::P2
    end

Represents the linear space of morphisms with codomain of type `P1` and domain of type `P2`.
Note that HomSpace is not a subtype of VectorSpace, i.e. we restrict the latter to denote
certain categories and their objects, and keep HomSpace distinct.
"""
struct HomSpace{S<:ElementarySpace, P1<:CompositeSpace{S}, P2<:CompositeSpace{S}}
    codomain::P1
    domain::P2
end

"""
    TensorSpace{S<:ElementarySpace}

A constant which is the union type of both elementary and product vector space and defined
as `HomSpace{S, ProductSpace{S, N₁}, ProductSpace{S, N₂}}`.
"""
const TensorSpace{S<:ElementarySpace} = Union{S, ProductSpace{S}}

"""
    TensorMapSpace{S<:ElementarySpace, N₁, N₂}

A constant which is an alias name of the HomSpace and defined as
`HomSpace{S, ProductSpace{S, N₁}, ProductSpace{S, N₂}}`.
"""
const TensorMapSpace{S<:ElementarySpace, N₁, N₂} =
    HomSpace{S, ProductSpace{S, N₁}, ProductSpace{S, N₂}}

"""
    →(dom::TensorSpace{S}, codom::TensorSpace{S}) where {S<:ElementarySpace}
    ←(codom::TensorSpace{S}, dom::TensorSpace{S}) where {S<:ElementarySpace}

Convenient constructor of HomSpace instance.
"""
→(dom::TensorSpace{S}, codom::TensorSpace{S}) where {S<:ElementarySpace} =
    HomSpace(ProductSpace(codom), ProductSpace(dom))

←(codom::TensorSpace{S}, dom::TensorSpace{S}) where {S<:ElementarySpace} =
    HomSpace(ProductSpace(codom), ProductSpace(dom))

# The index of a HomSpace is defined as the following: first, index spaces for codomain,
# then, index spaces for domain; dual is taken for each index space of domain, but the order
# are not reversed.
Base.getindex(W::TensorMapSpace{<:IndexSpace, N₁, N₂}, i) where {N₁, N₂} =
    i <= N₁ ? codomain(W)[i] : dual(domain(W)[i-N₁])

"""
    codomain(W::HomSpace)

Return the codomain of a HomSpace.
"""
codomain(W::HomSpace) = W.codomain

"""
    domain(W::HomSpace)

Return the domain of a HomSpace.
"""
domain(W::HomSpace) = W.domain

"""
    spacetype(W::HomSpace) -> Type{ElementarySpace}
    spacetype(::Type{<:HomSpace{S}})

Return the ElementarySpace type associated to a HomSpace instance or type.
"""
spacetype(W::HomSpace) = spacetype(typeof(W))
spacetype(::Type{<:HomSpace{S}}) where S = S

"""
    field(W::HomSpace) -> Field
    field(L::Type{<:HomSpace})

Return the field type over which a HomSpace instance or type is defined.
"""
field(W::HomSpace) = field(typeof(W))
field(L::Type{<:HomSpace}) = field(spacetype(L))

"""
    sectortype(W::HomSpace) -> Type{<:Sector}
    sectortype(L::Type{<:HomSpace)

Return the type of sector over which the HomSpace `W` is defined.
"""
sectortype(W::HomSpace) = sectortype(typeof(W))
sectortype(L::Type{<:HomSpace}) = sectortype(spacetype(L))

"""
    blocksectors(W::HomSpace)

Return an iterator over the different unique coupled sector labels, i.e. the intersection
of the different fusion outputs that can be obtained by fusing the sectors present in the
domain, as well as from the codomain.
"""
function blocksectors(W::HomSpace)
    sectortype(W) === Trivial &&
        return TrivialOrEmptyIterator(dim(domain(W)) == 0 || dim(codomain(W)) == 0)
    return intersect(blocksectors(codomain(W)), blocksectors(domain(W)))
end

"""
    dim(W::HomSpace)

Return the total dimension of a `HomSpace`, i.e. the number of linearly independent
morphisms that can be constructed within this space.
"""
function dim(W::HomSpace)
    d = 0
    for c in blocksectors(W)
        d += blockdim(codomain(W), c) * blockdim(domain(W), c)
    end
    return d
end

"""
    dual(W::HomSpace)

Return the dual of a HomSpace which contains the dual of morphisms in this space.
It corresponds to 180 degree rotation in the graphical representation.
"""
dual(W::HomSpace) = HomSpace(dual(W.domain), dual(W.codomain))

"""
    adjoint(W::HomSpace{<:EuclideanSpace})

Return the adjoint of a HomSpace which contains the dagger of morphisms in this space.
It corresponds to mirror operation and then reversing all arrows in the graphical
representation.
"""
Base.adjoint(W::HomSpace{<:EuclideanSpace}) = HomSpace(W.domain, W.codomain)

Base.hash(W::HomSpace, h::UInt) = hash(domain(W), hash(codomain(W), h))
Base.:(==)(W1::HomSpace, W2::HomSpace) =
    (W1.codomain == W2.codomain) && (W1.domain == W2.domain)

function Base.show(io::IO, W::HomSpace)
    if length(W.codomain) == 1
        print(io, W.codomain[1])
    else
        print(io, W.codomain)
    end
    print(io, " ← ")
    if length(W.domain) == 1
        print(io, W.domain[1])
    else
        print(io, W.domain)
    end
end

ext: jl


text:
q = ActiveQueue(6)
enqueue!(q, 1)
@test front(q) == 1
enqueue!(q, 2)
enqueue!(q, 3)
enqueue!(q, 4)
@test dequeue!(q) == 1
@test dequeue!(q) == 2
delete!(q, 4)
delete!(q, 3)
enqueue!(q, 5)
enqueue!(q, 6)
enqueue!(q, 4)
@test front(q) == 5
@test dequeue!(q) == 5
@test dequeue!(q) == 6
@test dequeue!(q) == 4
@test isempty(q)
enqueue!(q, 4)
enqueue!(q, 4)
delete!(q, 4)
@test isempty(q)


text:
using MPI
ENV["JULIA_MPI_BINARY"]="system"
using VertexEliminationOrder
using ArgParse

function parse_commandline(ARGS)
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--seed"
            help = "Base seed to use for sampling flow-cutter inputs."
            default = 42
            arg_type = Int64
        "--seed_bound"
            help = "Minimum difference of seed between MPI processes."
            default = 1000
            arg_type = Int64
        "--time"
            help = "The number of seconds to run flow cutter for."
            default = 30
            arg_type = Int64
        "--graph"
            help = "The name of the graph file .gr to open."
            default = "sycamore_53_20.gr"
            arg_type = String
        "--dir_graphfiles"
            help = "Path to your graph files .gr to read them."
            default = "$(@__DIR__)/../test/example_graphs/"
            arg_type = String
        "--nsample"
            help = "The number of sample of flow cutter random inputs per execution."
            default = 20
            arg_type = Int64

        "--sample_augment"
            help = "The number of execution needed before increasing nsample."
            default = 8
            arg_type = Int64
        "--max_imbalances"
            nargs ='*'
            help = "The max imabalances to select best cut from flow-cutter."
            default = [1.0, 0.8, 0.6]
            arg_type = Float64
    end

    return parse_args(ARGS, s)
end

parsed_args = parse_commandline(ARGS)

graph_name = parsed_args["graph"]
dirgraphfiles = parsed_args["dir_graphfiles"]
duration = parsed_args["time"]
max_imbalances = parsed_args["max_imbalances"]
base_seed = parsed_args["seed"]
seed_diff = parsed_args["seed_bound"]
base_nsample = parsed_args["nsample"]
sample_augment = parsed_args["sample_augment"]

println("- - - args - - -")
println("duration = $duration")
println("max_imbalances = $max_imbalances")
println("base_seed = $base_seed")
println("seed_diff = $seed_diff")
println("base_nsample = $base_nsample")
println("sample_augment = $sample_augment")
println("- - - result - - -")

MPI.Init()
comm = MPI.COMM_WORLD
my_rank = MPI.Comm_rank(comm)
root = 0
print("rank $my_rank has $(Threads.nthreads()) threads \n")

g = graph_from_gr(joinpath(dirgraphfiles, graph_name))
order, tw = order_tw_by_dissections_simple(
    g,
    duration;
    max_imbalances=max_imbalances ,
    seed=base_seed + my_rank * seed_diff,
    nsample=base_nsample,
    sample_augment=sample_augment
)
MPI.Barrier(comm)

"""returns the better (rank, treewidth) pair"""
function best(order1::Tuple{Int, Int}, order2::Tuple{Int, Int})
    order1[2] < order2[2] ? order1 : order2
end

rank_with_best, best_tw = MPI.Allreduce((my_rank, tw), best, comm)
MPI.Barrier(comm)

if my_rank == root
    best_order = similar(order)
    rreq = MPI.Irecv!(best_order, rank_with_best, rank_with_best, comm)
end

if my_rank == rank_with_best
    print("rank $my_rank sending elimination order to root \n")
    sreq = MPI.Isend(order, root, my_rank, comm)
end

MPI.Barrier(comm)

if my_rank == root
    @show best_tw
    @show best_order
end


text:
# This file was generated, do not modify it. # hide
using HTTP
using MLJ
using PyPlot
import DataFrames: DataFrame, describe
using UrlDownload
MLJ.color_off() # hide

url = "http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data"
header = ["Class", "Alcool", "Malic acid", "Ash", "Alcalinity of ash",
          "Magnesium", "Total phenols", "Flavanoids",
          "Nonflavanoid phenols", "Proanthcyanins", "Color intensity",
          "Hue", "OD280/OD315 of diluted wines", "Proline"]
data = urldownload(url, true, format=:CSV, header=header);


text:
# TODO: https://hal.inria.fr/inria-00558509v1/document
function liveness(ir)
  result = Dict(v => Set{Variable}() for v in keys(ir))
  result = merge(result, Dict(b.id => Set{Variable}() for b in blocks(ir)))
  queue = WorkQueue(1:length(blocks(ir)))
  while !isempty(queue)
    b = block(ir, pop!(queue))
    for br in branches(b)
      foreach(x -> x isa Variable && push!(result[b.id], x), arguments(br))
      br.condition isa Variable && push!(result[b.id], br.condition)
    end
    live = Set(result[b.id])
    for v in reverse(keys(b))
      delete!(live, v)
      varmap(x -> push!(live, x), b[v])
      union!(result[v], live)
    end
    setdiff!(live, arguments(b))
    for a in predecessors(b)
      isempty(setdiff(live, result[a.id])) && continue
      union!(result[a.id], live)
      push!(queue, a.id)
    end
  end
  return result
end

function interference(liveness)
  g = Graph{Variable}()
  for set in liveness
    set = collect(set)
    for i = 1:length(set), j = i+1:length(set)
      connect!(g, set[i], set[j])
    end
  end
  return g
end

interference(ir::IR) = interference(values(liveness(ir)))

function colouring(g::Graph, pre = Dict())
  nodes = sort(collect(keys(g)), by = x -> length(g[x]), rev = true)
  colours = Dict(n => get(pre, n, -1) for n in nodes)
  for n in nodes
    colour = 1
    while any(m -> colours[m] == colour, g[n])
      colour += 1
    end
    colours[n] = colour
  end
  return colours
end


text:
module BivMatFun

export fun2m

include("types.jl")
include("mpfr_tls.jl")
include("mp.jl")
include("utils.jl")

include("atoms.jl")
include("blocking.jl")
include("eigenvalues.jl")
include("fun2m.jl")
include("taylor.jl")

include("condest.jl")

end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
#================
Predictive
    Learning
        Models
================#
include("Distributions.jl")
@doc """
      |====== Lathe.models =====\n
      |____________/ Accessories ___________\n
      |_____models.Pipeline([steps])\n
      |____________/ Continuous models ___________\n
      |_____models.meanBaseline(y)\n
      |_____models.LinearRegression(x,y)\n
      |_____models.LeastSquare(x,y,Type)\n
      |_____models.PowerLog(prob1,prob2)\n
      |____________/ Categorical Models ___________\n
      |_____models.LogisticRegression(x,y)\n
      |_____models.majBaseline(y)\n
       """
module models
#==
Base
    Models
        Functions
==#
using Lathe
using Random
using DataFrames
#===========
Accessories
===========#
@doc """
      Pipelines can contain a predictable Lathe model with preprocessing that
      occurs automatically. This is done by putting X array processing methods
      into the iterable steps, and then putting your Lathe model in.\n
      --------------------\n
      ==PARAMETERS==\n
      [steps] <- An iterable list of methods to call for X modification. These mutations should
      have ALREADY BEEN MADE TO THE TRAIN X.\n
      pipl = Pipeline([StandardScalar(),LinearRegression(trainX,trainy)])\n
      --------------------\n
      ==Functions==\n
      predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
      """
function Pipeline(steps)
    predict(xt) = [object.predict(xt) for object in steps]
    (var)->(steps;model;predict)
end
#==============
========================================================
=======================================================================
            CONTINUOS MODELS               CONTINUOS MODELS
            CONTINUOS MODELS               CONTINUOS MODELS
======================================================================
======================================================================#
#==
Mean
    Baseline
==#
 # Model Type
 @doc """
       A mean baseline is great for getting a basic accuracy score in order
           to make a valid direction for your model.\n
         --------------------\n
         ==PARAMETERS==\n
        [y] <- Fill with your trainY values. Should be an array of shape (0,1) or (1,0)\n
        pipl = Pipeline([StandardScalar(),LinearRegression(trainX,trainy)])\n
        --------------------\n
        ==Functions==\n
        predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
                     """
function MeanBaseline(y)
    m = Lathe.stats.mean(m.y)
    predict(xt) =
    xt = [v = m for v in xt]
    (test)->(m;predict)
end
#==
Linear
    Regression
==#

@doc """
      Linear Regression is a well-known linear function used for predicting
      continuous features with a mostly linear or semi-linear slope.\n
      --------------------\n
      ==PARAMETERS==\n
     [y] <- Fill with your trainY values. Should be an array of shape (0,1) or (1,0)\n
     [x] <- Fill in with your trainX values. Should be an array of shape (0,1) or (1,0)\n
     --------------------\n
     ==Functions==\n
     predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
       """
function LinearRegression(x,y)
    # a = ((∑y)(∑x^2)-(∑x)(∑xy)) / (n(∑x^2) - (∑x)^2)
    # b = (x(∑xy) - (∑x)(∑y)) / n(∑x^2) - (∑x)^2
    if length(x) != length(y)
        throw(ArgumentError("The array shape does not match!"))
    end
    # Get our Summations:
    Σx = sum(x)
    Σy = sum(y)
    # dot x and y
    xy = x .* y
    # ∑dot x and y
    Σxy = sum(xy)
    # dotsquare x
    x2 = x .^ 2
    # ∑ dotsquare x
    Σx2 = sum(x2)
    # n = sample size
    n = length(x)
    # Calculate a
    a = (((Σy) * (Σx2)) - ((Σx * (Σxy)))) / ((n * (Σx2))-(Σx^2))
    # Calculate b
    b = ((n*(Σxy)) - (Σx * Σy)) / ((n * (Σx2)) - (Σx ^ 2))
    # The part that is super struct:
    predict(xt) = (xt = [i = a + (b * i) for i in xt])
    (test)->(a;b;predict)
end
#==
Linear
    Least
     Square
==#
@doc """
      Least Squares is ideal for predicting continous features.
      Many models use Least Squares as a base to build off of.\n
      --------------------\n
      x = [7,6,5,6,5]\n
      y  = [3.4.5.6.3]\n
      xtrain = [7,5,4,5,3,5,7,8]\n
      Type = :LIN\n
      model = Lathe.models.LeastSquare(x,y,Type)\n
      y_pred = Lathe.models.predict(model,xtrain)\n
      -------------------\n
      HYPER PARAMETERS\n
      Type <- Type determines which Linear Least Square algorithm to use,
      :LIN, :OLS, :WLS, and :GLS are the three options.\n
      - :LIN = Linear Least Square Regression\n
      - :OLS = Ordinary Least Squares\n
      - :WLS = Weighted Least Squares\n
      - :GLS = General Least Squares
      --------------------\n
      ==Functions==\n
      predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
       """
function LeastSquare(x,y,Type)
    if length(x) != length(y)
        throw(ArgumentError("The array shape does not match!"))
    end
    if Type == :LIN
        xy = x .* y
        sxy = sum(xy)
        n = length(x)
        x2 = x .^ 2
        sx2 = sum(x2)
        sx = sum(x)
        sy = sum(y)
        # Calculate the slope:
        slope = ((n*sxy) - (sx * sy)) / ((n * sx2) - (sx)^2)
     # Calculate the y intercept
        b = (sy - (slope*sx)) / n
    end
    predict(xt) =
    if Type == :LIN
        (xt = [z = (slope * x) + b for x in xt])
    end
    (test)->(slope;b;predict)
end
#==
Ridge
    Regression
==#
@doc """
      Ridge Regression is another regressor ideal for predicting linear,
          continuous features.\n
      --------------------\n
      array = [5,10,15]\n
      scaled_feature = Lathe.preprocess.OneHotEncode(array)\n
       """
function RidgeRegression(x,y)

end

#======================================================================
=======================================================================
            CATEGORICAL MODELS             CATEGORICAL MODELS
            CATEGORICAL MODELS             CATEGORICAL MODELS
======================================================================
======================================================================#
#==
Majority
    Class
        Baseline
==#
@doc """
      Majority class baseline is used to find the most often interpreted
      classification in an array.\n
      --------------------\n
      ==PARAMETERS==\n
     [y] <- Fill with your trainY values. Should be an array of shape (0,1) or (1,0)\n
     --------------------\n
     ==Functions==\n
     predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)\n
     counts() <- Returns a dictionary with the counts of all inserted keys.\n
     highest() <- Will return a Dictionary key with the count as well as the value for the most interpreted classification.
       """
function majClassBaseline(y)
    u=unique(y)
    d=Dict([(i,count(x->x==i,y)) for i in u])
    d = sort(collect(d), by=x->x[2])
    maxkey = d[length(d)]
    predict(xt) = [p = maxkey[1] for p in xt]
    counts() = d
    highest() = maxkey
    (var)->(y;maxkey;d;predict;counts;highest)
end
#==
Power Log
==#
@doc """
      A powerlog can be used to perform a one-tailed test, as well as get the proper sample size for a testing population.\n
      --------------------\n
      ==PARAMETERS==\n
     p1 <- A Float64 percentage representing the probability of scenario one.\n
     p2 <- A Float64 percentage representing the probability of scenario two. These two probability values should follow these guidelines: p1 = p1 + x = p2\n
     alpha = 0.05 <- Sets an alpha value\n
     --------------------\n
     Returns power, sample_size
       """
function PowerLog(p1::Float64,p2::Float64; alpha::Float64 = 0.05, rsq::Real = 0)
    pd = p2 - p1
    l1 = p1/(1-p1)
    l2 = p2/(1-p2)
    θ = l2 / l1
    or = θ
    λ = log(θ)
    λ2 = λ^2
    za = quantile(normal(),1-alpha)
    println("One-tailed test: alpha = ",alpha,", p1 = ",p1,", p2 = ",p2,", rsq = ",rsq,", odds ratio = ",or)
    δ = (1 + (1 + λ2)*exp(5 * λ2/4))/(1 + exp(-1*λ2/4))
    pwr = zeros(Float64,8)
    nn = zeros(Int64,8)
    i = 1
    for power = 0.6:.05:.95
        zb = quantile(normal(),power)

        N = ((za + zb*exp(-1 * λ2/4))^2 * (1 + 2*p1*δ))/(p1*λ2)
        N /= (1 - rsq)
        pwr[i] = power
        nn[i] = ceil(Int64,N)
        i += 1
    end
    return(pwr, nn)
end
#==
Multinomial
    Naive
        Bayes
==#
function MultinomialNB(x,y)

end
#==
Logistic
    Regression
==#
@doc """
      Majority class baseline is used to find the most often interpreted
      classification in an array.\n
      --------------------\n
      ==PARAMETERS==\n
     [X] <- Fill with your trainX values. Should be an array of shape (0,1) or (1,0)\n
     [y] <- Fill with your trainy values. Should be an array of shape (0,1) or (1,0)\n
     λ = .0001 <- Lambda Value\n
     fit_intercept = true <- Boolean determines whether to fit an intercept.\n
     max_iter = 1000 <- Determines the maximum number of iterations for the model to perform.\n
     --------------------\n
     ==Functions==\n
     predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)\n
       """
function LogisticRegression(X, y, λ=0.0001, fit_intercept=true, η=0.01, max_iter=1000)
    θ, 𝐉 = logistic_regression_sgd(X, y, 0.0001, true, 0.3, 3000);
    predict(xt) = yhat = predict_class(predict_proba(xt,0))
    cost = 𝐉
    (var) -> (predict;cost)
end
function sigmoid(z)
    return 1 ./ (1 .+ exp.(.-z))
end
function logistic_regression_sgd(X, y, λ, fit_intercept=true, η=0.01, max_iter=1000)

    # Initialize some useful values
    m = length(y); # number of training examples

    if fit_intercept
        # Add a constant of 1s if fit_intercept is specified
        constant = ones(m, 1)
        X = hcat(constant, X)
    else
        X # Assume user added constants
    end

    # Use the number of features to initialise the theta θ vector
    n = size(X)[2]
    θ = zeros(n)

    # Initialise the cost vector based on the number of iterations
    𝐉 = zeros(max_iter)

    for iter in range(1, stop=max_iter)

        # Calcaluate the cost and gradient (∇𝐉) for each iter
        𝐉[iter], ∇𝐉 = regularised_cost(X, y, θ, λ)

        # Update θ using gradients (∇𝐉) for direction and (η) for the magnitude of steps in that direction
        θ = θ - (η * ∇𝐉)
    end

    return (θ, 𝐉)
end
function regularised_cost(X, y, θ, λ)
    m = length(y)
    h = sigmoid(X * θ)
    positive_class_cost = ((-y)' * log.(h))
    negative_class_cost = ((1 .- y)' * log.(1 .- h))
    lambda_regularization = (λ/(2*m) * sum(θ[2 : end] .^ 2))
    𝐉 = (1/m) * (positive_class_cost - negative_class_cost) + lambda_regularization
    ∇𝐉 = (1/m) * (X') * (h-y) + ((1/m) * (λ * θ))
    ∇𝐉[1] = (1/m) * (X[:, 1])' * (h-y)
           return (𝐉, ∇𝐉)
end
function predict_proba(X, θ, fit_intercept=true)
    m = size(X)[1]

    if fit_intercept
        constant = ones(m, 1)
        X = hcat(constant, X)
    else
        X
    end

    h = sigmoid(X * θ)
    return h
end

function predict_class(proba, threshold=0.5)
    return proba .>= threshold
end
#==
Nueral
    Network
        Framework
==#
function calculate_activation_forward(A_pre , W , b , function_type)
    if (function_type == "sigmoid")
        Z , linear_step_cache = forward_linear(A_pre , W , b)
        A , activation_step_cache = sigmoid(Z)
    elseif (function_type == "relu")
        Z , linear_step_cache = forward_linear(A_pre , W , b)
        A , activation_step_cache = relu(Z)
    end
    cache = (linear_step_cache , activation_step_cache)
    return A , cache
end

function model_forward_step(X , params)
    all_caches = []
    A = X
    L = length(params)/2
    for l=1:L-1
        A_pre = A
        A , cache = calculate_activation_forward(A_pre , params[string("W_" , string(Int(l)))] , params[string("b_" , string(Int(l)))] , "relu")
        push!(all_caches , cache)
    end
    A_l , cache = calculate_activation_forward(A , params[string("W_" , string(Int(L)))] , params[string("b_" , string(Int(L)))] , "sigmoid")
    push!(all_caches , cache)
    return A_l , all_caches
end
function cost_function(AL , Y)
    cost = -mean(Y.*log.(AL) + (1 .- Y).*log.(1 .- AL))
    return cost
end
function forward_linear(A , w , b)
    Z = w*A .+ b
    cache = (A , w , b)
    return Z,cache
end
function init_param(layer_dimensions)
    param = Dict()
    for l=0:length(layer_dimensions)
        param[string("W_" , string(l-1))] = rand(layer_dimensions[l] ,
				layer_dimensions[l-1])*0.1
        param[string("b_" , string(l-1))] = zeros(layer_dimensions[l] , 1)
    end
    return param
end

function backward_linear_step(dZ , cache)
    A_prev , W , b = cache
    m = size(A_prev)[2]
    dW = dZ * (A_prev')/m
    db = sum(dZ , dims = 2)/m
    dA_prev = (W')* dZ
    return dW , db , dA_prev
end
function relu(X)
    rel = max.(0,X)
    return rel , X
end
function backward_relu(dA , cache_activation)
    return dA.*(cache_activation.>0)
end
function backward_sigmoid(dA , cache_activation)
    return dA.*(sigmoid(cache_activation)[1].*(1 .- sigmoid(cache_activation)[1]))
end
function backward_activation_step(dA , cache , activation)
    linear_cache , cache_activation = cache
    if (activation == "relu")
        dZ = backward_relu(dA , cache_activation)
        dW , db , dA_prev = backward_linear_step(dZ , linear_cache)
    elseif (activation == "sigmoid")
        dZ = backward_sigmoid(dA , cache_activation)
        dW , db , dA_prev = backward_linear_step(dZ , linear_cache)
    end
    return dW , db , dA_prev

end
function model_backwards_step(A_l , Y , caches)
    grads = Dict()
    L = length(caches)
    m = size(A_l)[2]
    Y = reshape(Y , size(A_l))
    dA_l = (-(Y./A_l) .+ ((1 .- Y)./( 1 .- A_l)))
    current_cache = caches[L]
    grads[string("dW_" , string(L))] , grads[string("db_" , string(L))] , grads[string("dA_" , string(L-1))] = backward_activation_step(dA_l , current_cache , "sigmoid")
    for l=reverse(0:L-2)
        current_cache = caches[l+1]
        grads[string("dW_" , string(l+1))] , grads[string("db_" , string(l+1))] , grads[string("dA_" , string(l))] = backward_activation_step(grads[string("dA_" , string(l+1))] , current_cache , "relu")

    end
    return grads
end

function update_param(parameters , grads , learning_rate)
    L = Int(length(parameters)/2)
    for l=0:(L-1)
        parameters[string("W_" , string(l+1))] -= learning_rate.*grads[string("dW_" , string(l+1))]
        parameters[string("b_",string(l+1))] -= learning_rate.*grads[string("db_",string(l+1))]
    end
    return parameters
end
@doc """
      'Network' specifically constructs a convolutional nueral network,\n
      though this network is still being worked on in both methodology\n
      and in a physical manner,\n thanks for understanding.
       """
function Network(X,Y,layers_dimensions,n_iter)
    params = init_param(layers_dimensions)
    costs = []
    iters = []
    accuracy = []
    for i=1:n_iter
        A_l , caches  = model_forward_step(X , params)
        cost = cost_function(A_l , Y)
        acc = check_accuracy(A_l , Y)
        grads  = model_backwards_step(A_l , Y , caches)
        params = update_param(params , grads , learning_rate)
        println("Iteration ->" , i)
        println("Cost ->" , cost)
        println("Accuracy -> " , acc)
        push!(iters , i)
        push!(costs , cost)
        push!(accuracy , acc)
    end
    predict(xt) = (xt + 4)
    (test)->(costs;params;predict)
end

#
#----------------------------------------------
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule wolfSSL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("wolfSSL")
JLLWrappers.@generate_main_file("wolfSSL", UUID("98c43586-9870-5ae5-ab22-acc77b9bbdb5"))
end  # module wolfSSL_jll


text:
using EuclidianNormalizingFlows
using Test

using Distributions
using InverseFunctions, ChangesOfVariables
using ForwardDiff

using EuclidianNormalizingFlows: JohnsonSU, JohnsonTrafo, JohnsonTrafoInv, johnsontrafo, johnsontrafo_inv, johnsontrafo_ladj, johnsontrafo_inv_ladj

@testset "JohnsonSU" begin
    K = rand(Normal(), 10^6)
    K_joh = johnsontrafo_inv.(K, -15,6.5,0,2.5);
    X = rand(JohnsonSU(-15,6.5,0,2.5), 10^6);

    @test isapprox(sum(abs.(sort(K_joh)))/10^6, sum(abs.(sort(X)))/10^6, rtol = 0.01)

    @test @inferred(johnsontrafo(0.5, 1, 2, 3, 4)) isa Float64
    @test @inferred(johnsontrafo_inv(0.5, 1, 2, 3, 4)) isa Float64

    @test johnsontrafo(0.3, 1, 3, -4, 0.5) ≈ 9.544817734776984
    @test johnsontrafo_inv(0.3, 1, 3, -4, 0.5) ≈ -4.1177281942392545

    Z = johnsontrafo_inv.(K, -2, 1, 0, 2.5)
    Z_reco = johnsontrafo.(Z, -2, 1, 0, 2.5)
    @test K ≈ Z_reco

    @test isapprox(@inferred(johnsontrafo_ladj(0.5, 4.2, 4, 2, 3)), log(abs(ForwardDiff.derivative(x -> johnsontrafo(x, 4.2, 4, 2, 3), 0.5))), rtol = 0.01)
    @test isapprox(@inferred(johnsontrafo_inv_ladj(0.5, 4.2, 4, 2, 3)), log(abs(ForwardDiff.derivative(x -> johnsontrafo_inv(x, 4.2, 4, 2, 3), 0.5))), rtol = 0.01)

    @test @inferred(JohnsonSU(4, 2, 3, 1)) isa JohnsonSU

    @test @inferred(JohnsonTrafo(4, 2, 3, 1)) isa JohnsonTrafo
    trafo = JohnsonTrafo(4, 2, 3, 1)
    @test @inferred(inverse(trafo)) isa JohnsonTrafoInv
    @test @inferred(inverse(inverse(trafo))) === trafo

    @test @inferred(trafo(4.2)) == johnsontrafo(4.2, 4, 2, 3, 1)

    @test @inferred(with_logabsdet_jacobian(JohnsonTrafo([4.0, 4.1], [3.0, 3.1], [2.0, 2.1], [1.0, 1.1]), [0.5, 0.6])) == (
        johnsontrafo.([0.5, 0.6], [4.0, 4.1], [3.0, 3.1], [2.0, 2.1], [1.0, 1.1]),
        sum(johnsontrafo_ladj.([0.5, 0.6], [4.0, 4.1], [3.0, 3.1], [2.0, 2.1], [1.0, 1.1])),
    )

    @test @inferred(with_logabsdet_jacobian(JohnsonTrafoInv([4.0, 4.1], [3.0, 3.1], [2.0, 2.1], [1.0, 1.1]), [0.5, 0.6])) == (
        johnsontrafo_inv.([0.5, 0.6], [4.0, 4.1], [3.0, 3.1], [2.0, 2.1], [1.0, 1.1]),
        sum(johnsontrafo_inv_ladj.([0.5, 0.6], [4.0, 4.1], [3.0, 3.1], [2.0, 2.1], [1.0, 1.1])),
    )

    let
        fwd = JohnsonTrafo([10.0, 11.0], [3.5, 3.6], [10.0, 11.0], [1.0, 1.1])
        @test fwd == deepcopy(fwd)
        @test isequal(fwd, deepcopy(fwd))
        @test hash(fwd) == hash(deepcopy(fwd))

        InverseFunctions.test_inverse(fwd, [2, 3])
        InverseFunctions.test_inverse(fwd, [-2,3])
        ChangesOfVariables.test_with_logabsdet_jacobian(fwd, [2, 3], ForwardDiff.jacobian)
        ChangesOfVariables.test_with_logabsdet_jacobian(fwd, [-2, -3], ForwardDiff.jacobian)

        rev = inverse(fwd)
        @test rev == deepcopy(rev)
        @test isequal(rev, deepcopy(rev))
        @test hash(rev) == hash(deepcopy(rev))

        InverseFunctions.test_inverse(rev, [2, 3])
        InverseFunctions.test_inverse(rev, [-2,3])
        ChangesOfVariables.test_with_logabsdet_jacobian(rev, [3, 4], ForwardDiff.jacobian)
        ChangesOfVariables.test_with_logabsdet_jacobian(rev, [-3, -4], ForwardDiff.jacobian)

        X = randn(2, 3)
        Y, ladjs = with_logabsdet_jacobian(fwd, X)
        @test hcat((getindex).(with_logabsdet_jacobian.(fwd, eachcol(X)), 1)...) == Y
        @test (getindex).(with_logabsdet_jacobian.(fwd, eachcol(X)), 2) == ladjs'
        X2, inv_ladjs = with_logabsdet_jacobian(rev, Y)
        @test X2 ≈ X
        @test inv_ladjs ≈ - ladjs
    end
end


text:
# -------------------------------------------------
#
# An explicit tensor constructor of the coefficients
# for \prod_i f_i (x_i) given the coefficients of
# each f(x_i).
#
# -------------------------------------------------

# This is a useful way to reduce a difficult O(N log^(d-1) N)
# space of integrations down to O(N d)
# for specifying coefficients of initial conditions

# Efficiency criticality: MEDIUM

# Accuracy criticality: LOW
# Most of the accuracy dependent on functions
# called from other scripts

function tensor_construct(D::Int, k::Int, n::Int, coeff_array::Array{Dict{CartesianIndex{1},
    Array{Array{T, 1}, 1}}, 1}; scheme = "sparse") where T <: Real
    tensor_construct(Val(D), k, n, coeff_array, Val(Symbol(scheme)))
end
function tensor_construct(::Val{D}, k::Int, n::Int, coeff_array::Array{Dict{CartesianIndex{1},
    Array{Array{T, 1}, 1}}, 1}, scheme::Val{Scheme}) where {D, T <: Real, Scheme}

    coeffs = Dict{CartesianIndex{D}, Array{Array{T,D},D}}()
    ls::NTuple{D, Int}    = ntuple(i-> (n+1), D)
    modes::NTuple{D, Int} = ntuple(i-> k, D)

    for level in CartesianIndices(ls)
        cutoff(scheme, level, n) && continue

        cells::NTuple{D, Int} = ntuple(i -> 1<<max(0, level[i]-2), D)
        level_coeffs = Array{Array{T,D}}(undef, cells)
        for cell in CartesianIndices(cells)
            cell_coeffs = Array{T}(undef, modes)
            for mode in CartesianIndices(modes)
                val = one(T)
                for d in 1:D
                    coeff = coeff_array[d]
                    tup::NTuple{1, Int} = (level[d],) # This is somehow still slow
                    l = CartesianIndex{1}(tup);
                    c = cell[d]
                    m = mode[d]
                    val *= coeff[l][c][m]
                end
                cell_coeffs[mode] = val
            end
            level_coeffs[cell] = cell_coeffs
        end
        coeffs[level] = level_coeffs
    end
    return coeffs
end

function tensor_construct(D::Int, k::Int, n::Int, vcoeff_array::Array{Array{T, 1}};
							scheme = "sparse") where T <: Real

	coeff_array = [V2D(1, k, n, vcoeff_array[i]) for i in 1:length(vcoeff_array)]
	return D2V(D, k, n, tensor_construct(D, k, n, coeff_array); scheme=scheme)

end


text:
module Colormaps

using PyCall
@pyimport importlib
using PyPlot
@pyimport matplotlib.colors as COL

data_dir = joinpath(Pkg.dir("Colormaps"), "data")
for file in filter(f -> ismatch(r"\.py$", f), readdir(data_dir))
    cmap_name = replace(rsplit(file, ".", limit=2)[1], ".", "_")
    mod = importlib.machinery[:SourceFileLoader](cmap_name,
                                                 joinpath(data_dir,file))[:load_module]()
    try
        plt[:cm][:get_cmap](cmap_name)
    catch
        cmap = COL.ListedColormap(mod[:cm_data], name=cmap_name)
        register_cmap(cmap_name, cmap)
    end
end

end # module


text:
#=
Pentagonal numbers are generated by the formula, Pn=n(3n−1)/2. The first ten pentagonal numbers are:

1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...

It can be seen that P4 + P7 = 22 + 70 = 92 = P8. However, their difference, 70 − 22 = 48, is not pentagonal.

Find the pair of pentagonal numbers, Pj and Pk, for which their sum and difference are pentagonal and D = |Pk − Pj| is minimised; what is the value of D?
=#
function pent(x)
  div(x*(3x-1),2)
end

function ispent(x)
  n = (sqrt(24x+1)+1)/6
  return n == round(Int,n) && n > 0
end

function calc()
  i = 0
  pentagonals = Int64[]
  while true
    i += 1
    p = pent(i)
    for j in pentagonals
      diff = abs(p-j)
      if ispent(diff)
        intsum = p + j
        if ispent(intsum)
          return diff
        end
      end
    end
    push!(pentagonals, p)
  end
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
"""
    struct ConvexRealSet
Convex real set (similar to a line segment).
Parameters: lower bound, upper bound: floating point numbers
            includelower, includeupper: boolean true or false to indicate whether
            the set has a closed boundary (set to true) or open (set to false).
"""
mutable struct ConvexRealSet
    lower::Float64
    includelower::Bool
    upper::Float64
    includeupper::Bool
    function ConvexRealSet(lo, up, incllo, inclup)
       this = new()
       this.upper = Float64(up)
       this.lower = Float64(lo)
       this.includelower = incllo
       this.includeupper = inclup
       this
    end
end

function ∈(s, xelem)
    x = Float64(xelem)
    if(x == s.lower)
        if(s.includelower)
            return true
        else
            return false
        end
    elseif(x == s.upper)
        if(s.includeupper)
            return true
        else
            return false
        end
    end
    s.lower < x && x < s.upper
end

⋃(aset, bset, x) = (∈(aset, x) || ∈(bset, x))

⋂(aset, bset, x) = (∈(aset, x) && ∈(bset, x))

-(aset, bset, x) = (∈(aset, x) && !∈(bset, x))

isempty(s::ConvexRealSet) = (s.lower > s.upper) ||
                           ((s.lower == s.upper) && !s.includeupper && !s.includelower)

const s1 = ConvexRealSet(0.0, 1.0, false, true)
const s2 = ConvexRealSet(0.0, 2.0, true, false)
const s3 = ConvexRealSet(1.0, 2.0, false, true)
const s4 = ConvexRealSet(0.0, 3.0, true, false)
const s5 = ConvexRealSet(0.0, 1.0, false, false)
const s6 = ConvexRealSet(0.0, 1.0, true, true)
const sempty = ConvexRealSet(0.0, -1.0, true, true)
const testlist = [0, 1, 2]

function testconvexrealset()
    for i in testlist
        println("Testing with x = $i.\nResults:")
        println("    (0, 1] ∪ [0, 2): $(⋃(s1, s2, i))")
        println("    [0, 2) ∩ (1, 2]: $(⋂(s2, s3, i))")
        println("    [0, 3) − (0, 1): $(-(s4, s5, i))")
        println("    [0, 3) − [0, 1]: $(-(s4, s6, i))\n")
    end
    print("The set sempty is ")
    println(isempty(sempty) ? "empty." : "not empty.")
end

testconvexrealset()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using UnPack

export  Bilayer,
        bilayer!

exppi(x) = cospi(x) + im*sinpi(x)

function phase(G,τ)
    if size(τ,2) == 1
        return exppi(G'*τ)
    else
        return exppi.(G'*τ)
    end
end

function bilayer!(kx,ky,V!,lattice₁,lattice₂,Moire,Matrix)
    @unpack fold, G₁, G₂, numfold, M = Moire
    @. Matrix = 0.
    G₂ = ( G₁/lattice₂.b .|> x -> round(x,RoundNearest) ) |> x -> x*lattice₂.b # 1/2は近い偶数になる
    @inbounds Threads.@threads for g in 1:numfold
        from = g*lattice₂.M+1 ; to = g*lattice₂.M+lattice₁.M
        V!(kx+G₁[g,1],ky+G₁[g,2],view(Matrix,1:lattice₁.M,from:to))
        rmul!( view(Matrix ,1:lattice₁.M,from:to) , phase( view(-G₂,g,:) , lattice₂.τ ) )
        lmul!( phase( view(G₁,g,:) , lattice₁.τ ) , view(Matrix ,1:lattice₁.M,from:to) )
    end
    @. Matrix += Matrix'
    from  = 1; to = lattice₁.M
    lattice₁.Ham!(kx,ky,view(Matrix,from:to,from:to))
    @inbounds Threads.@threads for g in 1:numfold
        from = g*lattice₂.M+1 ; to = g*lattice₂.M+lattice₁.M
        lattice₂.Ham!(kx+G₁[g,1],ky+G₁[g,2],view(Matrix,from:to,from:to))
    end
end

struct Bilayer
    Ham!
    function Bilayer(V,lattice₁,lattice₂,Moire,inplace)
        @unpack fold, G₁, G₂, numfold, M = Moire
        if inplace
            function _bilayer!(kx,ky,Matrix)
                @. Matrix = 0.
                G₂ = ( G₁/lattice₂.b .|> x -> round(x,RoundNearest) ) |> x -> x*lattice₂.b # 1/2は近い偶数になる
                @inbounds for g in 1:numfold
                    from = lattice₁.M+(g-1)*lattice₂.M+1 ; to = lattice₁.M+g*lattice₂.M
                    V(kx+G₁[g,1],ky+G₁[g,2],view(Matrix,1:lattice₁.M,from:to))
                    rmul!( view(Matrix ,1:lattice₁.M,from:to) , phase( view(-G₂,g,:) , lattice₂.τ ) )
                    lmul!( phase( view(G₁,g,:) , lattice₁.τ ) , view(Matrix ,1:lattice₁.M,from:to) )
                end
                @. Matrix += Matrix'
                from  = 1; to = lattice₁.M
                lattice₁.Ham!(kx,ky,view(Matrix,from:to,from:to))
                @inbounds for g in 1:numfold
                    from = lattice₁.M+(g-1)*lattice₂.M+1 ; to = lattice₁.M+g*lattice₂.M
                    lattice₂.Ham!(kx+G₁[g,1],ky+G₁[g,2],view(Matrix,from:to,from:to))
                end
            end
            return new(_bilayer!)

        else
            function _bilayer(kx,ky)
                Matrix = zeros(ComplexF64,lattice₁.M+numfold*lattice₂.M,lattice₁.M+numfold*lattice₂.M)
                G₂ = ( G₁/lattice₂.b .|> x -> round(x,RoundNearest) ) |> x -> x*lattice₂.b # 1/2は近い偶数になる
                @inbounds for g in 1:numfold
                    from = lattice₁.M+(g-1)*lattice₂.M+1 ; to = lattice₁.M+g*lattice₂.M
                    rmul!( view(Matrix ,1:lattice₁.M,from:to) , phase( view(-G₂,g,:) , lattice₂.τ ) )
                    lmul!( phase( view(G₁,g,:) , lattice₁.τ ) , view(Matrix ,1:lattice₁.M,from:to) )
                end
                @. Matrix += Matrix'
                from  = 1; to = lattice₁.M
                Matrix[from:to,from:to] = lattice₁.Ham!(kx,ky)
                @inbounds for g in 1:numfold
                    from = lattice₁.M+(g-1)*lattice₂.M+1 ; to = lattice₁.M+g*lattice₂.M
                    Matrix[from:to,from:to] = lattice₂.Ham!(kx+G₁[g,1],ky+G₁[g,2])
                end
                Matrix
            end
            return new(_bilayer)
        end
    end
end


text:
module TightBindingApproximation

using Printf: @sprintf
using TimerOutputs: @timeit
using Optim: optimize, LBFGS
using RecipesBase: RecipesBase, @recipe, @series
using LinearAlgebra: inv, dot, Hermitian, Diagonal, eigvals, cholesky, Eigen
using QuantumLattices: getcontent, iidtype, rcoord, icoord, expand, statistics, plain, creation, annihilation, atol, rtol, periods, decimaltostr
using QuantumLattices: AbstractLattice, AbstractPID, FID, NID, Index, CompositeOID, ID, Bonds, Hilbert, Metric, Operator, Operators, OIDToTuple, Table, Term, Boundary
using QuantumLattices: Internal, Fock, Phonon, Hopping, Onsite, Pairing, PhononKinetic, PhononPotential, MatrixRepresentation, BrillouinZone
using QuantumLattices: Engine, Parameters, AbstractGenerator, CompositeGenerator, Entry, Generator, Formulation, Action, Assignment, Algorithm

import LinearAlgebra: eigen, eigvals, ishermitian
import QuantumLattices: add!, contentnames, dimension, kind, matrix, update!, prepare!, run!

export TBAKind, AbstractTBA, TBAMatrix, TBAMatrixRepresentation, commutator
export TBA, EnergyBands, BerryCurvature
export SampleNode, deviation, optimize!

"""
    TBAKind{K}

The kind of a free quantum lattice system using the tight-binding approximation.
"""
struct TBAKind{K}
    TBAKind(k::Symbol) = new{k}()
end
@inline TBAKind{K}() where K = TBAKind(K)
@inline Base.promote_rule(::Type{TBAKind{K}}, ::Type{TBAKind{K}}) where K = TBAKind{K}
@inline Base.promote_rule(::Type{TBAKind{:TBA}}, ::Type{TBAKind{:BdG}}) = TBAKind{:BdG}
@inline Base.promote_rule(::Type{TBAKind{:BdG}}, ::Type{TBAKind{:TBA}}) = TBAKind{:BdG}

"""
    TBAKind(T::Type{<:Term})

Depending on the kind of a term type, get the corresponding TBA kind.
"""
@inline TBAKind(::Type{T}) where {T<:Term} = error("TBAKind error: not defined for $(kind(T)).")
@inline TBAKind(::Type{T}) where {T<:Union{Hopping, Onsite}} = TBAKind(:TBA)
@inline TBAKind(::Type{T}) where {T<:Union{Pairing, PhononKinetic, PhononPotential}} = TBAKind(:BdG)
@inline @generated function TBAKind(::Type{TS}) where {TS<:Tuple{Vararg{Term}}}
    exprs = []
    for i = 1:fieldcount(TS)
        push!(exprs, :(typeof(TBAKind(fieldtype(TS, $i)))))
    end
    return Expr(:call, Expr(:call, :reduce, :promote_type, Expr(:tuple, exprs...)))
end

"""
    Metric(::TBAKind, hilbert::Hilbert{<:Fock} -> OIDToTuple
    Metric(::TBAKind, hilbert::Hilbert{<:Phonon}) -> OIDToTuple

Get the oid-to-tuple metric for a free fermionic/bosonic system or a free phononic system.
"""
@inline @generated Metric(::TBAKind{:TBA}, hilbert::Hilbert{<:Fock}) = OIDToTuple(fieldnames(keytype(hilbert))..., :orbital, :spin)
@inline @generated Metric(::TBAKind{:BdG}, hilbert::Hilbert{<:Fock}) = OIDToTuple(:nambu, fieldnames(keytype(hilbert))..., :orbital, :spin)
@inline @generated Metric(::TBAKind, hilbert::Hilbert{<:Phonon}) = OIDToTuple(:tag, fieldnames(keytype(hilbert))..., :dir)

"""
    commutator(k::TBAKind, hilbert::Hilbert{<:Internal}) -> Union{AbstractMatrix, Nothing}

Get the commutation relation of the single-particle operators of a free quantum lattice system using the tight-binding approximation.
"""
@inline commutator(::TBAKind{:TBA}, ::Hilbert{<:Internal}) = nothing
@inline commutator(::TBAKind{:BdG}, ::Hilbert{<:Fock{:f}}) = nothing
@inline commutator(k::TBAKind{:BdG}, hilbert::Hilbert{<:Fock{:b}}) = Diagonal(kron([1, -1], ones(Int64, dimension(hilbert, k)÷2)))
@inline commutator(k::TBAKind{:BdG}, hilbert::Hilbert{<:Phonon}) = Hermitian(kron([0 -1im; 1im 0], Diagonal(ones(Int, dimension(hilbert, k)÷2))))
@inline dimension(hilbert::Hilbert, ::TBAKind{:BdG}) = sum(dimension, values(hilbert))

"""
    AbstractTBA{K, H<:AbstractGenerator, G<:Union{Nothing, AbstractMatrix}} <: Engine

Abstract type for free quantum lattice systems using the tight-binding approximation.
"""
abstract type AbstractTBA{K, H<:AbstractGenerator, G<:Union{Nothing, AbstractMatrix}} <: Engine end
@inline contentnames(::Type{<:AbstractTBA}) = (:H, :commutator)
@inline kind(tba::AbstractTBA) = kind(typeof(tba))
@inline kind(::Type{<:AbstractTBA{K}}) where K = K
@inline Base.valtype(::Type{<:AbstractTBA{K, H} where K}) where {H<:AbstractGenerator} = valtype(eltype(H))
@inline dimension(tba::AbstractTBA{K, <:CompositeGenerator}) where K = length(getcontent(getcontent(tba, :H), :table))
@inline update!(tba::AbstractTBA; k=nothing, kwargs...) = ((length(kwargs)>0 && update!(getcontent(tba, :H); kwargs...)); tba)
@inline Parameters(tba::AbstractTBA) = Parameters(getcontent(tba, :H))

"""
    TBAMatrix{T, H<:AbstractMatrix{T}, G<:Union{AbstractMatrix, Nothing}} <: AbstractMatrix{T}

Matrix representation of a free quantum lattice system using the tight-binding approximation.
"""
struct TBAMatrix{T, H<:AbstractMatrix{T}, G<:Union{AbstractMatrix, Nothing}} <: AbstractMatrix{T}
    H::H
    commutator::G
    function TBAMatrix(H::AbstractMatrix, commutator::Union{AbstractMatrix, Nothing})
        new{eltype(H), typeof(H), typeof(commutator)}(H, commutator)
    end
end
@inline Base.size(m::TBAMatrix) = size(m.H)
@inline Base.getindex(m::TBAMatrix, i::Integer, j::Integer) = m.H[i, j]
@inline ishermitian(m::TBAMatrix) = ishermitian(typeof(m))
@inline ishermitian(::Type{<:TBAMatrix}) = true

"""
    TBAMatrixRepresentation{K<:AbstractTBA, V, T} <: MatrixRepresentation

Matrix representation of the Hamiltonian of a tight-binding system.
"""
struct TBAMatrixRepresentation{K<:AbstractTBA, V, T} <: MatrixRepresentation
    k::V
    table::T
    gauge::Symbol
    function TBAMatrixRepresentation{K}(k, table, gauge::Symbol=:icoord) where {K<:AbstractTBA}
        @assert gauge∈(:rcoord, :icoord) "TBAMatrixRepresentation error: gauge must be :rcoord or :icoord."
        return new{K, typeof(k), typeof(table)}(k, table, gauge)
    end
end
@inline TBAMatrixRepresentation{K}(table, gauge::Symbol=:icoord) where {K<:AbstractTBA} = TBAMatrixRepresentation{K}(nothing, table, gauge)
@inline Base.valtype(mr::TBAMatrixRepresentation) = valtype(typeof(mr))
@inline Base.valtype(::Type{<:TBAMatrixRepresentation{K}}) where {K<:AbstractTBA} = Matrix{promote_type(valtype(K), Complex{Int})}
@inline Base.valtype(::Type{<:TBAMatrixRepresentation{K, Nothing}}) where {K<:AbstractTBA} = Matrix{valtype(K)}
@inline Base.valtype(R::Type{<:TBAMatrixRepresentation}, ::Type{<:Union{Operator, Operators}}) = valtype(R)
@inline Base.zero(mr::TBAMatrixRepresentation) = zeros(eltype(valtype(mr)), length(mr.table), length(mr.table))
@inline Base.zero(mr::TBAMatrixRepresentation, ::Union{Operator, Operators}) = zero(mr)
@inline (mr::TBAMatrixRepresentation)(m::Operator; kwargs...) = add!(zero(mr, m), mr, m; kwargs...)

"""
    add!(dest::Matrix,
        mr::TBAMatrixRepresentation{<:AbstractTBA{TBAKind(:TBA)}},
        m::Operator;
        kwargs...
        ) -> typeof(dest)
    add!(dest::Matrix,
        mr::TBAMatrixRepresentation{<:AbstractTBA{TBAKind(:BdG)}},
        m::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID}}}};
        atol=atol/5,
        kwargs...
        ) -> typeof(dest)
    add!(dest::Matrix,
        mr::TBAMatrixRepresentation{<:AbstractTBA{TBAKind(:BdG)}},
        m::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:NID}}}};
        atol=atol/5,
        kwargs...
        ) -> typeof(dest)

Get the matrix representation of an operator and add it to destination.
"""
function add!(dest::Matrix,
        mr::TBAMatrixRepresentation{<:AbstractTBA{TBAKind(:TBA)}},
        m::Operator;
        kwargs...
        )
    seq₁, seq₂ = mr.table[m[1].index'], mr.table[m[2].index]
    coord = mr.gauge==:rcoord ? rcoord(m) : icoord(m)
    phase = isnothing(mr.k) ? one(eltype(dest)) : convert(eltype(dest), exp(-1im*dot(mr.k, coord)))
    dest[seq₁, seq₂] += m.value*phase
    return dest
end
function add!(dest::Matrix,
        mr::TBAMatrixRepresentation{<:AbstractTBA{TBAKind(:BdG)}},
        m::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID}}}};
        atol=atol/5,
        kwargs...
        )
    seq₁, seq₂ = mr.table[m[1].index'], mr.table[m[2].index]
    coord = mr.gauge==:rcoord ? rcoord(m) : icoord(m)
    phase = isnothing(mr.k) ? one(eltype(dest)) : convert(eltype(dest), exp(-1im*dot(mr.k, coord)))
    seq₁==seq₂ || (atol = 0)
    dest[seq₁, seq₂] += m.value*phase+atol
    if m[1].index.iid.nambu==creation && m[2].index.iid.nambu==annihilation
        seq₁, seq₂ = mr.table[m[1].index], mr.table[m[2].index']
        sign = statistics(eltype(m))==:f ? -1 : +1
        dest[seq₁, seq₂] += sign*m.value*phase'+atol
    end
    return dest
end
function add!(dest::Matrix,
        mr::TBAMatrixRepresentation{<:AbstractTBA{TBAKind(:BdG)}},
        m::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:NID}}}};
        atol=atol/5,
        kwargs...
        )
    if m[1] == m[2]
        seq = mr.table[m[1].index]
        dest[seq, seq] += 2*m.value+atol
    else
        seq₁, seq₂ = mr.table[m[1].index], mr.table[m[2].index]
        coord = mr.gauge==:rcoord ? rcoord(m) : icoord(m)
        phase = isnothing(mr.k) ? one(eltype(dest)) : convert(eltype(dest), exp(-1im*dot(mr.k, coord)))
        dest[seq₁, seq₂] += m.value*phase
        dest[seq₂, seq₁] += m.value'*phase'
    end
    return dest
end

"""
    TBAMatrixRepresentation(tba::AbstractTBA, k=nothing, gauge::Symbol=:icoord)

Construct the matrix representation transformation of a free quantum lattice system using the tight-binding approximation.
"""
@inline function TBAMatrixRepresentation(tba::AbstractTBA, k=nothing, gauge::Symbol=:icoord)
    table = getcontent(getcontent(tba, :H), :table)
    return TBAMatrixRepresentation{typeof(tba)}(k, table, gauge)
end

"""
    matrix(tba::Union{AbstractTBA, Algorithm{<:AbstractTBA}}; k=nothing, gauge=:icoord, atol=atol/5, kwargs...) -> TBAMatrix

Get the matrix representation of a free quantum lattice system.
"""
@inline function matrix(tba::AbstractTBA; k=nothing, gauge=:icoord, atol=atol/5, kwargs...)
    H = getcontent(tba, :H)
    commutator = getcontent(tba, :commutator)
    isnothing(commutator) && (atol = 0)
    return TBAMatrix(Hermitian(TBAMatrixRepresentation(tba, k, gauge)(expand(H); atol=atol, kwargs...)), commutator)
end
@inline function matrix(tba::AbstractTBA{TBAKind(:Analytical)}; kwargs...)
    return TBAMatrix(Hermitian(getcontent(tba, :H)(; kwargs...)), getcontent(tba, :commutator))
end
@inline matrix(tba::Algorithm{<:AbstractTBA}; kwargs...) = matrix(tba.engine; kwargs...)

"""
    eigen(m::TBAMatrix) -> Eigen

Solve the eigen problem of a free quantum lattice system.
"""
@inline eigen(m::TBAMatrix{T, H, Nothing}) where {T, H<:AbstractMatrix{T}} = eigen(m.H)
function eigen(m::TBAMatrix{T, H, G}) where {T, H<:AbstractMatrix{T}, G<:AbstractMatrix}
    W = cholesky(m.H)
    K = eigen(Hermitian(W.U*m.commutator*W.L))
    @assert length(K.values)%2==0 "eigen error: wrong dimension of matrix."
    for i = 1:(length(K.values)÷2)
        K.values[i] = -K.values[i]
    end
    V = inv(W.U)*K.vectors*sqrt(Diagonal(K.values))
    return Eigen(K.values, V)
end

"""
    eigvals(m::TBAMatrix) -> Vector

Get the eigen values of a free quantum lattice system.
"""
@inline eigvals(m::TBAMatrix{T, H, Nothing}) where {T, H<:AbstractMatrix{T}} = eigvals(m.H)
function eigvals(m::TBAMatrix{T, H, G}) where {T, H<:AbstractMatrix{T}, G<:AbstractMatrix}
    values = eigen(m.H*m.commutator).values
    @assert length(values)%2==0 "eigvals error: wrong dimension of matrix."
    for i = 1:(length(values)÷2)
        values[i] = -values[i]
    end
    return values
end

"""
    TBA{K, L<:AbstractLattice, H<:AbstractGenerator, G<:Union{AbstractMatrix, Nothing}} <: AbstractTBA{K, H, G}

The usual tight binding approximation for quantum lattice systems.
"""
struct TBA{K, L<:AbstractLattice, H<:AbstractGenerator, G<:Union{AbstractMatrix, Nothing}} <: AbstractTBA{K, H, G}
    lattice::L
    H::H
    commutator::G
    function TBA{K}(lattice::AbstractLattice, H::AbstractGenerator, commutator::Union{AbstractMatrix, Nothing}) where K
        @assert isa(K, TBAKind) "TBA error: wrong kind."
        if !isnothing(commutator)
            values = eigvals(commutator)
            num₁ = count(isapprox(+1, atol=atol, rtol=rtol), values)
            num₂ = count(isapprox(-1, atol=atol, rtol=rtol), values)
            @assert num₁==num₂==length(values)//2 "TBA error: unsupported input commutator."
        end
        new{K, typeof(lattice), typeof(H), typeof(commutator)}(lattice, H, commutator)
    end
end
@inline contentnames(::Type{<:TBA}) = (:lattice, :H, :commutator)

"""
    TBA(lattice::AbstractLattice, hilbert::Hilbert, terms::Tuple{Vararg{Term}}; boundary::Boundary=plain)

Construct a tight-binding quantum lattice system.
"""
@inline function TBA(lattice::AbstractLattice, hilbert::Hilbert, terms::Tuple{Vararg{Term}}; boundary::Boundary=plain)
    tbakind = TBAKind(typeof(terms))
    table = Table(hilbert, Metric(tbakind, hilbert))
    commt = commutator(tbakind, hilbert)
    return TBA{tbakind}(lattice, Generator(terms, Bonds(lattice), hilbert; half=false, table=table, boundary=boundary), commt)
end

"""
    TBA(lattice::AbstractLattice, hamiltonian::Function, parameters::Parameters, commt::Union{AbstractMatrix, Nothing}=nothing)

Construct a tight-binding quantum lattice system by providing the analytical expressions of the Hamiltonian.
"""
@inline function TBA(lattice::AbstractLattice, hamiltonian::Function, parameters::Parameters, commt::Union{AbstractMatrix, Nothing}=nothing)
    return TBA{TBAKind(:Analytical)}(lattice, Formulation(hamiltonian, parameters), commt)
end

"""
    EnergyBands{P, L<:Union{Colon, Vector{Int}}} <: Action

Energy bands by tight-binding-approximation for quantum lattice systems.
"""
struct EnergyBands{P, L<:Union{Colon, Vector{Int}}} <: Action
    path::P
    levels::L
    gauge::Symbol
    function EnergyBands(path, levels::Union{Colon, Vector{Int}}=Colon(), gauge::Symbol=:icoord)
        new{typeof(path), typeof(levels)}(path, levels, gauge)
    end
end
@inline prepare!(eb::EnergyBands{P, Colon}, tba::AbstractTBA) where P = (zeros(Float64, length(eb.path)), zeros(Float64, length(eb.path), dimension(tba)))
@inline prepare!(eb::EnergyBands{P, Vector{Int}}, tba::AbstractTBA) where P = (zeros(Float64, length(eb.path)), zeros(Float64, length(eb.path), length(eb.levels)))
@inline Base.nameof(tba::Algorithm{<:AbstractTBA}, eb::Assignment{<:EnergyBands}) = @sprintf "%s_%s" repr(tba, ∉(keys(eb.action.path))) eb.id
function run!(tba::Algorithm{<:AbstractTBA}, eb::Assignment{<:EnergyBands})
    for (i, params) in enumerate(pairs(eb.action.path))
        eb.data[1][i] = length(params)==1 && isa(first(params), Number) ? first(params) : i-1
        update!(tba; gauge=eb.action.gauge, params...)
        @timeit tba.timer "matrix" (m = matrix(tba.engine; params...))
        @timeit tba.timer "eigen" (eb.data[2][i, :] = eigen(m).values[eb.action.levels])
    end
end

"""
    BerryCurvature{B<:BrillouinZone} <: Action

Berry curvature of energy bands with the spirit of a momentum space discretization method by [Fukui et al, JPSJ 74, 1674 (2005)](https://journals.jps.jp/doi/10.1143/JPSJ.74.1674).
"""
struct BerryCurvature{B<:BrillouinZone} <: Action
    brillouinzone::B
    levels::Vector{Int}
end
@inline function prepare!(bc::BerryCurvature, tba::AbstractTBA)
    N₁, N₂ = periods(eltype(bc.brillouinzone))
    x = collect(Float64, 0:(N₁-1))/N₁
    y = collect(Float64, 0:(N₂-1))/N₂
    z = zeros(Float64, length(y), length(x), length(bc.levels))
    n = zeros(Float64, length(bc.levels))
    return (x, y, z, n)
end
function run!(tba::Algorithm{<:AbstractTBA}, bc::Assignment{<:BerryCurvature})
    N₁, N₂ = length(bc.data[1]), length(bc.data[2])
    eigenvectors = zeros(ComplexF64, N₁, N₂, dimension(tba.engine), length(bc.action.levels))
    for momentum in bc.action.brillouinzone
        coord = expand(momentum, bc.action.brillouinzone.reciprocals)
        @timeit tba.timer "matrix" (m = matrix(tba.engine; k=coord))
        @timeit tba.timer "eigen" (eigenvectors[Int(momentum[1])+1, Int(momentum[2])+1, :, :] = eigen(m).vectors[:, bc.action.levels])
    end
    g = isnothing(tba.engine.commutator) ? Diagonal(ones(Int, dimension(tba.engine))) : inv(tba.engine.commutator)
    @timeit tba.timer "Berry curvature" for momentum in bc.action.brillouinzone
        i₁, j₁ = Int(momentum[1]), Int(momentum[2])
        i₂, j₂ = (i₁+1)%N₁, (j₁+1)%N₂
        vs₁ = eigenvectors[i₁+1, j₁+1, :, :]
        vs₂ = eigenvectors[i₂+1, j₁+1, :, :]
        vs₃ = eigenvectors[i₂+1, j₂+1, :, :]
        vs₄ = eigenvectors[i₁+1, j₂+1, :, :]
        for k = 1:length(bc.action.levels)
            p₁ = vs₁[:, k]'*g*vs₂[:, k]
            p₂ = vs₂[:, k]'*g*vs₃[:, k]
            p₃ = vs₃[:, k]'*g*vs₄[:, k]
            p₄ = vs₄[:, k]'*g*vs₁[:, k]
            bc.data[3][j₁+1, i₁+1, k] = angle(p₁*p₂*p₃*p₄)
            bc.data[4][k] += bc.data[3][j₁+1, i₁+1, k]/2pi
        end
    end
    @info (@sprintf "Chern numbers: %s" join([@sprintf "%s(%s)" cn level for (cn, level) in zip(bc.data[4], bc.action.levels)], ", "))
end
@recipe function plot(pack::Tuple{Algorithm{<:AbstractTBA}, Assignment{<:BerryCurvature}})
    levels, chernnumbers = pack[2].action.levels, pack[2].data[4]
    layout --> length(levels)
    aspect_ratio --> :equal
    framestyle --> :none
    colorbar --> :bottom
    plot_title --> nameof(pack[1], pack[2])
    plot_titlefontsize --> 10
    title := @sprintf "level %s (C = %s)" first(levels) decimaltostr(first(chernnumbers))
    titlefontsize --> 8
    for (i, (level, chn)) in enumerate(zip(levels, chernnumbers))
        @series begin
            seriestype := :heatmap
            title := @sprintf "level %s (C = %s)" level decimaltostr(chn)
            subplot := i
            pack[2].data[1], pack[2].data[2], pack[2].data[3][:, :, i]
        end
    end
    primary := false
    ()
end

"""
    SampleNode(reciprocals::AbstractVector{<:AbstractVector}, position::Vector, levels::Vector{Int}, values::Vector, ratio::Number)
    SampleNode(reciprocals::AbstractVector{<:AbstractVector}, position::Vector, levels::Vector{Int}, values::Vector, ratios::Vector=ones(length(levels)))

A sample node of a momentum-eigenvalues pair.
"""
struct SampleNode
    k::Vector{Float64}
    levels::Vector{Int64}
    values::Vector{Float64}
    ratios::Vector{Float64}
end
@inline function SampleNode(reciprocals::AbstractVector{<:AbstractVector}, position::Vector, levels::Vector{Int}, values::Vector, ratio::Number)
    return SampleNode(reciprocals, position, levels, values, ones(length(levels))*ratio)
end
function SampleNode(reciprocals::AbstractVector{<:AbstractVector}, position::Vector, levels::Vector{Int}, values::Vector, ratios::Vector=ones(length(levels)))
    @assert length(reciprocals)==length(position) "SampleNode error: mismatched reciprocals and position."
    @assert length(levels)==length(values)==length(ratios) "SampleNode error: mismatched levels, values and ratios."
    return SampleNode(mapreduce(*, +, reciprocals, position), levels, values, ratios)
end

"""
    deviation(tba::Union{AbstractTBA, Algorithm{<:AbstractTBA}}, samplenode::SampleNode) -> Float64
    deviation(tba::Union{AbstractTBA, Algorithm{<:AbstractTBA}}, samplesets::Vector{SampleNode}) -> Float64

Get the deviation of the eigenvalues between the sample points and model points.
"""
function deviation(tba::Union{AbstractTBA, Algorithm{<:AbstractTBA}}, samplenode::SampleNode)
    diff = eigvals(matrix(tba; k=samplenode.k))[samplenode.levels] .- samplenode.values
    return real(sum(conj(diff) .* diff .* samplenode.ratios))
end
function deviation(tba::Union{AbstractTBA, Algorithm{<:AbstractTBA}}, samplesets::Vector{SampleNode})
    result = 0.0
    for samplenode in samplesets
        result += deviation(tba, samplenode)
    end
    return result
end

"""
    optimize!(tba::Union{AbstractTBA, Algorithm{<:AbstractTBA}},
        samplesets::Vector{SampleNode},
        variables=keys(Parameters(tba));
        verbose=false,
        method=LBFGS()
        ) -> Tuple{typeof(tba), Optim.MultivariateOptimizationResults}

Optimize the parameters of a tight binding system whose names are specified by `variables` so that the total deviations of the eigenvalues between the model points and sample points minimize.
"""
function optimize!(tba::Union{AbstractTBA, Algorithm{<:AbstractTBA}}, samplesets::Vector{SampleNode}, variables=keys(Parameters(tba)); verbose=false, method=LBFGS())
    v₀ = collect(getfield(Parameters(tba), name) for name in variables)
    function diff(v::Vector)
        parameters = Parameters{variables}(v...)
        update!(tba; parameters...)
        verbose && println(parameters)
        return deviation(tba, samplesets)
    end
    op = optimize(diff, v₀, method)
    parameters = Parameters{variables}(op.minimizer...)
    update!(tba; parameters...)
    return tba, op
end

end # module


text:
export eval_internals, set_value!, test_expr, create_var

function eval_internals(x::AffineExpr, vars::Tuple)
  (rows, cols) = x.size
  value = zeros(rows, cols)
  for var in vars
    coeff = x.vars_to_coeffs_map[object_id(var)]
    value += reshape(coeff.value * vec(var.value), rows, cols)
  end
  value += reshape(x.constant.value, rows, cols)
  return value
end

function eval_internals(x::SumSquaresExpr)
  value = 0
  for i in 1:length(x.affines)
    value += x.multipliers[i]*vecnorm(evaluate(x.affines[i]))^2
  end
  value += x.scalar
  return value
end

function set_value!(x::AffineExpr, value::Number)
  x.value = convert(Float64, value)
end

function set_value!(x::AffineExpr, value::AbstractArray)
  x.value = convert(SparseMatrixCSC{Float64, Int64}, sparse(value))
end

function test_expr(x::SumSquaresExpr)
  val1 = evaluate(x)
  val2 = eval_internals(x)
  @assert abs(val1 .- val2) <= 0.0001
end

function test_expr(x::AffineExpr, vars::Tuple)
  val1 = eval_internals(x, vars)
  val2 = evaluate(x)
  @assert all(abs(val1 .- val2) .<= 0.0001)
end

function create_var(m::Int64, n::Int64)
  x = Variable(m, n)
  set_value!(x, rand(m, n))
  return x
end


text:
# The purpose of this script is to compare the filtering capability of the single snapshots
# estimates and the gas filter on a given

#function that Generates the dynamical fitnesses and estimates single snapshot
using StaticNets,JLD, DynNets
 function simWrap1(N,T, NTV;dynType = "SIN", degIOUncMeans =15.*ones(2N),degb = [10, 40] )

    dynFitDgp,indsTVnodes =  StaticNets.dgpDynamic(StaticNets.fooErgmDirBin1,"SIN",N,T;
                                        NTV = NTV,degIOUncMeans = degIOUncMeans,degb = degb  )
    dynExpDegsDgp = zero(dynFitDgp)
    for t =1:T  dynExpDegsDgp[:,t] = StaticNets.expValStats(StaticNets.fooErgmDirBin1,dynFitDgp[:,t]) end
    # sanmple dynamical networks from dgp fitnesses
    dynDegsSam = zeros(dynFitDgp)
    for t=1:T
        matSamp = StaticNets.sampl(StaticNets.ErgmDirBin1(dynFitDgp[:,t]),1;  parGroupsIO = dynFitDgp[:,t])
        dynDegsSam[:,t] = [sumSq(matSamp,2);sumSq(matSamp,1)]
    end
    estFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(dynDegsSam); identPost = false,identIter= true )
    return dynFitDgp,dynDegsSam,estFitSS,indsTVnodes
 end

## Simulate estimate and save
Nvals = [ 2000,3000] #[50,60,80,90,100,125,150,200,300,400]# 200##40 #
 minDeg,maxDeg = (10, 40)
 Nsample = 3
 T = 250
    dynFitDgp = Array{Array{Float64,2},3}(length(Nvals),Nsample,2)
    dynDegsSam = Array{Array{Float64,2},3}(length(Nvals),Nsample,2)
    estFitSS = Array{Array{Float64,2},3}(length(Nvals),Nsample,2)
    filFitGas = Array{Array{Float64,2},3}(length(Nvals),Nsample,2)
    indsTVnodes = Array{BitArray{1},3}(length(Nvals),Nsample,2)
    rmseSSandGas = Array{Array{Float64,3},2}(length(Nvals),2)
    storeGasPar = Array{Array{Array{Float64,1},1},3}(length(Nvals),Nsample,2)
    convGasFlag = BitArray{3}(length(Nvals),Nsample,2)
   valDens = [false true]
   for d=1:2
     for indN = 1:length(Nvals)
        N=Nvals[indN]
        NTV = round.(Int,N/2)
        dense = valDens[d] #false#false#true
        if  dense
            unifDeg = round(Int,N*0.1) #15
            degb = [round(Int,N*0.1), round(Int,N*0.8)]
        else
            unifDeg = round(Int,(maxDeg-  minDeg)/2)
             degb = [round(Int,minDeg), round(Int,maxDeg)]#
        end
        tmprmseSSandGas = zeros(2N,2,Nsample)
        prog=0
        @time for s=1:Nsample
            println((d,N,s))
            #sample the dgp
             dynFitDgp[indN,s,d],dynDegsSam[indN,s,d],estFitSS[indN,s,d],indsTVnodes[indN,s,d] =
                   simWrap1(N,T,NTV;degIOUncMeans =unifDeg*ones(2N),degb = degb)
             tmprmseSSandGas[:,1,s] =sqrt.(meanSq((dynFitDgp[indN,s,d] - estFitSS[indN,s,d]).^2,2))
             #estimate single snapshot sequence
              round(s/Nsample,1)>prog? (prog = round(s/Nsample,1);println(prog)):()
              #estimate score_driven_filter_or_dgp
              ## e ora stima il gasPar
              degsIO_T = dynDegsSam[indN,s,d]
              modGasDirBin1 = DynNets.SdErgmDirBin1(degsIO_T,"FISHER-DIAG")
               parGas,convGasFlag[indN,s,d]  =  DynNets.estimateTarg(modGasDirBin1;SSest = estFitSS[indN,s,d])
               #parGas = [[0.989525],[0.0680809]]
              storeGasPar[indN,s,d] = parGas
              display(parGas[2:3])
              gasFiltFit,~ = DynNets.score_driven_filter_or_dgp(modGasDirBin1,[parGas[1];parGas[2];parGas[3]])
              tmprmseSSandGas[:,2,s] =sqrt.(meanSq((dynFitDgp[indN,s,d] - gasFiltFit).^2,2))
              filFitGas[indN,s,d] = gasFiltFit
              end
              if false # some descriptive plots
                  s=1
               fig, ax = subplots(3,2)
               ax[1,1][:set_title]("N = $(N) const Deg = $(unifDeg)  deg bounds = $(degb)  ")
               ax[1,1][:plot](dynFitDgp[indN,s,d][indsTVnodes[indN,s,d],:]')
               ax[1,2][:plot](dynFitDgp[indN,s,d][.!indsTVnodes[indN,s,d],:]')
               ax[2,1][:plot](estFitSS[indN,s,d][indsTVnodes[indN,s,d],:]',".")
               ax[3,1][:plot](dynFitDgp[indN,s,d][:,:]',estFitSS[indN,s,d][indsTVnodes[indN,s,d],:]',".")
               ax[2,2][:plot](dynDegsSam[indN,s,d][indsTVnodes[indN,s,d],:]')
               ax[3,2][:plot](dynDegsSam[indN,s,d][indmax(maximum(dynDegsSam[indN,s,d],2)),:])
           end
        rmseSSandGas[indN,d] = tmprmseSSandGas
     end
 end

 using JLD2
 save_fold = "./data/estimatesTest/asympTest/"
 @save(save_fold*"FilterTestDynamicIncreasingSize_extraLarge_2_denseAndSparse_$(Nsample).jld",
              dynFitDgp,dynDegsSam,estFitSS,indsTVnodes,storeGasPark,rmseSSandGas,filFitGas,convGasFlag)

##

save_fold = "./data/estimatesTest/asympTest/"
@load(save_fold*"FilterTestDynamicIncreasingSize_extraLarge_2_denseAndSparse_$(Nsample).jld",
             dynFitDgp,dynDegsSam,estFitSS,indsTVnodes,storeGasPark,rmseSSandGas,filFitGas,convGasFlag)

## Test single path sampling and FilterTestDynamicIncreasingSize

T=300
    ## Simulate estimate and save
    minDeg,maxDeg = (10, 160)
    Nsample = 1
    Nvals =  200#[200,300,400]#40 #[40,60,80,90,100,125,150]#
    d=1
    indN = 1
    N= 200
    NTV = round.(Int,N/2)
    dense = false#false#true
    if  dense
        unifDeg = round(Int,N*0.1) #15
        degb = [round(Int,N*0.1), round(Int,N*0.8)]
    else
        unifDeg = round(Int,(maxDeg-  minDeg)/2)
         degb = [round(Int,minDeg), round(Int,maxDeg)]#
    end
    tmprmseSSandGas = zeros(2N,2,Nsample)
    prog=0

    #sample the dgp
    dynFitDgp_test,dynDegsSam_test,estFitSS_test,indsTVnodes_test  =
       simWrap1(N,T,NTV;degIOUncMeans =unifDeg*ones(2N),degb = degb)
    tmprmseSSandGas_test =sqrt.(meanSq((dynFitDgp_test - estFitSS_test).^2,2))
    #estimate single snapshot sequence
    #estimate score_driven_filter_or_dgp
    ## e ora stima il gasPar
    degsIO_T = dynDegsSam_test

 modGasDirBin1 = DynNets.SdErgmDirBin1(degsIO_T,"")#"")###
 parGas,convGasFlag_test  =  DynNets.estimateTarg(modGasDirBin1;SSest = estFitSS_test)

 modGasDirBin1_2 = DynNets.SdErgmDirBin1(degsIO_T,"FISHER-DIAG")#"")###
  parGas2,convGasFlag_test  =  DynNets.estimateTarg(modGasDirBin1_2;SSest = estFitSS_test)
 #gasFiltFit,~ = DynNets.score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[parGas[1];parGas[2];0.07])
 println((parGas[2],parGas[3]))
println((parGas2[2],parGas2[3]))
gasFiltFit,~ = DynNets.score_driven_filter_or_dgp( modGasDirBin1,[parGas[1];parGas[2];parGas[3]])
    tmprmseSSandGas_test =sqrt.(meanSq((dynFitDgp_test - gasFiltFit).^2,2))
    filFitGas_test = gasFiltFit
 indsNodesPlot = [1,N]
 #plot(degsIO_T[indsNodesPlot,:]')
 #close("all")
    fig = figure("pyplot_subplot_mixed",figsize=(10,10)) # Create a new blank figure
    subplot(221) # Create the 1st axis of a 2x2 arrax of axes
    title("min max degs =  $((minimum(degsIO_T[indsNodesPlot[1],:]), maximum(degsIO_T[indsNodesPlot[1],:])))") # Give the most recent axis a title
    plot(1:T,dynFitDgp_test[indsNodesPlot[1],:],"--",
       1:T,estFitSS_test[indsNodesPlot[1],:],".")
    plot(1:T,gasFiltFit[indsNodesPlot[1],:],"-")
    suptitle("Compare Fiter Single path 2 nodes")
    subplot(223)
    title("min max degs =  $((minimum(degsIO_T[indsNodesPlot[2],:]), maximum(degsIO_T[indsNodesPlot[2],:])))") # Give the most recent axis a title
    plot(1:T,dynFitDgp_test[indsNodesPlot[2],:],"--",
       1:T,estFitSS_test[indsNodesPlot[2],:],".")
    plot(1:T,gasFiltFit[indsNodesPlot[2],:],"-")

 #gasFiltFit,~ = DynNets.score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[parGas2[1];parGas2[2];0.07])
 gasFiltFit,~ = DynNets.score_driven_filter_or_dgp( modGasDirBin1_2,[parGas2[1];parGas2[2];parGas2[3]])
    tmprmseSSandGas_test =sqrt.(meanSq((dynFitDgp_test - gasFiltFit).^2,2))
    filFitGas_test = gasFiltFit
 indsNodesPlot = [1,N]
 #plot(degsIO_T[indsNodesPlot,:]')
 #close("all")
    fig = figure("pyplot_subplot_mixed",figsize=(10,10)) # Create a new blank figure
    subplot(222) # Create the 1st axis of a 2x2 arrax of axes
    title("min max degs =  $((minimum(degsIO_T[indsNodesPlot[1],:]), maximum(degsIO_T[indsNodesPlot[1],:])))") # Give the most recent axis a title
    plot(1:T,dynFitDgp_test[indsNodesPlot[1],:],"--",
       1:T,estFitSS_test[indsNodesPlot[1],:],".")
    plot(1:T,gasFiltFit[indsNodesPlot[1],:],"-")
    suptitle("Compare Fiter Single path 2 nodes")
    subplot(224)
    title("min max degs =  $((minimum(degsIO_T[indsNodesPlot[2],:]), maximum(degsIO_T[indsNodesPlot[2],:])))") # Give the most recent axis a title
    plot(1:T,dynFitDgp_test[indsNodesPlot[2],:],"--",
       1:T,estFitSS_test[indsNodesPlot[2],:],".")
    plot(1:T,gasFiltFit[indsNodesPlot[2],:],"-")

plot(maximum(degsIO_T,2),".")
plot(mean(degsIO_T,2),".")
plot(minimum(degsIO_T[indsTVnodes_test,:],2))

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
using DrWatson
@quickactivate "EPExperiment"

using EPExperiments

using FileIO
using JLD2
using LaTeXStrings
using StatsPlots
using StatsFuns: logsumexp

function get_sir_losses(chn, cost_fn)
    betas = chn[:β].value.data
    return cost_fn.(betas, )
end

function analyse_nuts_baseline(
    subfolder,
    experiment_name;
    thinning=1,
    burn_in=1001,
    num_replications=1
)
    # Load chain and experiment config
    experiment_folder = datadir(subfolder, experiment_name)
    results = load(joinpath(experiment_folder, "results.jld2"))
    experiment_config = results["experiment_config"]
    chn = results["nuts_samples"]

    # Get cost function
    if subfolder == "sir_nuts_baseline"
        # Indicates the number of intermediate distributions in AnIS. This is
        # important for adjusting the x-axis in the ESS plot.
        num_intermediates = 100

        @unpack cost_fn_name = experiment_config
        cost_fn = COST_FUNCTIONS[cost_fn_name]

        @unpack data_file = experiment_config
        input_data = load(datadir("sims", data_file))
        true_gamma = input_data["params"].γ
        @unpack i_0_true, β_true = input_data["params"]

        betas = Array(chn[:β])[burn_in:end,:]
        losses = cost_fn.(betas, true_gamma)

        # Plot the samples from the joint
        i_0s = Array(chn[:i₀])[burn_in:end,:]
        p = scatter(
            Array(chn[:β])[1:(burn_in-1),:],
            Array(chn[:i₀])[1:(burn_in-1),:],
            xlabel=L"\beta",
            ylabel=L"i_0",
            label="Burn-in samples",
            alpha=0.5,
            color=:black,
            markersize=2.5
        )
        scatter!(
            p,
            betas,
            i_0s,
            xlabel=L"\beta",
            ylabel=L"i_0",
            label="HMC samples",
            alpha=0.5,
            legend=false,
            xlims=[0,2.3],
            xticks=[0.0, 0.5, 1.0, 1.5, 2.0],
            ylims=[0,60],
            color=palette(:default)[3],
            fontfamily="serif-roman",
            framestyle=:semi,
            grid=false,
            thickness_scaling=1.5,
            markersize=2.5
        )
        savefig(p, joinpath(experiment_folder, "sir_joint_dist_mcmc_with_burn_in.png"))

        # Plot the maximum weight
        # veclosses = vec(losses)
        # ixs = Array(0:4500:size(veclosses, 1))
        # ixs[1] = 90
        # cummax_weight = Array{Float64,1}(undef, length(ixs))
        # for (i, ix) in enumerate(ixs)
        #     cummax_weight[i] = maximum(veclosses[1:ix])
        # end

        # plot_ixs = Array(0:50:10_000)
        # plot_ixs[1] = 1
        # p = plot(
        #     plot_ixs,
        #     cummax_weight,
        #     xscale=:log10,
        #     xlabel="Number of Samples",
        #     ylabel="Maximum log weight",
        #     legend=false
        # )
        # savefig(p, joinpath(experiment_folder, "max_weight.png"))
    elseif subfolder == "radon_nuts_baseline"
        num_intermediates = 200
        @unpack cost_fn_name = experiment_config
        @show cost_fn_name
        cost_fn = RADON_COST_FNS[cost_fn_name]

        # Compute weights by applying cost function to all samples
        # If we just use Array(chn) then the different chains get appended and
        # we "lose" the third dimension. Using chn.value.data preserves the three
        # dimensional shape.
        alpha = chn[burn_in:end,namesingroup(chn, :alpha),1:2000].value.data
        beta = chn[burn_in:end,namesingroup(chn, :beta),1:2000].value.data
        losses = Array{Float64}(undef, size(alpha, 1), size(alpha, 3))
        fill!(losses, -Inf)
        for sample_ix in 1:size(losses, 1)
            for chain_ix in 1:size(losses, 2)
                losses[sample_ix,chain_ix] = cost_fn(
                    alpha[sample_ix,:,chain_ix],
                    beta[sample_ix,:,chain_ix],
                    POST_PRED_COUNTY_IDX,
                    POST_PRED_FLOOR
                )
            end
        end
    elseif subfolder == "posterior_predictive_mh_baseline"
        num_intermediates = 100

        cost_fn(x) = pdf(MvNormal(-POST_PRED_Y_OBSERVED, 0.5*I), x)

        savefig(
            plot(chn), joinpath(experiment_folder, "traces.png")
        )

        xs = chn[burn_in:end,namesingroup(chn, :x),:].value.data
        println("Cost function evaluation time:")
        @time begin
        losses = Array{Float64}(undef, size(xs, 1), size(xs, 3))
        for sample_ix in 1:size(losses, 1)
            for chain_ix in 1:size(losses, 2)
                losses[sample_ix,chain_ix] = cost_fn(
                    xs[sample_ix,:,chain_ix]
                )
            end
        end
        end
    end

    # ESS based on Turing implementation
    # fxs = losses
    # @show size(fxs)
    # ess_vals = Array{Float64}(undef, num_replications)
    # chains_per_replication = Int(size(fxs, 2) / num_replications)
    # start_ix = 1
    # for i in 1:num_replications
    #     @show start_ix
    #     fxs_replication = fxs[:,start_ix:(start_ix+chains_per_replication-1)]
    #     fxs_chn = Chains(
    #         reshape(
    #             fxs_replication,
    #             size(fxs_replication, 1),
    #             1,
    #             size(fxs_replication, 2)
    #         ),
    #         ["fx"]
    #     )
    #     ess_vals[i] = ess(fxs_chn).nt[:ess][1]
    #     start_ix = start_ix + chains_per_replication
    # end
    # @show ess_vals
    # return ess_vals

    # Compute ESS based on those weights
    # num_samples = prod(size(losses))
    fxs = reshape(losses, :, num_replications)
    # println("ESS time:")
    # @time begin
    # #fake_ess = ess_batch(fake_log_weights; thinning=thinning)
    # fake_ess = Array{Float64,2}(undef, 100, 10)
    # end

    # Compute the estimates of the actual expecation
    @show size(fxs)
    println("Expectation estimation time:")
    @time begin
    exp_estimates = cumsum(fxs, dims=1) ./ Array(1:size(fxs, 1))
    end

    @show size(fake_ess)
    @show size(exp_estimates)
    @tagsave(joinpath(experiment_folder, "ess.jld2"), Dict(
        "fake_ess" => fake_ess,
        "num_intermediates" => num_intermediates,
        "exp_estimates" => exp_estimates
    ))

    return
    ixs = Array(1:thinning:num_samples) ./ num_intermediates
    ixs[1] = 1
    ess_plot = plot(
        ixs,
        fake_ess,
        xscale=:log10,
        yscale=:identity,
        label="MCMC",
        legend=true,
        xlabel="Number of Samples",
        ylabel="ESS"
    )
    savefig(ess_plot, joinpath(experiment_folder, "ess_plot_burn_in=$(burn_in).png"))
end

# analyse_nuts_baseline(
#     "sir_nuts_baseline",
#     "combined_baselines_5replications";
#     thinning=90,
#     num_replications=5
# )

# analyse_nuts_baseline(
#     "posterior_predictive_mh_baseline",
#     "new_model_combined_baselines_5mh_10replications";
#     thinning=90,
#     num_replications=10
# )

analyse_nuts_baseline(
    "radon_nuts_baseline",
    "combined_baselines_10replications";
    thinning=90,
    num_replications=10
)


text:
using DiffusionMap
using Test

@testset "DiffusionMap.jl" begin
    # Write your own tests here.
end


text:
function fib(n)
#
#  test polynomial suggested by Goedecker
#
    [-1.0; ones(n)]
end


text:
using Plots

# TO solution
x̄, ū = unpack(z̄, prob)
@show tf_nom = sum([ū[t][end] for t = 1:T-1]) # total time
t_nominal = range(0, stop = tf_nom, length = T)
plot(t_nominal, hcat(x̄...)', width = 10.0, color = :lawngreen,
	grid = false, label = false, axes = false)
plot(t_nominal[1:end-1], hcat(ū...)[1:1, :]',
	linetype = :steppost, width = 2.0, color = goldenrod_color)

# DPO solution
x, u = unpack(z, prob_dpo.prob.prob.nom)
@show tf_dpo = sum([u[t][end] for t = 1:T-1])
t_dpo = range(0, stop = tf_dpo, length = T)

x_sample = []
u_sample = []
t_sample = []
plt = plot()
for i = 1:2 * model.n
	z_sample = z[prob_dpo.prob.idx.sample[i]]
	_x_sample, _u_sample = unpack(z_sample, prob_dpo.prob.prob.nom)
	tf_sample = sum([_u_sample[t][end] for t = 1:T-1])
	@show tf_sample
	push!(x_sample, _x_sample)
	push!(u_sample, _u_sample)
	push!(t_sample, range(0, stop = tf_sample, length = T))
	# plt = plot!(t_sample[end], hcat(x_sample[end]...)',
	# 	color = :black, width = 5.0, label = "", grid = false)
	plt = plot!(t_sample[end], hcat(u_sample[end]..., u_sample[end][end])[1:1,:]',
		color = :black, width = 5.0, label = "", grid = false, linetype = :steppost)
end
# plt = plot!(t_dpo, hcat(x...)', width = 10.0, color = :lawngreen, label = "", grid = false)

plot!(t_dpo, hcat(u..., u[end])[1:1, :]',
	linetype = :steppost, width = 10.0, color = :lawngreen, label = "", grid = false)

display(plt)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PolicyRulesWithSubjects prescribes a test that applies to a request to an apiserver. The test considers the subject making the request, the verb being requested, and the resource to be acted upon. This PolicyRulesWithSubjects matches a request if and only if both (a) at least one member of subjects matches the request and (b) at least one member of resourceRules or nonResourceRules matches the request.

    IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects(;
        nonResourceRules=nothing,
        resourceRules=nothing,
        subjects=nothing,
    )

    - nonResourceRules::Vector{IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule} : `nonResourceRules` is a list of NonResourcePolicyRules that identify matching requests according to their verb and the target non-resource URL.
    - resourceRules::Vector{IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule} : `resourceRules` is a slice of ResourcePolicyRules that identify matching requests according to their verb and the target resource. At least one of `resourceRules` and `nonResourceRules` has to be non-empty.
    - subjects::Vector{IoK8sApiFlowcontrolV1alpha1Subject} : subjects is the list of normal user, serviceaccount, or group that this rule cares about. There must be at least one member in this slice. A slice that includes both the system:authenticated and system:unauthenticated user groups matches every request. Required.
"""
mutable struct IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects <: SwaggerModel
    nonResourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule} } # spec name: nonResourceRules
    resourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule} } # spec name: resourceRules
    subjects::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1Subject} } # spec name: subjects

    function IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects(;nonResourceRules=nothing, resourceRules=nothing, subjects=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects, Symbol("nonResourceRules"), nonResourceRules)
        setfield!(o, Symbol("nonResourceRules"), nonResourceRules)
        validate_property(IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects, Symbol("resourceRules"), resourceRules)
        setfield!(o, Symbol("resourceRules"), resourceRules)
        validate_property(IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects, Symbol("subjects"), subjects)
        setfield!(o, Symbol("subjects"), subjects)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects

const _property_map_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects = Dict{Symbol,Symbol}(Symbol("nonResourceRules")=>Symbol("nonResourceRules"), Symbol("resourceRules")=>Symbol("resourceRules"), Symbol("subjects")=>Symbol("subjects"))
const _property_types_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects = Dict{Symbol,String}(Symbol("nonResourceRules")=>"Vector{IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule}", Symbol("resourceRules")=>"Vector{IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule}", Symbol("subjects")=>"Vector{IoK8sApiFlowcontrolV1alpha1Subject}")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects)
    (getproperty(o, Symbol("subjects")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects }, name::Symbol, val)
end