text:
"""
`OrthogonalPolynomials` is the abstract supertype of all univariate orthogonal
polynomials.
"""
abstract type OrthogonalPolynomials{T} <: PolynomialBasis{T}
end

const OPS{T} = OrthogonalPolynomials{T}

abstract type OrthogonalPolynomial{T} <: Polynomial{T} end

approx_length(b::OPS, n::Int) = n

# Taking the derivative of a polynomial lowers its degree. This results in a
# rectangular differentiation matrix. With the reducedegree option set to false,
# the degree is maintained (resulting in a square differentiation matrix).
function derivative_dict(dict::OPS, order::Int; reducedegree = true, options...)
	@assert order >= 0
	reducedegree ? resize(dict, length(dict)-order) : dict
end

function antiderivative_dict(s::OPS, order::Int; options...)
	@assert order >= 0
	resize(s, length(s)+order)
end

size(o::OrthogonalPolynomials) = (o.n,)

p0(::OPS{T}) where {T} = one(T)

hasextension(b::OPS) = true

extension(::Type{T}, src::O, dest::O; options...) where {T,O <: OPS} = IndexExtension{T}(src, dest, 1:length(src))
restriction(::Type{T}, src::O, dest::O; options...) where {T,O <: OPS} = IndexRestriction{T}(src, dest, 1:length(dest))

include("recurrence.jl")

# Default evaluation of an orthogonal polynomial: invoke the recurrence relation
unsafe_eval_element(b::OPS, idx::PolynomialDegree, x) =
    recurrence_eval(b, idx, x)
function unsafe_eval_element_derivative(b::OPS, idx::PolynomialDegree, x, order)
	@assert order == 1
    recurrence_eval_derivative(b, idx, x)
end

hasmeasure(dict::OPS) = true

measure(ϕ::OrthogonalPolynomial) = measure(dictionary(ϕ))

weightfun(b::OPS, x) = weightfun(measure(b), x)

function gram1(T, dict::OPS, m; options...)
    isorthonormal(dict, m) && return IdentityOperator{T}(dict)
	if isorthogonal(dict, m)
		diagonal_gram(T, dict, m; options...)
	else
		default_gram(T, dict, m; options...)
	end
end

diagonal_gram(::Type{T}, dict::OPS, measure; options...) where {T} =
    default_diagonal_gram(T, dict::OPS, measure; options...)

function symmetric_jacobi_matrix(b::OPS)
    T = codomaintype(b)
    n = length(b)
    α, β = monic_recurrence_coefficients(b)
    J = zeros(T, n, n)
    for k in 1:n
        J[k,k] = α[k]
        if k > 1
            J[k,k-1] = sqrt(β[k])
        end
        if k < n
            J[k,k+1] = sqrt(β[k+1])
        end
    end
    J
end

function roots(b::OPS{T}) where {T<:Number}
    J = symmetric_jacobi_matrix(b)
    eigen(J).values
end

function roots(b::OPS{T}) where {T<:Union{BigFloat}}
    J = symmetric_jacobi_matrix(b)
    # assuming the user has imported GenericLinearAlgebra.jl
    sort(real(eigvals!(J)))
end

interpolation_grid(b::OPS) = roots(b)
hasinterpolationgrid(dict::OPS) = true
opsorthogonal(dict, measure) = length(dict) -issymmetric(dict) <= length(points(measure))

"Return the first moment, i.e., the integral of the weight function."
function first_moment(b::OPS)
    @warn "implement first_moment for $b"
end

"""
Compute the Gaussian quadrature rule using the roots of the orthogonal polynomial.
"""
function gauss_rule(b::OPS{T}) where {T <: Real}
    J = symmetric_jacobi_matrix(b)
    x,v = eigen(J)
    b0 = first_moment(b)
    # In the real-valued case it is sufficient to use the first element of the
    # eigenvector. See e.g. Gautschi's book, "Orthogonal Polynomials and Computation".
    w = b0 * v[1,:].^2
    x,w
end

# In the complex-valued case, we have to compute the weights by summing explicitly
# over the full eigenvector.
function gauss_rule(b::OPS{T}) where {T <: Complex}
    J = symmetric_jacobi_matrix(b)
    x,v = eigen(J)
    b0 = first_moment(b)
    w = similar(x)
    for i in 1:length(w)
        w[i] = b0 * v[1,i]^2 / sum(v[:,i].^2)
    end
    x,w
end

import FastGaussQuadrature: gaussjacobi
using GaussQuadrature: jacobi
gaussjacobi(n::Integer,α::T,β::T) where T<:AbstractFloat =
    jacobi(n, α, β)

"""
Compute the weights of Gaussian quadrature rule from the given roots of the
orthogonal polynomial and using the formula:
```
w_j = 1 / \\sum_{k=0}^{n-1} p_k(x_j)^2
```
This formula only holds for the orthonormal polynomials.
"""
function gauss_weights_from_points(b::OPS, roots)
    pk = zeros(eltype(roots), length(b))
    w = zeros(eltype(roots), length(b))
    for j in 1:length(w)
        recurrence_eval_orthonormal!(pk, b, roots[j])
        w[j] = 1/sum(pk.^2)
    end
    w
end

function leading_order_coefficient(b::OPS{T}, idx) where {T}
    @assert 1 <= idx <= length(b)
    γ = one(T)
    for k in 0:idx-2
        γ *= rec_An(b, k)
    end
    γ
end

"""
Stieltjes Algorithm

Given nodes (t_k) and weights (w_k) of a quadrature rule, the Stieltjes Algorithm
calculates the N first recurrence coefficients of the monic polynomials (p_k) satisfying

```
p_{k+1}(t) = (t-a_k)p_k(t)-b_kp_{k-1}(t)
p_0(t) = 1, p_{-1}(t) = 0
```

and orthogonal with respect to the inner product induced by the quadrature

```
\\sum w_k p_m(t_k)q_n(t_k) = \\delta_{m,n}
```

See e.g. Gautschi's book, "Orthogonal Polynomials and Computation" and its
accompanying code `https://www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html`

Note: The vectors used have indices starting at `1` (such that `a_n` above is
given by `a[n+1]`). The last value is `a_{n-1} = a[n]`.
"""
function stieltjes(N, x::Array{RT},w::Array{T}) where {T,RT}
    @assert real(T) == RT
    # ####
    # Remove data with zero weights done in, e.g, Gautschi's book, "Orthogonal Polynomials and Computation".
    # ####
    # indices = sortperm(w)
    # x = x[indices]
    # w = w[indices]
    # index = minimum(find(!(w==0)))
    # w = w[index:length(w)]
    # x = x[index:length(w)]
    #
    # indices = sortperm(x)
    # x = x[indices]
    # w = w[indices]

    Ncap=length(x)
    # ####
    if N < 1 || N > Ncap
        error("Too little quadrature points")
    end

    a = zeros(T,N)
    b = zeros(T,N)

    p1 = zeros(T,Ncap)
    p2 = ones(T,Ncap)
    # The following scaling has to be adopted in case the
    # orthogonal polynomials underflow or overflow. In the
    # former case, c has to be chosen sufficiently large
    # positive, in the latter case sufficiently small
    # positive. (The example below relates to Table2_9.m.)
    #
    # if N==320
    #     c=1e50;
    #     p2=c*p2;
    #     s0=c^2*s0;
    # end
    p0 = similar(p1)
    scratch = similar(p1)
    stieltjes!(a,b,N,x,w,p0,p1,p2,scratch)
    a,b
end

"""
In-place Stieltjes Algorithm

see `stieltjes`

`p0`, `p1`, `p2`, and `scratch` are vectors of size `N`
"""
function stieltjes!(a::Array{T},b::Array{T},N,x::Array{RT},w::Array{T},p0::Array{T},p1::Array{T},p2::Array{T},scratch::Array{T}) where {T, RT}
    @assert real(T) == RT
    tiny = 10*typemin(real(T))
    huge = .1*typemax(real(T))

    s0 = sum(w)
    a[1] = dot(x,w)/s0
    b[1] = s0
    for k=1:N-1
        copyto!(p0,p1)
        copyto!(p1,p2)

        p2 .= (x .- a[k]) .* p1 .- b[k] .* p0
        scratch .= (p2 .^ 2) .* w
        s1 = sum(scratch)
        scratch .= w .* (p2 .^ 2) .* x
        s2 = sum(scratch)
        scratch .= abs.(p2)
        if (maximum(abs.(scratch))>huge || abs(s2)>huge)
            error("impending overflow in stieltjes for k=$(k)")
        end
        if abs(s1)<tiny
            error("impending underflow in stieltjes for k=$(k)")
        end
        a[k+1] = s2/s1
        b[k+1] = s1/s0
        s0 = s1
    end
end

"""
Modified Chebyshev Algorithm

Implements the map of the modified moments

```
\\mu_k = \\int p_k(x) dλ(x), k=0,...,2N-1
```
where (p_k) are some monic orthogonal polynomials satisfying the recurrence relation
```
p_{k+1}(t) = a_k(t_k)p_k(t) - b_kp{k-1}(t)
p_0(t) = 1, p_{-1}(t) = 0
```
to the first recurrence coefficients (α_k), (β_k), with k=0,...,N-1, such that
the monic polynomials (q_k) satisfying
```
q_{k+1}(t) = α_k(t_k)q_k(t) - β_kq{k-1}(t)
q_0(t) = 1, q_{-1}(t) = 0
```
are orthogonal with respect to the measure λ.

By default a, and b are zero vectors and therefore the moments are assumed to be
```
\\mu_k = \\int x^k dλ(x), k=0,...,2N-1
```
This procedure is however not stable.

See e.g. Gautschi's book, "Orthogonal Polynomials and Computation". Algorithm 2.1, p77

Note: The vectors used have indices starting at `1` (such that `α_n` above is
given by `α[n+1]`). The last value is `α_{n-1} = α[n]`.
"""
function modified_chebyshev(m::Array{T}, a=zeros(T,length(m)), b=zeros(T,length(m))) where {T}
    L = length(m)
    n = L>>1
    α = Array{T}(undef, n)
    β = Array{T}(undef, n)
    σmone = Array{T}(undef, L)
    σzero = Array{T}(undef, L)
    σ = Array{T}(undef, L)

    modified_chebyshev!(α,β,m,a,b,σ,σzero,σmone,n,1)
    α, β
end

"""
Modified Chebyshev Algorithm

see `chebyshev`

`σ`, `σzero`, and `σmone` are vectors of size or larger then `2n`
`n` is the number of coefficients required and `os` is the offset, i.e.,
the initial index that can be used in `σ`, `σzero`, and `σmone`.
"""
function modified_chebyshev!(α,β,m,a,b,σ,σzero,σmone,n=length(α),os=1)
    # Initialisation
    α[1] = a[1]+m[2]/m[1]
    β[1] = m[1]
    fill!(σmone,0)
    copyto!(σzero,os,m,os,2n)
    # continue
    for k=1:n-1
        for l in k+os:2n-k-1+os
            σ[l] = σzero[l+1]-(α[k]-a[l]).*σzero[l]-β[k].*σmone[l]+b[l].*σzero[l-1]
        end
        α[k+1] = a[k+1]
        α[k+1] += σ[k+2]/σ[k+1]
        α[k+1] = a[k+1]+σ[k+2]/σ[k+1]-σzero[k+1]/σzero[k]
        β[k+1] = σ[k+1]/σzero[k]
        copyto!(σmone,os,σzero,os,2n)
        copyto!(σzero,os,σ,os,2n)
    end
    nothing
end

"""
Calculates the first `n` recurrence coefficients (`α_k`) and (`β_k`) of the monic orthogonal polynomials
up to the tolerance `tol`
given a quadrature rule `my_quadrature_rule` with a relative degree of exactness
```
d(M)/M = δ + O(1/M) as M goes to infinity
```
with d(M) the degree of exactness. Thus δ is 2 for Gaussian quadrature.

See equation page 101 from Gautschi's book, "Orthogonal Polynomials and Computation"
"""
function adaptive_stieltjes(n,my_quadrature_rule::Function; tol = 1e-12, δ = 1, maxits = 20, quadrature_size=false)
    M = 1 + floor(Int, (2n-1)/δ)
    nodes, weights = my_quadrature_rule(M)
    ELT = eltype(nodes)

    α0,β0 = stieltjes(n,nodes,weights)

    M = M+1
    nodes, weights = my_quadrature_rule(M)
    α1,β1 = stieltjes(n,nodes,weights)

    no_its = 2
    while reduce(&, abs.(β1-β0) .> tol.*abs.(β1))
        if no_its > maxits
            @warn("accuracy of Stieltjes is not obtained, degree of Quadrature is $(M) and error is $(maximum(abs.(β1-β0)./abs.(β1)))")
            break
        end
        M = M+2^floor(Int,no_its/5)
        no_its = no_its + 1
        copyto!(α0,α1)
        copyto!(β0,β1)

        nodes, weights = my_quadrature_rule(M)
        α1,β1 = stieltjes(n,nodes,weights)
    end
    quadrature_size ? (α1, β1, M) : (α1, β1)
end

"""
Transforms the `N` recurrence coefficients (α_k) and (β_k) of monic
orthogonal polynomials (p_k), which satisfy the three-term recurrence relation
```
p_{k+1}(t) = (t-α_k)p_k(t)-β_kp_{k-1}(t)
```
to the `N-1` first recurrence coefficients (a_k), (b_k), and (c_k) the associated
monic polynomials (q_k), such that

```
q_{k+1}(t) = (a_kt+b_k)q_k(t)-c_kq_{k-1}(t)
```
"""
function monic_to_orthonormal_recurrence_coefficients(α::Array{T}, β::Array{T}) where {T}
    n = length(α)
    a = Array{T}(undef, n-1)
    b = Array{T}(undef, n-1)
    c = Array{T}(undef, n-1)
    monic_to_orthonormal_recurrence_coefficients!(a,b,c,α,β)
    a,b,c
end

"""
See `monic_to_orthonormal_recurrence_coefficients`
"""
function monic_to_orthonormal_recurrence_coefficients!(a::Array{T},b::Array{T},c::Array{T},α::Array{T},β::Array{T}) where {T}
    a .= 1 ./ sqrt.(view(β,2:length(β)))
    b .= -1.0.*view(α,1:length(α)-1)./sqrt.(view(β,2:length(β)))
    c .= sqrt.(view(β,1:length(β)-1)./view(β,2:length(β)))
    a,b,c
end


text:
module Brushes

using ..Ahorn, Maple

displayName = "Brushes"
group = "Brushes"

drawingLayers = Ahorn.Layer[]

toolsLayer = nothing
targetLayer = nothing
material = '0'

brushes = Ahorn.Brush[
    Ahorn.Brush(
        "Pencil",
        hcat(1)
    ),

    Ahorn.Brush(
        "Dither",
        [
            1 0;
            0 1
        ]
    ),
    Ahorn.Brush(
        "Ahorn",
        [
            0 0 1 0 0 0 0;
            0 1 0 1 0 0 0;
            0 0 0 1 1 0 1;
            0 0 1 1 1 0 0;
            0 1 1 1 1 1 0;
            0 1 1 1 1 1 0;
            1 1 1 1 1 1 0;
            1 1 1 1 1 0 0;
            0 1 1 1 0 0 0;
        ]
    )
]

selectedBrush = brushes[1]

hoveringBrush = nothing
phantomBrushes = Dict{Tuple{Integer, Integer}, Ahorn.Brush}()

function applyPhantomBrushes()
    if !isempty(phantomBrushes)
        Ahorn.History.addSnapshot!(Ahorn.History.RoomSnapshot("Brush ($(selectedBrush.name), $material)", Ahorn.loadedState.room))
    end

    for (pos, brush) in phantomBrushes
        x, y = pos

        Maple.updateTileSize!(Ahorn.loadedState.room, '0', '0')
        roomTiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
        Ahorn.applyBrush!(brush, roomTiles, material, x, y)
    end

    if !isempty(phantomBrushes)
        empty!(phantomBrushes)

        Ahorn.redrawLayer!(targetLayer)
        Ahorn.redrawLayer!(toolsLayer)
    end
end

function drawBrushes(layer::Ahorn.Layer, room::Maple.Room)
    if !isa(hoveringBrush, Nothing)
        x, y, brush = hoveringBrush

        Ahorn.drawBrush(brush, layer, x, y)
    end

    # Make sure we should actually render the phantom brushes
    if Ahorn.mouseButtonHeld(0x1)
        for (pos, brush) in phantomBrushes
            x, y = pos

            Ahorn.drawBrush(brush, layer, x, y)
        end
    end
end

function cleanup()
    global hoveringBrush = nothing
    empty!(phantomBrushes)

    Ahorn.redrawLayer!(toolsLayer)
end

function setMaterials!(layer::Ahorn.Layer)
    validTiles = Ahorn.validTiles(layer)
    tileNames = Ahorn.tileNames(layer)

    Ahorn.setMaterialList!([tileNames[mat] for mat in validTiles], row -> row[1] == get(tileNames, material, nothing))
end

function mouseMotion(x::Number, y::Number)
    if hoveringBrush === nothing || x != hoveringBrush[1] || y != hoveringBrush[2]
        global hoveringBrush = (x, y, deepcopy(selectedBrush))

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function middleClick(x::Number, y::Number)
    tiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    tileNames = Ahorn.tileNames(targetLayer)
    target = get(tiles.data, (y, x), '0')

    global material = target
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = material
    Ahorn.selectMaterialList!(tileNames[target])
end

function leftClick(x::Number, y::Number)
    layer = Ahorn.layerName(targetLayer)
    Ahorn.History.addSnapshot!(Ahorn.History.RoomSnapshot("Brush $(selectedBrush.name)", Ahorn.loadedState.room))

    Maple.updateTileSize!(Ahorn.loadedState.room, '0', '0')
    roomTiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    Ahorn.applyBrush!(selectedBrush, roomTiles, material, x, y)

    Ahorn.redrawLayer!(targetLayer)
end

function selectionMotion(x1::Number, y1::Number, x2::Number, y2::Number)
    box, boy = selectedBrush.offset

    startX = x1
    startY = y1

    pixels = rotr90(selectedBrush.pixels, selectedBrush.rotation)

    bh, bw = size(pixels)
    ox, oy = mod(startX, bw), mod(startY, bh)

    bx, by = div(x2, bw) * bw + ox - box + 1, div(y2, bh) * bh + oy - boy + 1

    if !haskey(phantomBrushes, (bx, by))
        phantomBrushes[(bx, by)] = deepcopy(selectedBrush)

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function selectionFinish(rect::Ahorn.Rectangle)
    applyPhantomBrushes()
end

function toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer)
    layerName = Ahorn.layerName(targetLayer)
    tileNames = Ahorn.tileNames(targetLayer)
    global material = get(Ahorn.persistence, "brushes_material_$(layerName)", tileNames["Air"])[1]

    wantedBrush = get(Ahorn.persistence, "brushes_brushes_brush", brushes[1].name)
    Ahorn.updateTreeView!(subTools, [brush.name for brush in brushes], row -> row[1] == wantedBrush)

    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")
    Ahorn.updateLayerList!(["fgTiles", "bgTiles"], row -> row[1] == wantedLayer)

    Ahorn.redrawingFuncs["tools"] = drawBrushes
    Ahorn.redrawLayer!(toolsLayer)
end

function layerSelected(list::Ahorn.ListContainer, materials::Ahorn.ListContainer, selected::String)
    global targetLayer = Ahorn.getLayerByName(drawingLayers, selected)
    Ahorn.persistence["brushes_layer"] = selected

    tileNames = Ahorn.tileNames(targetLayer)
    layerName = Ahorn.layerName(targetLayer)
    global material = get(Ahorn.persistence, "brushes_material_$(layerName)", tileNames["Air"])[1]
    setMaterials!(targetLayer)
end

function materialSelected(list::Ahorn.ListContainer, selected::String)
    tileNames = Ahorn.tileNames(targetLayer)
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = tileNames[selected]
    global material = tileNames[selected]
end

function materialFiltered(list::Ahorn.ListContainer)
    tileNames = Ahorn.tileNames(targetLayer)
    Ahorn.selectRow!(list, row -> row[1] == Ahorn.tileNames[material])
end

function layersChanged(layers::Array{Ahorn.Layer, 1})
    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")

    global drawingLayers = layers
    global toolsLayer = Ahorn.getLayerByName(layers, "tools")
    global targetLayer = Ahorn.selectLayer!(layers, wantedLayer, "fgTiles")
end

function subToolSelected(list::Ahorn.ListContainer, selected::String)
    for brush in brushes
        if brush.name == selected
            Ahorn.persistence["brushes_brushes_brush"] = selected
            global selectedBrush = brush
        end
    end
end

function keyboard(event::Ahorn.eventKey)
    shouldRedraw = false

    if event.keyval == Ahorn.keyval("l")
        selectedBrush.rotation = mod(selectedBrush.rotation + 1, 4)

        shouldRedraw |= true

    elseif event.keyval == Ahorn.keyval("r")
        selectedBrush.rotation = mod(selectedBrush.rotation - 1, 4)

        shouldRedraw |= true
    end

    if shouldRedraw && hoveringBrush !== nothing
        global hoveringBrush = (hoveringBrush[1], hoveringBrush[2], deepcopy(selectedBrush))

        Ahorn.redrawLayer!(toolsLayer)
    end
end

end

ext: jl

max_stars_repo_name: cdparks/Ahorn

max_issues_repo_name: cdparks/Ahorn

max_forks_repo_name: cdparks/Ahorn


text:
using Test
using LinearAlgebra
using LFPSQP
import ForwardDiff.Dual

# projected CG tests
include("test_cg.jl")

# gradient and Hessian generator function tests
include("test_autodiff.jl")

# test retraction validity (for equality constraints)
include("test_retractions.jl")

# test linesearches
include("test_linesearch.jl")

# test structs and associated functions to handle inequality constraints
include("test_inequalities.jl")


text:
@testset "Spaces" begin
  k = FlintQQ

  Qx, x = PolynomialRing(FlintQQ, "x")
  K1, a1 = NumberField(x^2 - 2, "a1")
  K2, a2 = NumberField(x^3 - 2, "a2")

  K1t, t = PolynomialRing(K1, "t")

  L, b = NumberField(t^2 + a1)

  for K in [k, K1, K2, L]
    V = @inferred quadratic_space(K, 2)
    @test_throws ArgumentError quadratic_space(K, -1)

    V = @inferred quadratic_space(K, identity_matrix(FlintZZ, 2))
    @test (@inferred gram_matrix(V)) == identity_matrix(K, 2)
    @test (@inferred dim(V)) == 2
    @test (@inferred rank(V)) == 2
    @test @inferred isregular(V)
    @test (@inferred involution(V)) === identity

    V = @inferred quadratic_space(K, zero_matrix(K, 2, 2))
    @test (@inferred gram_matrix(V)) == zero_matrix(K, 2, 2)
    @test (@inferred gram_matrix(V)) isa dense_matrix_type(K)
    @test (@inferred dim(V)) == 2
    @test (@inferred rank(V)) == 0
    @test @inferred isquadratic(V)
    @test @inferred !ishermitian(V)
    @test (@inferred fixed_field(V)) === K

    @test_throws ArgumentError quadratic_space(K, zero_matrix(K, 1, 2))

    M = identity_matrix(K, 2)
    M[1, 2] = 2
    M[2, 1] = 1
    @test_throws ArgumentError quadratic_space(K, M)
  end

end


text:
# write_data and write_info
"""
    write_data(data::Data, name_file::AbstractString)

Do this will write all information in `data` to file `name_file` in `.data` formation

# Example
```julia-repl
data_cell = genr_cell([10 10 3])
str = Si3N4()
data_atom = genr_atom(data_cell, str)
write_data(data_atom, "test.data")
```
"""
function write_data(data::Data, name_file::AbstractString)
    for info in fieldnames(typeof(data))
        write_info(getfield(data, info), name_file)
    end
end

function write_info(info::Data_Basic, name_file::AbstractString)
    # List Setting
    list_01 = "atoms bonds angles dihedrals impropers"
    list_02 = ["atom types", "bond types", "angle types", "dihedral types", "improper types"]
    list = vcat(split(list_01), list_02)
    num_paras = length(list)

    # Writting Para info
    io = open(name_file, "w")
    write(io, join(["Lammps .data file creat at ", Dates.now(), """ by Julia Package "lmp_str"\n\n"""]))

    fields = fieldnames(typeof(info))

    for para in 1 : num_paras
        write(io, join([string(getfield(info, fields[para])), "\t\t", list[para], "\n"]))
    end

    # Writing Box info

    para = num_paras + 1
    box_range = getfield(info, fields[para])
    write(io, join([string(box_range[1,1]), "\t\t", string(box_range[1,2]), "\t\t xlo xhi\n"]))
    write(io, join([string(box_range[2,1]), "\t\t", string(box_range[2,2]), "\t\t ylo yhi\n"]))
    write(io, join([string(box_range[3,1]), "\t\t", string(box_range[3,2]), "\t\t zlo zhi\n"]))

    para += 1
    box_tilt = getfield(info, fields[para])
    if !(length(findall(x->x==0, box_tilt)) == length(box_tilt))
        for tilt in box_tilt
            write(io, join([string(tilt), "\t"]))
        end
        write(io, "xy xz yz\n")
    end

    close(io)
end

function write_info(info::Str, name_file::AbstractString)
    # Reading Input
    fields = fieldnames(typeof(info))
    io = open(name_file, "a")

    # Judgement and Output
    if in(:para_mass, fields)
        write(io, "\n\nMasses\n\n")
        for atom = 1 : info.num_atom_types
            write(io, join([string(atom), "\t", string(info.para_mass[atom]), "\t\t# ", info.atom_name[atom], "\n"]))
        end
    end

    if in(:para_pair, fields)
        write(io, "\n\nPair Coeffs\n\n")
        for atom = 1 : info.num_atom_types
            para_now = info.para_pair[atom, :]
            write(io, join([string(atom)]))
            for para in para_now
                write(io, join(["\t\t", string(para)]))
            end
            write(io, "\n")
        end
    end

    if in(:para_bond, fields)
        if info.num_bond_types != 0
            write(io, "\n\nBond Coeffs\n\n")
            for bond = 1 : info.num_bond_types
                para_now = info.para_bond[bond, :]
                write(io, join([string(bond)]))
                for para in para_now
                    write(io, join(["\t\t", string(para)]))
                end
                write(io, "\n")
            end
        end
    end

    if in(:para_angle, fields)
        write(io, "\n\nAngle Coeffs\n\n")
        if info.num_angle_types != 0
            for angle = 1 : info.num_angle_types
                para_now = info.para_angle[angle, :]
                write(io, join([string(angle)]))
                for para in para_now
                    write(io, join(["\t\t", string(para)]))
                end
                write(io, "\n")
            end
        end
    end

    if in(:para_dihedral, fields)
        if info.num_dihedral_types != 0
            write(io, "\n\nDihedral Coeffs\n\n")
            for dihedral = 1 : info.num_dihedral_types
                para_now = info.para_dihedral[dihedral, :]
                write(io, join([string(dihedral)]))
                write(io, join(["\t\t", string(para_now[1])]))
                write(io, join(["\t\t", string(Int(para_now[2]))]))
                write(io, join(["\t\t", string(Int(para_now[3]))]))
                write(io, join(["\t\t", string(para_now[4])]))
                write(io, "\n")
            end
        end
    end

    if in(:para_improper, fields)
        if info.num_improper_types != 0
            write(io, "\n\nImproper Coeffs\n\n")
            for improper = 1 : info.num_improper_types
                para_now = info.para_improper[improper, :]
                write(io, join([string(improper)]))
                for para in para_now
                    write(io, join(["\t\t", string(para)]))
                end
                write(io, "\n")
            end
        end
    end

    write(io, "\n\n")

    close(io)
end

function write_info(info::Vector{Str}, name_file::AbstractString)
    # Reading Input
    fields = fieldnames(typeof(info[1]))
    io = open(name_file, "a")

    # Judgement and Output
    if in(:para_mass, fields)
        write(io, "\n\nMasses\n\n")
        id = 1
        for str = 1 : length(info)
            for atom = 1 : info[str].num_atom_types
                write(io, join([string(id), "\t", string(info[str].para_mass[atom]), "\t\t# ", info[str].atom_name[atom], "\n"]))
                id += 1
            end
        end
    end
    write(io, "\n\n")

    close(io)
end

function write_info(info::Vector{Atom}, name_file::AbstractString)
    # Reading Input
    io = open(name_file, "a")
    fields = fieldnames(typeof(info[1]))
    num_fields = length(fields)

    # Writing Output
    write(io, "Atoms # full\n\n")

    for atom in info
        for para in 1 : num_fields-1
            write(io, join([string(getfield(atom, fields[para])), " "]))
        end
        para = num_fields
        coord = getfield(atom,fields[para])
        for dim = 1 : 3
            write(io, join([string(coord[dim])," "]))
        end
        write(io, "\n")
    end
    write(io, "\n\n")

    close(io)
end

function write_info(info::Vector{T}, name_file::AbstractString) where T <: Union{Bond, Angle, Dihedral, Improper}
    # Reading Input
    io = open(name_file, "a")
    fields = fieldnames(typeof(info[1]))
    num_fields = length(fields)

    # Writing Output
    para = string(typeof(info[1]))
    write(io, join([para, "s\n\n"]))
    #para[findall(x->in('.', x), para)+1 : end]
    for item in info
        for para in 1 : num_fields - 1
            write(io, join([string(getfield(item, fields[para])), " "]))
        end
        atom = getfield(item, fields[end])
        for para in atom
            write(io, join([string(para), " "]))
        end
        write(io, "\n")
    end
    write(io, "\n\n")

    close(io)
end

function write_info(info::Int64, name_file::AbstractString)
    # do nothing for vec_bond or vec_angle = 0
    # It is very important although no code inside
end


text:
### Create data file
function write_data(output_path, dataset, Cnst, w_upper)

    open(output_path, "w") do io
        write(io, "param n := $(dataset.n);\n")
        write(io, "param m := $(dataset.m);\n")
        write(io, "param K := $(dataset.n_communities);\n")
        write(io, "\n")
        write(io, "param\tA:")

        for i=1:dataset.n
            write(io, "\t$i")
        end
        write(io, ":=")
        for i=1:dataset.n
            write(io, "\n\t\t$i")
            for j=1:dataset.n
                write(io, "\t$(dataset.A[i,j])")
            end
        end
        write(io, ";\n")

        write(io, "param\tk :=")
        for i=1:dataset.n
            write(io, "\n\t\t$i\t$(dataset.k[i])")
        end
        write(io, ";\n")

        write(io, "param\tCnst := $Cnst;\n")
        write(io, "param\tw_upper := $w_upper;\n")
    end;
end

### Create run file
function write_run_file(run_file, mod_file, data_file, log_file)
    open(run_file, "w+") do io
        write(io, "reset;\n")
        write(io, "model $mod_file;\n")
        write(io, "data $data_file;\n")
        write(io, "option log_file '$log_file';\n")
        write(io, "expand L, Constraint;\n")
        write(io, "option solver couenne;\n")
        write(io, "solve;\n")
        write(io, "display n;\n")
        write(io, "display m;\n")
        write(io, "display K;\n")
        write(io, "display A;\n")
        write(io, "display k;\n")
        write(io, "display Cnst;\n")
        write(io, "display w;\n")
        write(io, "display z;\n")
        write(io, "display L;\n")
    end
end

### Read log file and extract information
function extract_solution(dataset_name, log_file, sol_file)
    io = open(log_file)

    n = m = K = -1
    lincuts_root = lincuts_total = separation_time = total_solve_time = bb_time =\
        obj_lb = obj_ub = gap = bb_nodes = cnst = -1
    z = w = nothing
    z_mat = nothing  # Variable to track if z is written in a matrix or vector format
    status = "?"
    processing_x = false
    processing_w = false

    while !eof(io)
        line = readline(io)

        if startswith(line, "n = ")
            n = parse(Int, strip(split(line, "n = ")[2]))
        end
        if startswith(line, "m = ")
            m = parse(Int, strip(split(line, "m = ")[2]))
        end
        if startswith(line, "K = ")
            K = parse(Int, strip(split(line, "K = ")[2]))
        end
        if startswith(line, "couenne: Optimal")
            status = "OPTIMAL"
        end
        if startswith(line, "couenne: Optimization interrupted on limit.")
            status = "TIME_LIMIT"
        end
        if startswith(line, "Linearization cuts added at root node:")
            lincuts_root = parse(Int, strip(split(line, ":")[2]))
        end
        if startswith(line, "Linearization cuts added in total:")
            lincuts_total = parse(Int, strip(split(split(line, ":")[2], "(")[1]))
        end
        if startswith(line, "Linearization cuts added in total:")
            separation_time = parse(Float64, strip(split(split(line, "separation time:")[2], "s)")[1]))
        end
        if startswith(line, "Total solve time:")
            total_solve_time = parse(Float64, strip(split(split(line, ":")[2], "(")[1])[1:end-1])
        end
        if startswith(line, "Total solve time:")
            bb_time = parse(Float64, strip(split(split(line, "(")[2], "s in branch-and-bound)")[1]))
        end
        if startswith(line, "Lower bound:")
            obj_lb = parse(Float64, strip(split(line, ":")[2]))
        end
        if startswith(line, "Upper bound:")
            obj_ub = parse(Float64, strip(split(split(line, ":")[2], "(")[1]))
        end
        if startswith(line, "Upper bound:")
            gap = strip(split(split(line, "(gap: ")[2], ")")[1])
            if gap == "--"
                gap = NaN
            else
                gap = parse(Float64, gap[1:end-1])
            end
        end
        if startswith(line, "Branch-and-bound nodes:")
            bb_nodes = parse(Int, strip(split(line, ":")[2]))
        end
        if startswith(line, ";") && processing_x
            processing_x = false
        end
        if processing_x
            if z_mat
                if startswith(line, ":")
                else
                    elements = parse.(Float64, split(strip(line)))
                    println(elements)
                    i = Int(elements[1])
                    z[i,:] = elements[2:end]
                end
            else
                i,j,val = parse.(Float64, split(strip(line)))
                i = Int(i)
                j = Int(j)
                z[i,j] = val
            end
        end

        if startswith(line, "z :=")
            z = zeros(Float64, n, K)
            processing_x = true
            z_mat = false
        end
        if startswith(line, "z [*,*]")
            z = zeros(Float64, n, K)
            processing_x = true
            z_mat = true
        end
        if startswith(line, ";") && processing_w
            processing_w = false
        end
        if processing_w
            i,g,val = parse.(Float64, split(strip(line)))
            i = Int(i)
            g = Int(g)
            w[i,g] = val
        end
        if startswith(line, "w :=")
            w = zeros(Float64, K, K)
            processing_w = true
        end
        if startswith(line, "Cnst =")
            cnst = parse(Float64, strip(split(line, "=")[2]))
        end
    end
    close(io)

    println("Status: $status")
    println("Linearization cuts at root: $lincuts_root")
    println("Linearization cuts total: $lincuts_total")
    println("Separation time: $separation_time")
    println("Total solve time: $total_solve_time")
    println("B&B time: $bb_time")
    println("LB: $obj_lb")
    println("UB: $obj_ub")
    println("gap: $gap")
    println("B&B nodes: $bb_nodes")
    println("z = ")
    println(z)
    println("w = ")
    println(w)
    println("Cnst = $cnst")

    # Convert array type to Int
    z_int = convert.(Int, z)

    # Save solution to file
    open(sol_file, "w") do io
        write(io, "\tSolution\n")
        write(io, "dataset_name\t$dataset_name\n")
        write(io, "n\t$n\n")
        write(io, "m\t$m\n")
        write(io, "K\t$K\n")
        write(io, "Cnst\t$cnst\n")
        write(io, "method\tminlp\n")
        write(io, "LB\t$obj_lb\n")
        write(io, "UB\t$obj_ub\n")
        write(io, "status\t$status\n")
        write(io, "solvetime\t$total_solve_time\n")
        write(io, "bb_solvetime\t$bb_time\n")
        write(io, "nodecount\t$bb_nodes\n")
        write(io, "lincuts_root\t$lincuts_root\n")
        write(io, "lincuts_total\t$lincuts_total\n")
        write(io, "separation_time\t$separation_time\n")
        write(io, "gap\t$gap\n")
        write(io, "w\t$w\n")
        write(io, "z\t$z_int\n")
        write(io, "z_float\t$z\n") # TODO: convert to int/float?
    end
end

function MINLP_AMPL(dataset, dataset_name, mod_file, out_folder, ampl_path)

    println("Processing dataset:", dataset_name)

    data_file = string(out_folder, dataset_name, ".dat")

    Cnst = getObjectiveConstant(dataset)
    w_upper = getWUpperBound(dataset)
    write_data(data_file, dataset, Cnst, w_upper)

    method = split(basename(mod_file), ".")[1]

    ### Create run file
    run_file = string(out_folder, method, "-", dataset_name, ".run")
    log_file = string(out_folder, method, "-", dataset_name, ".log")

    write_run_file(run_file, mod_file, data_file, log_file)

    ### Run ampl couenne
    println("Running couenne")
    run(`$ampl_path $run_file`)

    log_file = string(out_folder, method, "-", dataset_name, ".log")

    # Extract solution details from log
    sol_file = string(out_folder, method, "-", dataset_name, ".out")
    extract_solution(dataset_name, log_file, sol_file)
    println()
end

max_stars_repo_name: vidalt/Exact-SBM

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM


text:
module MCTSParser

using Merlin
using ProgressMeter

include("game.jl")
include("players.jl")
include("mcts.jl")
include("game/tictactoe.jl")

end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	storage_asymmetric(EP::Model, inputs::Dict, Reserves::Int)

Sets up variables and constraints specific to storage resources with asymmetric charge and discharge capacities. See ```storage()``` in ```storage.jl``` for description of constraints.
"""
function storage_asymmetric(EP::Model, inputs::Dict, Reserves::Int)
	# Set up additional variables, constraints, and expressions associated with storage resources with asymmetric charge & discharge capacity
	# (e.g. most chemical, thermal, and mechanical storage options with distinct charge & discharge components/processes)
	# STOR = 2 corresponds to storage with distinct power and energy capacity decisions and distinct charge and discharge power capacity decisions/ratings

	println("Storage Resources with Asmymetric Charge/Discharge Capacity Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)

	STOR_ASYMMETRIC = inputs["STOR_ASYMMETRIC"]

	### Constraints ###

	# Storage discharge and charge power (and reserve contribution) related constraints for symmetric storage resources:
	if Reserves == 1
		EP = storage_asymmetric_reserves(EP, inputs)
	else
		# Maximum charging rate must be less than charge power rating
		@constraint(EP, [y in STOR_ASYMMETRIC, t in 1:T], EP[:vCHARGE][y,t] <= EP[:eTotalCapCharge][y])
	end

	return EP
end

@doc raw"""
	storage_asymmetric_reserves(EP::Model, inputs::Dict)

Sets up variables and constraints specific to storage resources with asymmetric charge and discharge capacities when reserves are modeled. See ```storage()``` in ```storage.jl``` for description of constraints.
"""
function storage_asymmetric_reserves(EP::Model, inputs::Dict)

	dfGen = inputs["dfGen"]
	T = inputs["T"]

	STOR_ASYMMETRIC = inputs["STOR_ASYMMETRIC"]

	STOR_ASYM_REG = intersect(STOR_ASYMMETRIC, inputs["REG"]) # Set of asymmetric storage resources with REG reserves
	STOR_ASYM_NO_REG = setdiff(STOR_ASYMMETRIC, STOR_ASYM_REG) # Set of asymmetric storage resources without REG reserves

	if !isempty(STOR_ASYM_REG)
		# Storage units charging can charge faster to provide reserves down and charge slower to provide reserves up
		# Maximum charging rate plus contribution to regulation down must be less than charge power rating
		@constraint(EP, [y in STOR_ASYM_REG, t in 1:T], EP[:vCHARGE][y,t]+EP[:vREG_charge][y,t] <= EP[:eTotalCapCharge][y])
	else
		# Storage units charging can charge faster to provide reserves down and charge slower to provide reserves up
		# Maximum charging rate plus contribution to regulation down must be less than charge power rating
		@constraint(EP, [y in STOR_ASYM_NO_REG, t in 1:T], EP[:vCHARGE][y,t] <= EP[:eTotalCapCharge][y])
	end

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
module StableRNGs

export StableRNG

using Random: Random, AbstractRNG, Sampler, SamplerType

import Random: rand, seed!

# implementation of LehmerRNG based on the constants found at the
# MIT licensed code by Melissa E. O'Neill at
# https://gist.github.com/imneme/aeae7628565f15fb3fef54be8533e39c

"""
    LehmerRNG
    StableRNG

Simple RNG with stable streams, usually suitable for testing.
Use only the alias `StableRNG`, as the name `LehmerRNG` is not
part of the API.

Construction: `StableRNG(seed::Integer)`.

Seeding: `Random.seed!(rng::StableRNG, seed::Integer)`.
"""
mutable struct LehmerRNG <: AbstractRNG
    state::UInt128

    LehmerRNG(seed::Integer) = seed!(new(), seed)

    function LehmerRNG(; state::UInt128)
        isodd(state) || throw(ArgumentError("state must be odd"))
        new(state)
    end
end

const StableRNG = LehmerRNG

function seed!(rng::LehmerRNG, seed::Integer)
    seed >= 0 || throw(ArgumentError("seed must be non-negative"))
    seed <= typemax(UInt64) ||
        # this constraint could be loosened a bit if requested
        throw(ArgumentError("seed must be <= $(typemax(UInt64))"))

    seed = ((seed % UInt128) << 1) | one(UInt128) # must be odd
    rng.state = seed
    rng
end

Base.show(io::IO, rng::LehmerRNG) =
    print(io, LehmerRNG, "(state=0x", string(rng.state, base=16, pad=32), ")")

function Base.copy!(dst::LehmerRNG, src::LehmerRNG)
    dst.state = src.state
    dst
end

Base.copy(src::LehmerRNG) = LehmerRNG(state=src.state)

Base.:(==)(x::LehmerRNG, y::LehmerRNG) = x.state == y.state

Base.hash(rng::LehmerRNG, h::UInt) = hash(rng.state, 0x93f376feff2bc48e % UInt ⊻ h)

## Sampling

function rand(rng::LehmerRNG, ::SamplerType{UInt64})
    rng.state *= 0x45a31efc5a35d971261fd0407a968add
    (rng.state >> 64) % UInt64
end

for T = [Bool, Base.BitInteger64_types...]
    T === UInt64 && continue
    @eval rand(rng::LehmerRNG, ::SamplerType{$T}) = rand(rng, UInt64) % $T
end

rand(rng::LehmerRNG, ::SamplerType{UInt128}) =
     rand(rng, UInt64) | ((rand(rng, UInt64) % UInt128) << 64)

rand(rng::LehmerRNG, ::SamplerType{Int128}) = rand(rng, UInt128) % Int128

Random.rng_native_52(::LehmerRNG) = UInt64

## within a range

# adapted verion of Random.SamplerRangeFast for native 64 bits generation
# we don't use "near division-less algorithm", as it doesn't work for Int128,
# and we want to maintain as little code as possible here

using Base: BitUnsigned, BitInteger
using Random: LessThan, Masked, uniform

struct SamplerRangeFast{U<:BitUnsigned,T<:BitInteger} <: Sampler{T}
    a::T      # first element of the range
    bw::UInt  # bit width
    m::U      # range length - 1
    mask::U   # mask generated values before threshold rejection
end

SamplerRangeFast(r::AbstractUnitRange{T}) where T<:BitInteger =
    SamplerRangeFast(r, T <: Base.BitInteger64 ? UInt64 : UInt128)

function SamplerRangeFast(r::AbstractUnitRange{T}, ::Type{U}) where {T,U}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    m = (last(r) - first(r)) % unsigned(T) % U
    #                        ^--- % unsigned(T) to not propagate sign bit
    bw = (sizeof(U) << 3 - leading_zeros(m)) % UInt # bit-width
    mask = ((1 % U) << bw) - (1 % U)
    SamplerRangeFast{U,T}(first(r), bw, m, mask)
end

function rand(rng::LehmerRNG, sp::SamplerRangeFast{UInt64,T}) where T
    a, bw, m, mask = sp.a, sp.bw, sp.m, sp.mask
    x = rand(rng, LessThan(m, Masked(mask, uniform(UInt64))))
    (x + a % UInt64) % T
end

function rand(rng::LehmerRNG, sp::SamplerRangeFast{UInt128,T}) where T
    a, bw, m, mask = sp.a, sp.bw, sp.m, sp.mask
    x = bw <= 64  ?
        rand(rng,
             LessThan(m % UInt64,
                      Masked(mask % UInt64, uniform(UInt64)))) % UInt128 :
        rand(rng, LessThan(m, Masked(mask, uniform(UInt128))))
    x % T + a
end

for T in Base.BitInteger_types
    # eval because of ambiguities with `where T <: BitInteger`
    @eval Sampler(::Type{LehmerRNG}, r::AbstractUnitRange{$T}, ::Random.Repetition) =
        SamplerRangeFast(r)
end

end # module

max_stars_repo_name: rfourquet/StableRNGs

max_issues_repo_name: rfourquet/StableRNGs

max_forks_repo_name: rfourquet/StableRNGs


text:
module DefaultConfigHexTest

using Hashids
using Test

@testset "Default Config Hex" begin

config = Hashids.configure()

@testset "encodehex and decodehex back" begin

@test decodehex(config, encodehex(config, "deadbeef")) == "deadbeef"
@test decodehex(config, encodehex(config, "abcdef123456")) == "abcdef123456"
@test decodehex(config, encodehex(config, "ABCDDD6666DDEEEEEEEEE")) == lowercase("ABCDDD6666DDEEEEEEEEE")
@test decodehex(config, encodehex(config, "507f1f77bcf86cd799439011")) == "507f1f77bcf86cd799439011"
@test decodehex(config, encodehex(config, "f00000fddddddeeeee4444444ababab")) == "f00000fddddddeeeee4444444ababab"
@test decodehex(config, encodehex(config, "abcdef123456abcdef123456abcdef123456")) == "abcdef123456abcdef123456abcdef123456"
@test decodehex(config, encodehex(config, "f000000000000000000000000000000000000000000000000000f")) == "f000000000000000000000000000000000000000000000000000f"
@test decodehex(config, encodehex(config, "fffffffffffffffffffffffffffffffffffffffffffffffffffff")) == "fffffffffffffffffffffffffffffffffffffffffffffffffffff"

end  # of "encodehex and decodehex back"

@testset "encodehex" begin

@test encodehex(config, "deadbeef") == "wpVL4j9g"
@test encodehex(config, "abcdef123456") == "kmP69lB3xv"
@test encodehex(config, "ABCDDD6666DDEEEEEEEEE") == "47JWg0kv4VU0G2KBO2"
@test encodehex(config, "507f1f77bcf86cd799439011") == "y42LW46J9luq3Xq9XMly"
@test encodehex(config, "f00000fddddddeeeee4444444ababab") == "m1rO8xBQNquXmLvmO65BUO9KQmj"
@test encodehex(config, "abcdef123456abcdef123456abcdef123456") == "wBlnMA23NLIQDgw7XxErc2mlNyAjpw"
@test encodehex(config, "f000000000000000000000000000000000000000000000000000f") == "VwLAoD9BqlT7xn4ZnBXJFmGZ51ZqrBhqrymEyvYLIP199"
@test encodehex(config, "fffffffffffffffffffffffffffffffffffffffffffffffffffff") == "nBrz1rYyV0C0XKNXxB54fWN0yNvVjlip7127Jo3ri0Pqw"

end  # of "encodehex"

end  # of "Default Config Hex"

end  # module


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

mutable struct UpdateFirewallRuleWithAccountParameters <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    properties::Any # spec type: Union{ Nothing, UpdateFirewallRuleProperties } # spec name: properties

    function UpdateFirewallRuleWithAccountParameters(;name=nothing, properties=nothing)
        o = new()
        validate_property(UpdateFirewallRuleWithAccountParameters, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(UpdateFirewallRuleWithAccountParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type UpdateFirewallRuleWithAccountParameters

const _property_map_UpdateFirewallRuleWithAccountParameters = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("properties")=>Symbol("properties"))
const _property_types_UpdateFirewallRuleWithAccountParameters = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("properties")=>"UpdateFirewallRuleProperties")
Base.propertynames(::Type{ UpdateFirewallRuleWithAccountParameters }) = collect(keys(_property_map_UpdateFirewallRuleWithAccountParameters))
Swagger.property_type(::Type{ UpdateFirewallRuleWithAccountParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpdateFirewallRuleWithAccountParameters[name]))}
Swagger.field_name(::Type{ UpdateFirewallRuleWithAccountParameters }, property_name::Symbol) =  _property_map_UpdateFirewallRuleWithAccountParameters[property_name]

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

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


text:
using Test

using CUDAdrv, CUDAnative
include(joinpath(@__DIR__, "..", "test", "array.jl"))   # real applications: use CuArrays.jl

function vadd(a, b, c)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    c[i] = a[i] + b[i]
    return
end

dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)
c = similar(a)

d_a = CuTestArray(a)
d_b = CuTestArray(b)
d_c = CuTestArray(c)

len = prod(dims)
@cuda threads=len vadd(d_a, d_b, d_c)
c = Array(d_c)
@test a+b ≈ c


text:
int main() {
    // Arithmetic operations with cross-types
    int x;
    double y;
    x%y;

    return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# generates precompile.jl files
using CSV, Parsers, SentinelArrays, SnoopCompile

const dir = joinpath(dirname(pathof(CSV)), "..", "test", "testfiles")

timings = @snoopi tmin=0.001 begin

f = CSV.File(joinpath(dir, "precompile.csv"))
f = CSV.File(joinpath(dir, "precompile_small.csv"))

end # @snoopi

pc = SnoopCompile.parcel(timings)
precom = joinpath(dirname(pathof(CSV)), "..", "precompile")
SnoopCompile.write(precom, pc)
mv(joinpath(precom, "precompile_CSV.jl"), joinpath(dirname(pathof(CSV)), "precompile.jl"); force=true)
mv(joinpath(precom, "precompile_Parsers.jl"), joinpath(dirname(pathof(Parsers)), "precompile.jl"); force=true)
mv(joinpath(precom, "precompile_SentinelArrays.jl"), joinpath(dirname(pathof(SentinelArrays)), "precompile.jl"); force=true)


text:
@testset "Flux Scattering Transform methods" begin
    ifGpu = identity
    i = 40
    s = 6 // 5
    nExtraDims = 2
    xExtraDims = 3
    k = 4
    N = 2
    i = 25
    s = 2 // 1
    nExtraDims = 2
    xExtraDims = 2
    k = 3
    N = 1
    import ScatteringTransform: stopAtExactly_WithRate_
    @testset "testing pooling" begin
        subsampRates = [3 // 2, 2, 5 // 2, 6 // 5]
        windowSize = [2, 3, 4]
        NdimsExtraDims = [(N, extra) for extra in 2:4, N in 1:2 if !(N == 2 && extra == 4)]
        @testset "bizzare pooling sizes 1D i=$i, s=$s, nExtraDims=$nExtraDims, xExtraDims=$xExtraDims, k=$k N=$N" for i = 25:3:40, s in subsampRates, xExtraDims in 2:nExtraDims, k in windowSize, (N, nExtraDims) in NdimsExtraDims
            selectXDims = ntuple(x -> 2, xExtraDims)
            x = ifGpu(randn(ntuple(x -> i, N)..., rand(2:10, xExtraDims)...))
            r = RationPool(ntuple(x -> s, N), k, nExtraDims=nExtraDims)
            @test length(r.m.k) == length(r.resSize) + nExtraDims - 2
            Nneed = ndims(r.m) + 2
            Sx = r(x)
            @test size(Sx) == (poolSize(r, size(x))..., size(x)[N+1:end]...)
            # does it take the mean we expect?
            chosenLocs = stopAtExactly_WithRate_(i - k + 1, s)
            loc = chosenLocs[5]
            neighborhood = Int.((1:k) .- ceil(k / 2))
            # the averaging neighborhood is a little bit hard to map, but it should
            # be centered correctly (possibly 1 larger, possibly 1 smaller)
            neighborhoodCent = ntuple(x -> loc .+ neighborhood, N)
            neighborhoodLowered = ntuple(x -> (loc - 1) .+ neighborhood, N)
            neighborhoodRaised = ntuple(x -> (loc + 1) .+ neighborhood, N)
            cent = Sx[ntuple(x -> 5, N)..., selectXDims...] ≈ mean(x[neighborhoodCent..., selectXDims...])
            lowered = Sx[ntuple(x -> 5, N)..., selectXDims...] ≈ mean(x[neighborhoodLowered..., selectXDims...])
            raised = Sx[ntuple(x -> 5, N)..., selectXDims...] ≈ mean(x[neighborhoodRaised..., selectXDims...])
            @test cent || lowered || raised
            # which entries matter for the 5th entry, according to the gradient?
            ∇ = gradient(x -> r(x)[ntuple(x -> 5, N)..., selectXDims...], x)
            δCent = zeros(size(x))
            δCent[neighborhoodCent..., selectXDims...] .= 1 / k^N
            δLowered = zeros(size(x))
            δLowered[neighborhoodLowered..., selectXDims...] .= 1 / k^N
            δRaised = zeros(size(x))
            δRaised[neighborhoodRaised..., selectXDims...] .= 1 / k^N
            ∇Locs = findall(∇[1] .!= 0)
            @test (∇Locs == findall(δCent .!= 0)) || (∇Locs == findall(δLowered .!= 0)) || (∇Locs == findall(δRaised .!= 0))
            @test ∇[1] ≈ δCent || ∇[1] ≈ δLowered || ∇[1] ≈ δRaised
        end
    end

    @test ScatteringTransform.nPathDims(1) == 1
    @test ScatteringTransform.nPathDims(2) == 1
    @test ScatteringTransform.nPathDims(3) == 2

    nFilters = [1, 9, 9, 9]
    @testset "1D basics" begin
        init = 10 .+ randn(64, 3, 2)
        sst = stFlux(size(init), 2, poolBy=3 // 2, outputPool=(2,))
        res = sst(init)
        @test length(res.output) == 2 + 1
        @test size(res.output[1]) == (32, 3, 2)
        @test minimum(abs.(res.output[1])) > 0
        @test size(res.output[2]) == (22, 3 * nFilters[2], 2)
        @test minimum(abs.(res.output[2])) > 0
        @test size(res.output[3]) == (14, nFilters[3], 3 * nFilters[2], 2)
        @test minimum(abs.(res.output[3])) > 0
        totalSize = 32 * 3 + 22 * 3 * nFilters[2] + 14 * nFilters[3] * 3 * nFilters[2]
        smooshed = ScatteringTransform.flatten(res)
        @test size(smooshed) == (totalSize, 2)
        sst1 = stFlux(size(init), 2, poolBy=3 // 2, outputPool=(2,), flatten=true)
        res1 = sst1(init)
        @test res1 isa Array{Float32,2}
        @test size(res1) == (32 * 3 + 22 * 3 * nFilters[2] + 14 * nFilters[3] * nFilters[2] * 3, 2)
        @test res1[1:32*3, 1] ≈ reshape(res[0][:, :, 1], (32 * 3,))

        sst = stFlux(size(init), 3, poolBy=3 // 2, outputPool=(2,))
        res = sst(init)
        @test length(res.output) == 3 + 1
        @test size(res.output[1]) == (32, 3, 2)
        @test minimum(abs.(res.output[1])) > 0
        @test size(res.output[2]) == (22, 3 * nFilters[2], 2)
        @test minimum(abs.(res.output[2])) > 0
        @test size(res.output[3]) == (14, nFilters[3], 3 * nFilters[2], 2)
        @test minimum(abs.(res.output[3])) > 0
        @test size(res.output[4]) == (10, nFilters[4], nFilters[3], 3 * nFilters[2], 2)
        @test minimum(abs.(res.output[4])) > 0
        totalSize = 32 * 3 + 22 * 3 * nFilters[2] + 14 * nFilters[3] * 3 * nFilters[2] + 10 * 3 * prod(nFilters)
        smooshed = ScatteringTransform.flatten(res)
        @test size(smooshed) == (totalSize, 2)
        sst1 = stFlux(size(init), 3, poolBy=3 // 2, outputPool=(2,), flatten=true)
        res1 = sst1(init)
        @test res1 isa Array{Float32,2}
        @test size(res1) == (32 * 3 + 22 * 3 * nFilters[2] + 14 * nFilters[3] * nFilters[2] * 3 + 10 * 3 * prod(nFilters), 2)
        @test res1[1:32*3, 1] ≈ reshape(res[0][:, :, 1], (32 * 3,))
    end

    nFilters = [1, 10, 9]
    @testset "1D integer pooling" begin
        stEx = stFlux((131, 1, 1), 2, poolBy=3)

        scat = stEx(ifGpu(randn(131, 1, 1)))
        @test size(stEx.mainChain[1].fftPlan[1]) == (2 * 131, 1, 1)
        @test size(stEx.mainChain[1].fftPlan[2]) == (2 * 131, 1, 1)
        @test size(stEx.mainChain[4].fftPlan[1]) == (2 * 44, nFilters[2], 1)
        @test size(stEx.mainChain[4].fftPlan[2]) == (2 * 44, nFilters[2], 1)
        @test size(stEx.mainChain[7].fftPlan[1]) == (2 * 15, nFilters[3], nFilters[2], 1)
        @test size(stEx.mainChain[7].fftPlan[2]) == (2 * 15, nFilters[3], nFilters[2], 1)
        @test stEx.mainChain[1].bc == FourierFilterFlux.Sym()
        @test stEx.mainChain[4].bc == FourierFilterFlux.Sym()
        @test stEx.mainChain[7].bc == FourierFilterFlux.Sym()
        @test stEx.outputPool == ((2,), (2,), (2,))
        @test ndims(stEx) == 1

        resultSize = ((66, 1, 1), (22, nFilters[2], 1), (8, nFilters[3], nFilters[2], 1))
        @test stEx.outputSizes == resultSize
        @test ([size(s) for s in scat.output]...,) == resultSize
        @test stEx.outputSizes == ([size(s) for s in scat.output]...,)
    end

    # rolling and flattening does nothing
    @testset "roll and flatten 1D" begin
        init = ifGpu(randn(64, 1, 2))
        sst = stFlux(size(init), 2, poolBy=3 // 2)
        res = sst(init)
        smooshed = ScatteringTransform.flatten(res)
        if ifGpu != identity
            @test typeof(smooshed) <: CuArray
        end
        reconst = roll(smooshed, sst)
        @test all(reconst .≈ res)
        @test typeof(reconst.output) <: Tuple
    end

    # normalization in 2D
    x = randn(10, 4, 3, 5, 7)
    Nd = 2
    xp = ScatteringTransform.normalize(x, 2)
    for w in eachslice(xp, dims=ndims(x))
        @test norm(w, 2) ≈ 3 * 5
    end

    # normalization in 1D
    x = randn(10, 3, 5, 7)
    xp = ScatteringTransform.normalize(x, 1)
    for w in eachslice(xp, dims=ndims(x))
        @test norm(w, 2) ≈ 3 * 5
    end
end
@testset "pathLocs" begin
    using ScatteringTransform: parseOne
    exs = Colon()
    # layer 0
    @test (34:53, 1, :) == parseOne((0, 34:53), 1, exs)
    @test (34:53, 34:53, 1, :) == parseOne((0, (34:53, 34:53)), 2, exs)
    @test (:, :, :) == parseOne((0, :), 1, exs)
    @test (:, :, :, :) == parseOne((0, :), 2, exs)
    @test (4:35, 1, 1:10) == parseOne((0, (4:35, 1, 1:10)), 1, exs)
    @test (4:35, 1, 1, 1:10) == parseOne((0, (4:35, 1, 1, 1:10)), 2, exs)
    # layer 1
    @test (:, 34:53, :) == parseOne((1, 34:53), 1, exs)
    @test (:, 5, :) == parseOne((1, 5), 1, exs)
    @test (:, :, 34:53, :) == parseOne((1, 34:53), 2, exs)
    @test (:, :, 5, :) == parseOne((1, 5), 2, exs)
    @test (:, :, :) == parseOne((1, :), 1, exs)
    @test (:, :, :, :) == parseOne((1, :), 2, exs)
    @test (4:35, 1, 1:10) == parseOne((1, (4:35, 1, 1:10)), 1, exs)
    @test (4:35, 1, 1, 1:10) == parseOne((1, (4:35, 1, 1, 1:10)), 2, exs)
    # layer 2
    @test (:, :, 34:53, :) == parseOne((2, 34:53), 1, exs)
    @test (:, :, :, 34:53, :) == parseOne((2, 34:53), 2, exs)
    @test (:, 34:53, 5, :) == parseOne((2, (34:53, 5)), 1, exs)
    @test (:, :, 34:53, 5, :) == parseOne((2, (34:53, 5)), 2, exs)
    @test (:, :, :, :) == parseOne((2, :), 1, exs)
    @test (:, :, :, :, :) == parseOne((2, :), 2, exs)
    @test (4:35, 1, 1:10, 1, :) == parseOne((2, (4:35, 1, 1:10)), 1, exs)
    @test (4:35, 1, 1, 1:10, 1, :) == parseOne((2, (4:35, 1, 1, 1:10)), 2, exs)

    # actual pathLocs construction
    @test ((:, :, :), (:, 3:4, :), (:, 5, :, :)) == pathLocs(1, 3:4, 0, :, 2, (5, :)).indices
    @test ((:, :, :), nothing, (:, 5, :, :)) == pathLocs(0, :, 2, (5, :)).indices

    # pathLocs usage
    ex = ScatteredOut((randn(50, 1, 1), randn(34, 13, 1), randn(23, 11, 13, 1)))
    p = pathLocs(0, :, 2, (5, :))
    @test minimum((ex.output[1], ex.output[3][:, 5, :, :]) .≈ ex[p])
    p = pathLocs(2, (4:9, 5))
    @test ex.output[3][:, 4:9, 5, :] ≈ ex[p]
    p = pathLocs(1, :)
    @test ex.output[2] ≈ ex[p]
    p = pathLocs()
    @test minimum(ex[p] .≈ ex.output)

    # setindex (some more in the nonZero paths section)
    p = pathLocs(2, (4:9, 5))
    newVal = randn(23, 6, 1)
    ex[p] = newVal
    @test ex[p] ≈ newVal

    p = pathLocs(2, (4:9, 5), 1, (3:5,))
    newVal = (randn(34, 3, 1), randn(23, 6, 1))
    ex[p] = newVal
    @test all(ex[p] .≈ newVal)

    # this is still broken; some sort of broadcast shenanigans
    p = pathLocs(0, 40:46)
    ex[p] .= 1
    #@test minimum(ex[p] .≈ 1) #TODO broken

    # single entries behave a bit strangely
    p = pathLocs(0, 3)
    newVal = randn()
    ex[p] = newVal
    @test minimum(ex[p] .≈ newVal)

    # getindex using ints and arrays
    @test ex[0] == ex.output[1]
    @test ex[1] == ex.output[2]
    @test ex[0:1] == ex.output[1:2]

    # pathLoc getindex adjoint
    p = pathLocs(2, (4:9, 5))
    y, back = pullback(x -> x[p], ex)
    res = back(y)
    @test res[1][p] == ex[p]
    anti_p = pathLocs(2, ([1:3..., 10, 11], [1:4... 6:13...]))
    @test res[1][anti_p] == zeros(23, 5, 1, 12, 1)
    @test maximum(maximum.(res[1][0:1])) == 0

    y, back = pullback(getindex, ex, 1)
    res = back(y)
    @test res[1][1] == ex[1]
    @test all(res[1][0] .== 0)
    @test all(res[1][2][1] .== 0)
    @test res[2] == nothing

    y, back = pullback(getindex, ex, 1:2)
    res = back(y)
    @test res[1][1:2] == ex[1:2]
    @test all(res[1][0] .== 0)
    @test res[2] == nothing

    # catting paths (may result in extra indices grabbed)
    pathsByHand = (pathLocs(0, 40:46, exs=1),
        pathLocs(1, (24:26, 4), exs=2),
        pathLocs(1, (29:31, 6), exs=1:2),
        pathLocs(2, (11:14, 10, 1), exs=1),
        pathLocs(2, (11:15, 3, 5:7), exs=1))
    joint = cat(pathsByHand...)

    # nonZeroPaths
    ex = ScatteredOut((zeros(50, 1, 2), zeros(34, 13, 2), zeros(23, 11, 13, 2)))
    for (ii, p) in enumerate(pathsByHand)
        ex[p] = (-1)^ii * ii * ones(size(ex[p]))
    end
    paths = nonZeroPaths(ex, wholePath=false, allTogetherInOne=true)
    @test ex[paths][1] == ex[pathLocs(0, 40:46, exs=1)]
    @test ex[paths][2] == [fill(-3.0, 3)...; fill(2.0, 3)...; fill(-3.0, 3)...] #not
    #easy to describe the location
    @test findall(paths.indices[2]) == [CartesianIndex(29, 6, 1), CartesianIndex(30, 6, 1), CartesianIndex(31, 6, 1), CartesianIndex(24, 4, 2), CartesianIndex(25, 4, 2), CartesianIndex(26, 4, 2), CartesianIndex(29, 6, 2), CartesianIndex(30, 6, 2), CartesianIndex(31, 6, 2)]
    @test findall(paths.indices[3]) == [CartesianIndex(11, 10, 1, 1), CartesianIndex(12, 10, 1, 1), CartesianIndex(13, 10, 1, 1), CartesianIndex(14, 10, 1, 1), CartesianIndex(11, 3, 5, 1), CartesianIndex(12, 3, 5, 1), CartesianIndex(13, 3, 5, 1), CartesianIndex(14, 3, 5, 1), CartesianIndex(15, 3, 5, 1), CartesianIndex(11, 3, 6, 1), CartesianIndex(12, 3, 6, 1), CartesianIndex(13, 3, 6, 1), CartesianIndex(14, 3, 6, 1), CartesianIndex(15, 3, 6, 1), CartesianIndex(11, 3, 7, 1), CartesianIndex(12, 3, 7, 1), CartesianIndex(13, 3, 7, 1), CartesianIndex(14, 3, 7, 1), CartesianIndex(15, 3, 7, 1)]
    @test ex[paths][3] == [fill(4.0, 4)...; fill(-5.0, 15)...]

    # setindex using a boolean array instead of a normal one
    mostlyNull = ScatteredOut((zeros(50, 1, 2), zeros(34, 13, 2), zeros(23, 11, 13, 2)))
    newValues = randn(35)
    mostlyNull[paths] #todo this should be a tuple not an array
    mostlyNull[paths] = newValues
    @test cat(mostlyNull[paths]..., dims=1) ≈ newValues
    # testing adding one location at a time
    nullEx = ScatteredOut((zeros(50, 1, 1), zeros(34, 13, 1), zeros(23, 11, 13, 1)))
    nullEx[0][3, 1, 1] = 50
    nullEx[2][3:20, 1, 1, 1] .= 25
    nullEx[2][5, 2:5, 1, 1] .= 26
    nullEx[2][5, 9, 12:13, 1] .= 2
    addFrom = nonZeroPaths(nullEx, allTogetherInOne=true, wholePath=false)
    addTo = addNextPath(addFrom)
    @test findfirst(addTo.indices[1] != addFrom.indices[1]) == nothing
    @test addTo.indices[2] == nothing
    @test findfirst(addTo.indices[3]) == nothing

    addToNext = addNextPath(addTo, addFrom)
    @test addToNext.indices[2] == nothing
    @test findfirst(addToNext.indices[1] != addFrom.indices[1]) == nothing
    @test findfirst(addToNext.indices[3] .!= addFrom.indices[3]) == CartesianIndex(4, 1, 1, 1)
    @test findall(addToNext.indices[3]) == findall(addFrom.indices[3])[1:1]
    ii = 2
    while addTo != addToNext && ii <= 24
        addTo = addToNext
        addToNext = addNextPath(addTo, addFrom)
        @test findall(addToNext.indices[3]) == findall(addFrom.indices[3])[1:ii]
        ii += 1
    end
    @test all(addToNext.indices .== addFrom.indices)
end
Sx = ScatteredOut((randn(16, 1, 1), randn(11, 32, 1), randn(7, 27, 32, 1)))
@testset "Scattering gradients" begin
    # testing Zygote.@adjoint function getindex(F::T, i::Integer) where T <: Scattered
    @testset "Integer Access" begin
        # second layer
        ∇ = gradient(Sx -> Sx[2][2, 3, 5, 1], Sx)[1]
        finalLayer = zeros(7, 27, 32, 1)
        finalLayer[2, 3, 5, 1] = 1
        almostAllZeros = ScatteredOut((zeros(16, 1, 1), zeros(11, 32, 1), finalLayer))
        @test ∇ == almostAllZeros
        # first layer
        ∇ = gradient(Sx -> Sx[1][5, 8, 1], Sx)[1]
        firstLayer = zeros(11, 32, 1)
        firstLayer[5, 8, 1] = 1
        almostAllZeros = ScatteredOut((zeros(16, 1, 1), firstLayer, zeros(7, 27, 32, 1)))
        @test ∇ == almostAllZeros
        # zeroth layer
        ∇ = gradient(Sx -> Sx[0][3, 1, 1], Sx)[1]
        zeroLayer = zeros(16, 1, 1)
        zeroLayer[3, 1, 1] = 1
        almostAllZeros = ScatteredOut((zeroLayer, zeros(11, 32, 1), zeros(7, 27, 32, 1)))
        @test ∇ == almostAllZeros
    end
    @testset "PathLocs Access" begin
        p = pathLocs(0, (5))
        ∇ = gradient(Sx -> Sx[p][1], Sx)[1]
        @test ∇[0][5] ≈ 1
        p = pathLocs(0, (:))
        ∇ = gradient(Sx -> Sx[p][4], Sx)[1]
        @test ∇[0][4] ≈ 1
        y, ∂ = pullback(Sx -> Sx[p], Sx)
        @test ∂(y)[1][p] ≈ y
    end
    sampleData = (randn(10, 1, 3), randn(10, 5, 1, 3), randn(10, 5, 5, 1, 3))
    y, Δ = pullback(x -> ScatteredOut(x), sampleData)
    @test Δ(y)[1] == sampleData
end


text:
using DataFrames: DataFrame
# using Plots: scatter, plot!, default
# default(framestyle = :box, minorticks = 4)

# check whether the data are valid for a fit
function check_data(xdata, ydata, errdata)
    if any(ismissing, xdata) || any(ismissing, ydata) || any(ismissing, errdata)
        error("Error: Data contain `missing` values.")
    end
    if any(isinf, xdata) || any(isinf, ydata) || any(isnan, xdata) || any(isnan, ydata) || any(isnan, errdata)
        error("Error: Data contain `Inf` or `NaN` values.")
    end
    if any(iszero, errdata)
        error("Error: Data contain 0 in the errors.")
    end
end

"""
    Data(x::T, y::T, err::T) where {T<:Vector{Real}}
    Data(df::DataFrame)

Fields: `x, y, err, ndata`

This defines a type for data with three columns:` x, y, err`; `ndata` is the number of data rows.
Different `Data` sets can be concatenated as `vcat(dat1, dat2, dat3)`.

Only symmetric errors (of `y`) are supported.
"""
struct Data
    x::Vector{Float64}
    y::Vector{Float64}
    err::Vector{Float64}
    ndata::Int
    function Data(x, y, err) # Data(x::Vector{T}, y::Vector{T}, err::Vector{T}) where {T<:Real}
        check_data(x, y, err)
        ndata = length(x)
        new(x, y, err, ndata)
    end
end

Data(df::DataFrame) = Data(df[:,1], df[:,2], df[:,3])

import Base.vcat, Base.getindex
vcat(dat1::Data, dat2::Data) = Data(vcat(dat1.x, dat2.x), vcat(dat1.y, dat2.y), vcat(dat1.err, dat2.err))
vcat(dat1::Data, dat...) = reduce(vcat, [dat1, dat...])
getindex(data::Data, in) = Data(map(x-> getindex(x, in), (data.x, data.y, data.err))...)

"""
    @plt_data(data, kws...)
    @plt_data!(data, kws...)

Convenient mascros to make an errorbar plot of the `data`; all combinations of
keyword settings for `scatter` in `Plots` can be used for the optional arguments `kws...`
"""
macro plt_data(data, kws...)
    _plt = quote
        if isempty($kws)
            Plots.scatter($data.x, $data.y, yerror = $data.err,
            xlab = "x", ylab = "y", label = "Data")
        else
            Plots.scatter($data.x, $data.y, yerror = $data.err; $(kws...) )
        end
    end
    return esc(_plt)
end
macro plt_data!(data, kws...)
    _plt = quote
        if isempty($kws)
            Plots.scatter!($data.x, $data.y, yerror = $data.err,
            xlab = "x", ylab = "y", label = "Data")
        else
            Plots.scatter!($data.x, $data.y, yerror = $data.err; $(kws...) )
        end
    end
    return esc(_plt)
end

@doc raw"""
    chisq(dist::Function, data, par; fitrange = ())

defines the ``\chi^2`` function: `fun` the function to be fitted to the data given by `data`.
The parameters are collected into `par`, given as an array or a tuple.
* `data` can be either given as the `Data` type, or of the form `(xdata, ydata [, err])`.
If no `err` is given explicitly, the errors are assumed to be 1 for all data points.
* `fitrange`: default to the whole data set; may given as, e.g., `2:10`,
which means only fitting to the 2nd to the 10th data points.
"""
function chisq(dist::Function, data::Data, par; fitrange = ())
    fitrange = (isempty(fitrange) ? (1:data.ndata) : fitrange)
    # __dist__(x, par) = (length(func_argnames(dist)) > 2 ? dist(x, par...) : dist(x, par) )
    res = 0.0
    @simd for i = fitrange[1]:fitrange[end]
        @inbounds res += ( (data.y[i]- dist(data.x[i], par))/data.err[i] )^2
    end
    return res
end
# function chisq(dist::Function, data::Data, par::Tuple; fitrange = ())
#     fitrange = (isempty(fitrange) && 1:data.ndata)
#     res = 0.0
#     @simd for i = fitrange
#         @inbounds res += ( (data.y[i]- dist(data.x[i], par))/data.err[i] )^2
#     end
#     return res
# end
function chisq(dist::Function, data, par; fitrange = ())
    _x = data[1]; _y = data[2]
    _n = length(_x)
    _err = (length(data) == 2 ? ones(_n) : data[3])
    fitrange = (isempty(fitrange) ? (1:data.ndata) : fitrange)
    # __dist__(x, par) = (length(func_argnames(dist)) > 2 ? dist(x, par...) : dist(x, par) )
    res = 0.0
    @simd for i = fitrange[1]:fitrange[end]
        @inbounds res += ( (_y[i]- dist(_x[i], par))/_err[i] )^2
    end
    return res
end

"""
    @plt_best(dist, fit, data, kws...)
    @plt_best!(dist, fit, data, kws...)

A convenient macro for comparing the best-fit result with the data; all combinations of
keyword settings for `plot` in `Plots` can be used for the optional arguments `kws...`
The ordering of `dist`, `fit`, and `data` does not matter.
"""
macro plt_best(dist, fit, data, kws...)
    _expr = quote
        _arr = [$dist, $fit, $data]
        _fit = _arr[findfirst(x->typeof(x) <:AbstractFit, _arr)]
        _dist = _arr[findfirst(x->typeof(x) <:Function, _arr)]
        _data = _arr[findfirst(x->typeof(x) <:Data, _arr)]

        _npts = 100
        _paras = args(_fit)
        _dis(x) = (typeof(_fit) == ArrayFit ? _dist(x, _paras) : _dist(x, _paras...) )
        _xrange = _data.x
        _wv = LinRange(_xrange[1], _xrange[end], _npts)

        Plots.scatter($data.x, $data.y, yerror = $data.err, label = "Data")
        if isempty($kws)
            Plots.plot!(_wv, _dis.(_wv), xlab = "x", ylab = "y", label = "Best fit", lw = 1.5 )
        else
            Plots.plot!(_wv, _dis.(_wv); $(kws...)  )
        end
    end
    esc( _expr )
end
macro plt_best!(dist, fit, data, kws...)
    _expr = quote
        _arr = [$dist, $fit, $data]
        _fit = _arr[findfirst(x->typeof(x) <:AbstractFit, _arr)]
        _dist = _arr[findfirst(x->typeof(x) <:Function, _arr)]
        _data = _arr[findfirst(x->typeof(x) <:Data, _arr)]

        _npts = 100
        _paras = args(_fit)
        _dis(x) = (typeof(_fit) == ArrayFit ? _dist(x, _paras) : _dist(x, _paras...) )
        _xrange = _data.x
        _wv = LinRange(_xrange[1], _xrange[end], _npts)

        Plots.scatter!($data.x, $data.y, yerror = $data.err, label = "Data")
        if isempty($kws)
            Plots.plot!(_wv, _dis.(_wv), xlab = "x", ylab = "y", label = "Best fit", lw = 1.5 )
        else
            Plots.plot!(_wv, _dis.(_wv); $(kws...)  )
        end
    end
    esc( _expr )
end


text:
# check the compatibility to the original docopt

using DocOpt
using Compat

immutable Token
    kind::Symbol
    value
end

immutable UserError; end

const user_error = "user-error"

function readuntil_or_error(io, delim)
    v = readuntil(io, delim)
    if !endswith(v, delim)
        error("'$delim' is expected")
    end
    # strip delim
    v[1:end-length(delim)]
end

function expects(io, s)
    for i in 1:length(s)
        @assert read(io, Char) == s[i]
    end
end

function read_char(io)
    c = read(io, Char)
    while c == '#'
        readuntil(io, '\n')
        c = read(io, Char)
    end
    c
end

# tokenizer of testcases of docopt
function token(io)
    # read char
    eof(io) && return Token(:eof, ())
    char = read_char(io)
    while isspace(char)
        eof(io) && return Token(:eof, ())
        char = read_char(io)
    end

    value = ()
    if char == 'r'
        kind = :docstring
        expects(io, "\"\"\"")
        value = readuntil_or_error(io, "\"\"\"")
    elseif char == '$'
        kind = :argument
        expects(io, " ")
        value = readuntil_or_error(io, '\n')
    elseif char == '{'
        kind = :lbrace
    elseif char == '}'
        kind = :rbrace
    elseif char == '['
        kind = :lbracket
    elseif char == ']'
        kind = :rbracket
    elseif char == ':'
        kind = :colon
    elseif char == ','
        kind = :comma
    elseif char == '"'
        kind = :string
        value = readuntil_or_error(io, '"')
    elseif char == 't'
        kind = :bool
        expects(io, "rue")
        value = true
    elseif char == 'f'
        kind = :bool
        expects(io, "alse")
        value = false
    elseif char == 'n'
        kind = :null
        expects(io, "ull")
        value = nothing
    elseif isdigit(char)
        kind = :number
        value = parse(Int, string(char))
    else
        error("unknown token")
    end
    Token(kind, value)
end

# parsers take an IO stream and a current token, then return a parsed value and a next token

function parse_testcases(io, t)
    testcases = Any[]
    while t.kind !== :eof
        testcase, t = parse_testcase(io, t)
        push!(testcases, testcase)
    end
    testcases, t
end

function parse_testcase(io, t)
    @assert t.kind === :docstring
    docstring = t.value
    t = token(io)
    pairs, t = parse_pairs(io, t)
    (docstring, pairs), t
end

function parse_pairs(io, t)
    @assert t.kind === :argument
    pairs = Any[]
    while t.kind === :argument
        arg = t.value
        res, t = parse_result(io, token(io))
        push!(pairs, (arg, res))
    end
    pairs, t
end

function parse_result(io, t)
    if t.kind === :string
        @assert t.value == user_error
        return UserError, token(io)
    end
    @assert t.kind === :lbrace
    result = Dict{String,Any}()
    t = token(io)
    while t.kind !== :rbrace
        (key, value), t = parse_keyvalue(io, t)
        result[key] = value
        t.kind === :comma && (t = token(io))
    end
    @assert t.kind === :rbrace
    result, token(io)
end

function parse_keyvalue(io, t)
    @assert t.kind === :string
    key = t.value
    t = token(io)
    @assert t.kind === :colon
    t = token(io)
    value, t = parse_value(io, t)
    (key, value), t
end

function parse_value(io, t)
    if t.kind ∈ [:string, :bool, :null, :number]
        return t.value, token(io)
    elseif t.kind === :lbracket
        values = Any[]
        t = token(io)
        while t.kind !== :rbracket
            value, t = parse_value(io, t)
            push!(values, value)
            t.kind === :comma && (t = token(io))
        end
        return values, token(io)
    end
    error("invalid token where a value is expected: $t")
end

let
    # language-agnostic test cases from:
    #   https://github.com/docopt/docopt/blob/master/testcases.docopt
    open("test/testcases.docopt") do f
        t = token(f)
        testcases, _ = parse_testcases(f, t)
        for (docstring, pairs) in testcases
            println("-"^60)
            print(docstring)
            for (arg, expected) in pairs
                println(arg)
                args = split(replace(arg, r"prog\s*", ""), r"\s+")
                args = filter!(a -> !isempty(a), args)
                local result
                try
                    result = docopt(docstring, args, exit_on_error=false)
                catch ex
                    if isa(ex, DocOpt.DocOptExit)
                        result = UserError
                    else
                        rethrow()
                    end
                end
                if result == expected
                    println("=> OK")
                else
                    println("=> NOT OK")
                    @show result expected
                end
                println()
            end
        end
    end
end


text:
module Utils

using Distributions, Discretizers, Optim, CSV, Plots
import Printf
import KernelDensity; const KDE = KernelDensity

"""
Function to perform the MCMC metropolis algorithm.
"""
function mcmc_metropolis(x0::AbstractArray, logPfunc::Function, Lchain::Integer;
                         propStd::Union{AbstractFloat,AbstractArray}=0.1,
                         step::Integer=500, burn::Integer=500, printFreq::Integer=10000,
                         prior=:none, verbose::Bool=true, scaleProp::Bool=true)

	if length(size(x0)) > 1 # restart from old chain
        if verbose; println("Restarting from old chain"); end
		xOld = x0[end,:]
		n = length(xOld)
	else
	    xOld = x0
	    n = length(x0)
	end
    chain = zeros(Float64, Lchain,n)
    chain[1,:] = xOld
    acc = 0

    if prior == :none
    	logpOld = logPfunc(xOld)
    else
        logpOld = logPfunc(xOld)
        for  (ip,prr) in enumerate(prior)
        	logpOld += log(Distributions.pdf(prr,xOld[ip]))
        end
    end

    for ii=2:Lchain
        if scaleProp
            proposal = MvNormal(propStd.*xOld)
        else
            if isa(propStd, Array)
                proposal = MvNormal(propStd)
            else
                proposal = MvNormal(propStd.*ones(size(xOld)))
            end
        end
        xNew = xOld + Distributions.rand(proposal)

        # Reject if any x are negative
        if any(x->x<0, xNew); continue; end

        # Evaluate new log-likelihood
        if prior == :none
	    	logpNew = logPfunc(xNew)
	    else
	        logpNew = 0.0
	        for  (ip,prr) in enumerate(prior)
	        	logpNew += log(Distributions.pdf(prr,xNew[ip]))
	        end
	        if !isinf(logpNew); logpNew += logPfunc(xNew); end
	    end
        a = exp(logpNew - logpOld)

        if Base.rand() < a
            xOld = xNew
            logpOld = logpNew
            acc += 1
        end
        chain[ii,:] = xOld

        if ii % printFreq == 0
            if verbose
                Printf.@printf("Completed iteration %i out of %i. \n", ii,Lchain)
            end
        end
    end

    if verbose
        Printf.@printf("Acceptance ratio of %.2f (%i out of %i).\n", acc/Lchain,acc,Lchain)
    end

	if length(size(x0)) > 1
		chainRed = [x0; chain[step:step:end,:]] # Append new chain to old
	else
		chainRed = chain[burn:step:end,:]
	end
	return chainRed

end

"""
Function to find the MAP from the posterior distribution.
Finds independently for each parameter.
"""
function find_MAP(chain; idx=1)

    if length(size(chain))>1
        data = chain[:,idx]
    else
        data = copy(chain)
    end

    # Perform the KDE
    kdeObj = kde_wrpr(data)

    # Find MAP via crude optimisation
    func = x->KDE.pdf(kdeObj,x[1])
    x = collect(range(min(data...),stop=max(data...),length=10000))
    y = map(func,x)
    MAP = x[y.==maximum(y)][1]

    return MAP

end

"""
Function to find the MAP from the posterior distribution.
Finds independently for each parameter.
"""
function credibleintervals(chain; idx=1, spread=0.682)

    if length(size(chain))>1
        data = chain[:,idx]
    else
        data = copy(chain)
    end

    # Perform the KDE
    kdeObj = kde_wrpr(data)
    func = z->KDE.pdf(kdeObj,z)
    x = collect(range(min(data...),stop=max(data...),length=10000))
    dx = x[2]-x[1]
    y = map(func,x)

    # Move out from MAP in both directions
    # Stop once intervals contain "spread" of probability
    MAP = find_MAP(data)
    xMin = last(x[x.<MAP])
    xMax = first(x[x.>MAP])
    tot = 0.5*dx*(KDE.pdf(kdeObj,xMin)+KDE.pdf(kdeObj,xMax))
    flag = false
    while tot<spread
        xMinN = xMin-dx
        xMaxN = xMax+dx
        if xMinN<0.0
            flag = true
            break
        end
        if KDE.pdf(kdeObj,xMinN)>KDE.pdf(kdeObj,xMaxN)
            tot += 0.5*dx*(KDE.pdf(kdeObj,xMin)+KDE.pdf(kdeObj,xMinN))
            xMin = xMinN
        else
            tot += 0.5*dx*(KDE.pdf(kdeObj,xMax)+KDE.pdf(kdeObj,xMaxN))
            xMax = xMaxN
        end
    end

    if flag
        sort!(data)
        idx = Int(round(spread*length(data)))
        xMin = 0.0
        xMax = data[idx]
    end

    return xMin,xMax

end

"""
Wrapper around the KDE function to ensure that all calls to KDE use the same settings
"""
function kde_wrpr(data)
    KDE.kde(data, boundary=(0.0,2*max(data...)))
end

"""
Function to plot the full chains and posterior distributions.
"""
function plot_chain(chain, MAP=:none, Ints=:none)

	nVar = size(chain)[2]
	plots = Array{Any,1}(undef,nVar*2)

    # Set MAP to the mean if not provided
    if isequal(MAP,:none)
        MAP = mean(chain,dims=1)
        lbl = "Mean"
    else
        lbl = "MAP"
    end

	for ii=1:nVar
		plots[ii] = plot(chain[:,ii], legend=false, title=Printf.@sprintf("Parameter %i",ii));

		smpls = chain[:,ii]
		kdeObj = kde_wrpr(smpls)
		x = collect(range(Base.minimum(smpls),stop=Base.maximum(smpls),length=100))
		y = map(z->KDE.pdf(kdeObj,z),x)
		if ii ==1
			tmp = plot(x,y, label="Posterior", legend=:top);
			plot!(tmp, [MAP[ii],MAP[ii]],[0,Base.maximum(y)], label=lbl, line=(2,:black));
		else
			tmp = plot(x,y, legend=false);
			plot!(tmp, [MAP[ii],MAP[ii]],[0,Base.maximum(y)], line=(2,:black));
		end
        if !isequal(Ints,:none)
            plot!(tmp, [Ints[ii][1],Ints[ii][1]],[0,Base.maximum(y)], line=(2,:black,:dash), label="");
            plot!(tmp, [Ints[ii][2],Ints[ii][2]],[0,Base.maximum(y)], line=(2,:black,:dash), label="");
        end
		plots[ii+nVar] = tmp;

	end

	if nVar == 2
		plt = plot(plots[1],plots[2],plots[3],plots[4],
			       layout=(2,nVar), size=(800,800))
	elseif nVar == 3
		plt = plot(plots[1],plots[2],plots[3],plots[4],plots[5],plots[6],
			       layout=(2,nVar), size=(1200,800))
	elseif nVar == 4
		plt = plot(plots[1],plots[2],plots[3],plots[4],plots[5],plots[6],plots[7],plots[8],
			       layout=(2,nVar), size=(1600,800))
	elseif nVar == 5
		plt = plot(plots[1],plots[2],plots[3],plots[4],plots[5],plots[6],plots[7],plots[8],plots[9],plots[10],
			       layout=(2,nVar), size=(2000,800))
	elseif nVar == 6
		plt = plot(plots[1],plots[2],plots[3],plots[4],plots[5],plots[6],plots[7],plots[8],plots[9],plots[10],plots[11],plots[12],
			       layout=(2,nVar), size=(2400,800))
	else
		Printf.@printf("nVar of %i not currently catered for.", nVar)
		plt = plot()
	end

	return plt

end

"""
Function to load in the experimental data, filtering out missings/NaNs and
applying upper cut-off in terms of standard deviations from the mean.
"""
function load_data(File::String, Folder::String, cutOff::Number=Inf, format=:Int)

    if occursin(".csv",File)
        rnaData = CSV.read(Folder*File, datarow=1)[:,1]
    else
	    rnaData = CSV.read(Folder*File*".csv", datarow=1)[:,1]
    end

	filter!(x -> !isnan(x), rnaData)

    fltData = filter(x->x<cutOff, rnaData)

    if format==:Int
	    return Integer.(round.(fltData))
    else
        return fltData
    end

end

"""
1D integer version of genpdf, intended for use when the underlying data is discrete.
"""
function genpdf(data::Array{Int64,1}, nbin=:auto::Union{Int,Symbol})

    # Set bins to all the integers between lo and hi
    lo, hi = extrema(data)
    edges = collect(lo-1:hi) .+ 0.5

    disc = LinearDiscretizer(edges)
    counts = get_discretization_counts(disc, data)

    x = edges[1:end-1]+0.5*diff(edges)
    y = (counts./(diff(edges)*float(size(data)[1]))) # Normalise
    return (x,y)

end

"""
Function to remove the spurious zeros that can arise during MCMC.
"""
function rmv_zeros(chain)

    l = length(chain[:,1])
    cnt = 0
    for ii=1:l
        if iszero(chain[ii,3])
            cnt+=1
        end
    end

    chain_red = Array{Float64,2}(undef, l-cnt,3)
    idx = 1
    for ii=1:l
        if !iszero(chain[ii,3])
            chain_red[idx,:] = chain[ii,:]
            idx+=1
        end
    end
    return chain_red

end

"""
Return log-likelihood of getting data X from zero-inflated negative binomial with prms
"""
function log_likelihood_zi_negbinom(X, prms, verbose=false)

    if verbose
        println(prms)
    end

    m = MixtureModel([DiscreteUniform(0,0), NegativeBinomial(prms[1:2]...)], [prms[3],1-prms[3]])
    sum(log.([pdf(m,X[ii]) for ii=1:length(X)]))
end

"""
Obtain confidence intervals for the statistic specified in statFunc using bootstrapping
"""
function bootstrap_intervals(data, statFunc::Function, N=10000)

	statistic = statFunc(data)
	nSamp = length(data)

	statVals = Array{Float64}(undef, N)
	for ii=1:N
		# Resample from data with replacement
		idxs = Int.(round.(rand(nSamp).*nSamp .+0.5))
		samps = data[idxs]
		statVals[ii] = statFunc(samps)
	end

	ints = std(statVals)
	return statistic-ints,statistic,statistic+ints

end

end # module

max_stars_repo_name: rdbrackston/tx-analysis

max_issues_repo_name: rdbrackston/tx-analysis

max_forks_repo_name: rdbrackston/tx-analysis


text:
using Images: mapwindow, Fill

function part1(input)
    lut, image = parse_input(input)
    return sum(enhance(image, lut, 2))
end

function part2(input)
    lut, image = parse_input(input)
    return sum(enhance(image, lut, 50))
end

function parse_input(input)
    parts = split(readchomp(input), "\n\n")
    lut = collect(first(parts)) .== '#'
    image = reduce(hcat, (x -> Int.(collect(x) .== '#')).(split(last(parts), "\n")))
    return lut, image
end

function enhance(image, lut, n)
    for i = 1:n
        border = first(lut) & ((i + 1) % 2)
        padded_image = fill(border, size(image) .+ 2)
        padded_image[2:end-1, 2:end-1] .= image
        image = mapwindow(padded_image, (3, 3), border = Fill(border)) do x
            lut[1 + sum((1 .<< (8:-1:0)) .* vec(x))]
        end
    end
    return image
end


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

JLLWrappers.@generate_main_file_header("GnuTLS")
JLLWrappers.@generate_main_file("GnuTLS", UUID("0951126a-58fd-58f1-b5b3-b08c7c4a876d"))
end  # module GnuTLS_jll


text:
"""
SCOMPONENTS
-----------
    compute the strongly connected components of a graph

    ci=scomponents(A) returns an index for the component number of every
    vertex in the graph A.  The total number of components is maximum(ci).
    If the input is undirected, then this algorithm outputs just the
    connected components.  Otherwise, it output the strongly connected components.

    The implementation is from Tarjan's 1972 paper: Depth-first search and
    linear graph algorithms. In SIAM's Journal of Computing, 1972, 1,
    pp.146-160.

Functions
---------
- `cc = scomponents(A::MatrixNetwork)`
- `cc = scomponents{T}(A::SparseMatrixCSC{T,Int64})`
- `sci = strong_components_map(A::MatrixNetwork)`
- `sci = strong_components_map{T}(A::SparseMatrixCSC{T,Int64})`
- `sc_rich = enrich(cc::Strong_components_output) # check ?enrich for more`

Example
-------
~~~
A = load_matrix_network("cores_example")
cc = scomponents(A)
scomponents(A).number
scomponents(A).sizes
scomponents(A).map
strong_components_map(A)     # if you just want the map
enrich(scomponents(A)) # produce additional enriched output

# Can work on [ei,ej]
ei = [1;2;3]
ej = [2;4;1]
cc = scomponents(ei,ej)

# Can work on sparse matrix A
A = sprand(5,5,0.5)
cc = scomponents(A)
~~~
"""
function scomponents end

###########################
##    Type Definitions    #
###########################
mutable struct Strong_components_output
    map::Vector{Int64}    # The indicator map
    sizes::Vector{Int64}  # Array of sizes corresponding to map
    number::Int64         # Int64
    A::MatrixNetwork      # MatrixNetwork
end

mutable struct Strong_components_rich_output
    reduction_matrix::SparseMatrixCSC{Int64,Int64} # the reduction matrix (restriction matrix)
    transitive_order::SparseMatrixCSC{Int64,Int64}
    transitive_map::SparseMatrixCSC{Int64,Int64}
end

####################
##    Functions    #
####################

"""
Return information on the strongly connected components of a graph that
is the minimum required computation.
Example:
``strong_components_map(MatrixNetwork(sprand(5,4,0.1)))``
"""
function strong_components_map(A::MatrixNetwork)

    # TODO, remove the dt variable here
    n=A.n
    sci=zeros(Int64,n)
    cn=1
    root=zeros(Int64,n)
    dt=zeros(Int64,n)
    t=0
    cs=zeros(Int64,n)
    css=0 # component stack
    rs=zeros(Int64,2*n)
    rp = A.rp
    ci = A.ci
    rss=0 # recursion stack holds two nums (v,ri)

    # start dfs at 1
    for sv=1:n
        v=sv
        if root[v]>0
            continue
        end
        rss=rss+1
        rs[2*rss-1]=v
        rs[2*rss]=rp[v] #add v to the stack
        root[v]=v
        sci[v]=-1
        dt[v]=t
        t=t+1
        css=css+1
        cs[css]=v #add w to component stack
        while rss>0
            v=rs[2*rss-1] # pop v from the stack
            ri=rs[2*rss]
            rss=rss-1
            while ri<rp[v+1]
                w=ci[ri]
                ri=ri+1
                if root[w]==0
                    root[w]=w; sci[w]=-1;  dt[w]=t; t=t+1;
                    css=css+1 # add w to component stack
                    cs[css]=w
                    rss=rss+1
                    rs[2*rss-1]=v
                    rs[2*rss]=ri # add v to the stack
                    v=w
                    ri=rp[w]
                    continue
                end
            end
            for ri=rp[v]:rp[v+1]-1
                w=ci[ri]
                if sci[w]==-1
                    if dt[root[v]]>dt[root[w]]
                       root[v]=root[w]
                    end
                end
            end
            if root[v]==v
                while css>0
                    w=cs[css]
                    css=css-1
                    sci[w]=cn
                    if w==v
                        break
                    end
                end
                cn=cn+1
            end
        end
    end
    return sci
end

########################################################
##    Conversion Functions for strong_components_map   #
########################################################
# CSC:
strong_components_map(A::SparseMatrixCSC{T,Int64}) where {T} =
                                        strong_components_map(MatrixNetwork(A))
# Triplet:
strong_components_map(ei::Vector{Int64},ej::Vector{Int64}) =
                                        strong_components_map(MatrixNetwork(ei,ej))
# CSR
strong_components_map(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64) where {T} =
                                        strong_components_map(MatrixNetwork(n,rp,ci,vals))

######################
##    scomponents    #
######################

function scomponents(A::MatrixNetwork)
    mapping = strong_components_map(A)
    number = 0
    if length(mapping) > 0
    	number = maximum(mapping)
    end
	sizes = zeros(Int64,number)
    for i = 1:length(mapping)
        sizes[mapping[i]] += 1
    end

    return Strong_components_output(mapping, sizes, number, A)
end

###############################
##    Conversion Functions    #
###############################

# CSC:
scomponents(A::SparseMatrixCSC{T,Int64}) where {T} = scomponents(MatrixNetwork(A))
# Triplet:
scomponents(ei::Vector{Int64},ej::Vector{Int64}) = scomponents(MatrixNetwork(ei,ej))
# CSR
scomponents(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64) where {T} = scomponents(MatrixNetwork(n,rp,ci,vals))

"""
This function adds the following helpers variables
* reduction_matrix - a reduction matrix to project down to the component graph
* transitive_order - a transitive ordering of the components
* transitive_map - a map to components that respects the transitive ordering
* largest
"""
function enrich(rval::Strong_components_output)
    ci = rval.map
    sizes = rval.sizes
    ncomp = maximum(ci)
    R = sparse(collect(1:rval.A.n),ci,1,rval.A.n,ncomp)
    A = SparseMatrixCSC(rval.A.n,rval.A.n,rval.A.rp,rval.A.ci,rval.A.vals)
    CG = R'*A'*R
    return Strong_components_rich_output(R,A,CG)
end


text:
"""
$(SIGNATURES)

Takes an `MOI.TerminationStatusCode` and a `MOI.ResultStatusCode` and returns
the tuple `(valid_result::Bool, feasible::Bool)`. The value `valid_result` is
`true` if the pair of codes prove that either the subproblem solution was solved
to global optimality or the subproblem solution is infeasible. The value of
`feasible` is true if the problem is feasible and false if the problem is infeasible.
"""
function is_globally_optimal(t::MOI.TerminationStatusCode, r::MOI.ResultStatusCode)

    if t == MOI.INFEASIBLE && r == MOI.INFEASIBILITY_CERTIFICATE
        valid_result = true
        feasible = false

    elseif t == MOI.INFEASIBLE && r == MOI.NO_SOLUTION
        valid_result = true
        feasible = false

    elseif t == MOI.INFEASIBLE && r == MOI.UNKNOWN_RESULT_STATUS
        valid_result = true
        feasible = false

    elseif t == MOI.OPTIMAL && r == MOI.FEASIBLE_POINT
        valid_result = true
        feasible = true

    elseif t == MOI.INFEASIBLE_OR_UNBOUNDED && r == MOI.NO_SOLUTION
        valid_result = true
        feasible = false

    else
        valid_result = false
        feasible = false
    end

    return valid_result, feasible
end

"""
$(SIGNATURES)

Retrieves the lower and upper duals for variable bounds from the
`relaxed_optimizer` and sets the appropriate values in the
`_lower_lvd` and `_lower_uvd` storage fields.
"""
function _set_dual!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    opt = m.relaxed_optimizer
    foreach(t -> (m._lower_uvd[t[2]] = MOI.get(opt, MOI.ConstraintDual(), t[1]);), m._relaxed_variable_lt)
    foreach(t -> (m._lower_lvd[t[2]] = MOI.get(opt, MOI.ConstraintDual(), t[1]);), m._relaxed_variable_gt)
    return nothing
end

function interval_objective_bound(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    interval_objective_bound = bound_objective(m)
    if interval_objective_bound > m._lower_objective_value
        m._lower_objective_value = interval_objective_bound
        fill!(m._lower_lvd, zero(T))
        fill!(m._lower_uvd, zero(T))
        return true
    end
    return false
end

"""
$(SIGNATURES)

A fallback lower bounding problem that consists of an natural interval extension
calculation. This is called when the optimizer used to compute the lower bound
does not return a termination and primal status code indicating that it
successfully solved the relaxation to a globally optimal point.
"""
function fallback_interval_lower_bound!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    feas = true
    wp = _working_problem(m)
    if !cp_condition(m)
        feas = foreach_until(x -> _is_feas(Val{:leq}, m, x), wp._saf_leq, feas)
        feas = foreach_until(x -> _is_feas(Val{:eq}, m, x), wp._saf_eq, feas)
        feas = foreach_until(x -> _is_feas(m, x), wp._sqf_leq, feas)
        feas = foreach_until(x -> _is_feas(m, x), wp._sqf_eq,  feas)
        feas = foreach_until(x -> _is_feas(m, x), wp._nonlinear_constr, feas)
    end

    if feas
        interval_objective_used = interval_objective_bound(m)
        @__dot__ m._current_xref = _mid(BranchVar, m, 1:N)
        unsafe_check_fill!(isnan, m._current_xref, zero(T), N)
    else
        m._lower_objective_value = -Inf
    end
    m._lower_feasibility = feas

    return
end

function _update_box!(m::GlobalOptimizer{N,T,S}, v::Val{:cont}) where {N,T<:AbstractFloat,S}
    wp = _working_problem(m)
    lb = _lower_bound(FullVar, m)
    ub = _upper_bound(FullVar, m)
    opt = m.mip_optimizer
    foreach(x -> MOI.set(opt, MOI.ConstraintSet(), x[1], ET(lb[x[2]])),           m._relaxed_variable_eq)
    foreach(x -> MOI.set(opt, MOI.ConstraintSet(), x[1], LT(ub[x[2]])),           m._relaxed_variable_lt)
    foreach(x -> MOI.set(opt, MOI.ConstraintSet(), x[1], GT(ub[x[2]])),           m._relaxed_variable_gt)
    foreach(x -> MOI.set(opt, MOI.ConstraintSet(), x[1], IT(lb[x[2]], ub[x[2]])), m._relaxed_variable_it)
    foreach(x -> MOI.set(opt, MOI.ConstraintSet(), x[1], IT(zero(T), one(T))),    m._relaxed_variable_zo)
    foreach(x -> MOI.set(opt, MOI.ConstraintSet(), x[1], IT(lb[x[2]], ub[x[2]])), m._relaxed_variable_int)
    return nothing
end

function _update_box!(m::GlobalOptimizer{N,T,S}, v::Val{:mixed}) where {N,T<:AbstractFloat,S}
    wp = _working_problem(m)
    for x in m._relaxed_variable_zo
        l = _lower_bound(FullVar, m, x[2])
        u = _upper_bound(FullVar, m, x[2])
        if u < one(T)
            MOI.set(wp.mip_optimizer, MOI.ConstraintSet(), x[1], ET(zero(T)))
        elseif l > zero(T)
            MOI.set(wp.mip_optimizer, MOI.ConstraintSet(), x[1], ET(one(T)))
        else
            MOI.set(wp.mip_optimizer, MOI.ConstraintSet(), x[1], ZO())
        end
    end
    for x in m._relaxed_variable_int
        l = _lower_bound(FullVar, m, x[2])
        u = _upper_bound(FullVar, m, x[2])
        if iszero(l - floor(u))
            MOI.set(wp.mip_optimizer, MOI.ConstraintSet(), x[1], ET(l))
        elseif iszero(u - ceil(l))
            MOI.set(wp.mip_optimizer, MOI.ConstraintSet(), x[1], ET(u))
        else
            MOI.set(wp.mip_optimizer, MOI.ConstraintSet(), x[1], MOI.Semiinteger(l,u))
        end
    end
    return nothing
end

function _reset_objective!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    wp = _working_problem(m)
    MOI.set(wp.mip_optimizer, MOI.ObjectiveFunction{SAF}(), wp._objective)
    return nothing
end

"""
$(SIGNATURES)

Constructs and solves the relaxation using the default EAGO relaxation scheme
and optimizer on node `y`.
"""
function lower_problem!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    if !_relaxation_available(m)
        _update_box!(m, Val(:cont))
        _relax_problem!(m, Val(:cont))
    else
        _reset_objective!(m)
    end

    # solve LPs, adding linear cuts
    MOI.set(m.relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(m.relaxed_optimizer)
    valid_flag, feasible_flag = is_globally_optimal(m._lower_termination_status,
                                                    m._lower_result_status)

    # add continuous linear cuts
    linear_cut_condition = true
    while valid_flag && feasible_flag && linear_cut_condition
        _relax_problem!(m, Val(:cont))
        MOI.optimize!(m.relaxed_optimizer)
        valid_flag, feasible_flag = is_globally_optimal(m._lower_termination_status,
                                                        m._lower_result_status)
        linear_cut_condition = false # TODO: ADD CUT CHECK...
    end

    # activate integer variables and solve mixed integer formulation
    is_integer_prob = _integer_variable_num(m) > 0
    if valid_flag && feasible_flag && is_integer_prob
        _update_box!(m, Val(:mixed))      # need to update qp and conic relaxations as well
        _relax_problem!(m, Val(:mixed))
        MOI.optimize!(m.relaxed_optimizer)
    end

    if valid_flag && feasible_flag
        _set_dual!(m)
        m._lower_feasibility = true
        m._lower_objective_value = MOI.get(d, MOI.ObjectiveValue())
        result_count = MOI.get(d, MOI.ResultCount())
        for i = m._lower_result_count_max:(result_count-1)
            push!(m._lower_solution, zeros(m._working_problem._variable_num))
        end
        for j = 1:result_count
            for i = 1:_variable_num(_working_problem(m))
                 m._lower_solution[j][i] = MOI.get(d, MOI.VariablePrimal(j), m._relaxed_variable_index[i])
             end
         end
         m._lower_result_count = result_count
    elseif valid_flag && !feasible_flag
        m._lower_feasibility  = false
        m._lower_objective_value = -Inf
    else
        fallback_interval_lower_bound!(m)
    end

    return nothing
end

"""
$(SIGNATURES)

Checks if additional cuts should be added by the `add_cut!` routine.
"""
cut_condition(t::ExtensionType, m::GlobalOptimizer) = false

"""
$(SIGNATURES)

Add additionals cuts to subproblem while `cut_condition` is satisfied subproblem.
"""
add_cut!(t::ExtensionType, m::GlobalOptimizer) = nothing

#=
m._working_problem._relaxed_evaluator.is_post = m.subgrad_tighten
if !m._obbt_performed_flag
    if m._nonlinear_evaluator_created
        set_node!(m._working_problem._relaxed_evaluator, n)
        set_node_flag!(m)
        set_reference_point!(m)
        fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)
    end
    update_relaxed_problem_box!(m)
end
m._working_problem._relaxed_evaluator.interval_intersect = false

if !m._obbt_performed_flag
    relax_constraints!(m, 1)
end

# Optimizes the object
MOI.set(m.relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
MOI.optimize!(m.relaxed_optimizer)
if !intrepret_relaxed_solution(m, m.relaxed_optimizer)
    fallback_interval_lower_bound!(m, n)
end
=#