text:
## Unsafe scalar indexing

# Extract relevant fields of a DataArray to a tuple
# The extracted tuple can be passed to `unsafe_isna`,
# `unsafe_getindex_notna`, `unsafe_setna!`, `unsafe_setnotna!`, and
# `unsafe_dasetindex!`. This has a meaningful performance impact within
# very tight loops.
daextract(da::DataArray) = (da.data, da.na.chunks)
daextract(pda::PooledDataArray) = (pda.refs, pda.pool)
daextract(a) = nothing

# Check for NA
unsafe_isna(da::DataArray, extr, idx::Real) = Base.unsafe_bitgetindex(extr[2], idx)
unsafe_isna(pda::PooledDataArray, extr, idx::Real) = extr[1][idx] == 0
unsafe_isna(a, extr, idx::Real) = false
unsafe_getindex_notna(da::DataArray, extr, idx::Real) = getindex(extr[1], idx)
unsafe_getindex_notna(pda::PooledDataArray, extr, idx::Real) = getindex(extr[2], extr[1][idx])
unsafe_getindex_notna(a, extr, idx::Real) = Base.unsafe_getindex(a, idx)

# Set NA or data portion of DataArray

unsafe_bitsettrue!(chunks::Vector{UInt64}, idx::Real) =
    chunks[Base._div64(@compat(Int(idx))-1)+1] |= (@compat(UInt64(1)) << Base._mod64(@compat(Int(idx))-1))
unsafe_bitsetfalse!(chunks::Vector{UInt64}, idx::Real) =
    chunks[Base._div64(@compat(Int(idx))-1)+1] &= ~(@compat(UInt64(1)) << Base._mod64(@compat(Int(idx))-1))

unsafe_setna!(da::DataArray, extr, idx::Real) = unsafe_bitsettrue!(extr[2], idx)
unsafe_setna!(da::PooledDataArray, extr, idx::Real) = setindex!(extr[1], 0, idx)
unsafe_setnotna!(da::DataArray, extr, idx::Real) = unsafe_bitsetfalse!(extr[2], idx)
unsafe_setnotna!(da::PooledDataArray, extr, idx::Real) = nothing

# Fast setting of NA values in DataArrays
# These take the data and chunks (extracted as da.data and
# da.na.chunks), a value, and a linear index. They assume
# a certain initialization pattern:
#
# - For DataArrays, da.na should be falses
# - For PooledDataArrays, pda.refs should be zeros
unsafe_dasetindex!(data::Array, na_chunks::Vector{UInt64}, val::NAtype, idx::Real) =
    unsafe_bitsettrue!(na_chunks, idx)
unsafe_dasetindex!(data::Array, na_chunks::Vector{UInt64}, val, idx::Real) =
    setindex!(data, val, idx)
unsafe_dasetindex!(da::DataArray, extr, val::NAtype, idx::Real) =
    unsafe_setna!(da, extr, idx)
unsafe_dasetindex!(da::PooledDataArray, extr, val::NAtype, idx::Real) = nothing
unsafe_dasetindex!(da::DataArray, extr, val, idx::Real) = setindex!(extr[1], val, idx)
unsafe_dasetindex!(pda::PooledDataArray, extr, val, idx::Real) =
    setindex!(extr[1], getpoolidx(pda, val), idx)
unsafe_dasetindex!(a::AbstractArray, extr, val, idx::Real) = setindex!(a, val, idx)

## PooledDataArray helper functions

# Append newpool to pool. Return indices of newpool in pool.
function combine_pools!(pool, newpool)
    seen = Dict{eltype(pool),Int}()
    sizehint!(seen, length(pool)+length(newpool))

    # Create mapping from pool elements to indices
    i = 0
    for elem in pool
        seen[elem] = (i += 1)
    end

    # Find pool elements in existing array, or add them
    poolidx = Array(Int, length(newpool))
    for j = 1:length(newpool)
        poolidx[j] = Base.@get!(seen, newpool[j], (push!(pool, newpool[j]); i += 1))
    end
    poolidx
end

## General indexing functions

# Indexing with NA throws an error
function Base.to_index(A::DataArray)
    any(A.na) && throw(NAException("cannot index an array with a DataArray containing NA values"))
    Base.to_index(A.data)
end

if isdefined(Base, :checkindex) && isdefined(Base, :AbstractUnitRange)
    Base.checkindex(::Type{Bool}, ::AbstractUnitRange, ::NAtype) =
        throw(NAException("cannot index an array with a DataArray containing NA values"))
elseif isdefined(Base, :checkindex)
    Base.checkindex(::Type{Bool}, ::UnitRange, ::NAtype) =
        throw(NAException("cannot index an array with a DataArray containing NA values"))
else
    Base.checkbounds(::Type{Bool}, sz::Int, I::AbstractDataVector{Bool}) = length(I) == sz
    function Base.checkbounds{T<:Real}(::Type{Bool}, sz::Int, I::AbstractDataArray{T})
        anyna(I) && throw(NAException("cannot index into an array with a DataArray containing NAs"))
        extr = daextract(I)
        b = true
        for i = 1:length(I)
            @inbounds v = unsafe_getindex_notna(I, extr, i)
            b &= Base.checkbounds(Bool, sz, v)
        end
        b
    end
end

import Base: index_shape, index_lengths, setindex_shape_check

if isdefined(Base, :OneTo)
    _index_shape(x...) = Base.to_shape(index_shape(x...))
else
    _index_shape = index_shape
end

# Fallbacks to avoid ambiguity
Base.setindex!(t::AbstractDataArray, x, i::Real) =
    throw(MethodError(setindex!, typeof(t), typeof(x), typeof(i)))
Base.getindex(t::AbstractDataArray, i::Real) =
    throw(MethodError(getindex, typeof(t), typeof(i)))

## getindex: DataArray

# Scalar case
function Base.getindex(da::DataArray, I::Real)
    if getindex(da.na, I)
        return NA
    else
        return getindex(da.data, I)
    end
end
@nsplat N function Base.getindex(da::DataArray, I::NTuple{N,Real}...)
    if getindex(da.na, I...)
        return NA
    else
        return getindex(da.data, I...)
    end
end

if VERSION > v"0.5-"
    Base.unsafe_getindex(x::Number, i::Int) = (@inbounds r = x[i]; r)
end

# Vector case
@generated function _getindex!(dest::DataArray, src::DataArray, I::Union{Real, AbstractArray, Colon}...)
    N = length(I)
    quote
        $(Expr(:meta, :inline))
        idxlens = index_lengths(src, I...) # TODO: unsplat?
        srcextr = daextract(src)
        destextr = daextract(dest)
        srcsz = size(src)
        k = 1
        @nloops $N i d->(1:idxlens[d]) d->(@inbounds j_d = getindex(I[d], i_d)) begin
            offset_0 = @ncall $N sub2ind srcsz j
            if unsafe_isna(src, srcextr, offset_0)
                unsafe_dasetindex!(dest, destextr, NA, k)
            else
                unsafe_dasetindex!(dest, destextr, unsafe_getindex_notna(src, srcextr, offset_0), k)
            end
            k += 1
        end
        dest
    end
end

function _getindex{T}(A::DataArray{T}, I::@compat Tuple{Vararg{Union{Int,AbstractVector}}})
    shape = _index_shape(A, I...)
    _getindex!(DataArray(Array(T, shape), falses(shape)), A, I...)
end

@nsplat N function Base.getindex(A::DataArray, I::NTuple{N,(@compat Union{Real,AbstractVector})}...)
    checkbounds(A, I...)
    _getindex(A, Base.to_indexes(I...))
end

# Dispatch our implementation for these cases instead of Base
function Base.getindex(A::DataArray, I::AbstractVector)
    checkbounds(A, I)
    _getindex(A, (Base.to_index(I),))
end
function Base.getindex(A::DataArray, I::AbstractArray)
    checkbounds(A, I)
    _getindex(A, (Base.to_index(I),))
end

## getindex: PooledDataArray

# Scalar case
function Base.getindex(pda::PooledDataArray, I::Real)
    if getindex(pda.refs, I) == 0
        return NA
    else
        return pda.pool[getindex(pda.refs, I)]
    end
end
@nsplat N function Base.getindex(pda::PooledDataArray, I::NTuple{N,Real}...)
    if getindex(pda.refs, I...) == 0
        return NA
    else
        return pda.pool[getindex(pda.refs, I...)]
    end
end

# Vector case
@nsplat N function Base.getindex(A::PooledDataArray, I::NTuple{N,(@compat Union{Real,AbstractVector})}...)
    PooledDataArray(RefArray(getindex(A.refs, I...)), copy(A.pool))
end

# Dispatch our implementation for these cases instead of Base
Base.getindex(A::PooledDataArray, I::AbstractVector) =
    PooledDataArray(RefArray(getindex(A.refs, I)), copy(A.pool))
Base.getindex(A::PooledDataArray, I::AbstractArray) =
    PooledDataArray(RefArray(getindex(A.refs, I)), copy(A.pool))

## setindex!: DataArray

function Base.setindex!(da::DataArray, val::NAtype, i::Real)
    da.na[i] = true
    return da
end

function Base.setindex!(da::DataArray, val, ind::Real)
    da.data[ind] = val
    da.na[ind] = false
    return da
end

## setindex!: PooledDataArray

function Base.setindex!(pda::PooledDataArray, val::NAtype, ind::Real)
    pda.refs[ind] = 0
    return pda
end

function Base.setindex!(x::PooledDataArray, val, ind::Real)
    x.refs[ind] = getpoolidx(x, val)
    return x
end

## setindex!: both DataArray and PooledDataArray

@ngenerate N typeof(A) function Base.setindex!(A::AbstractDataArray, x,
                                               J::NTuple{N,(@compat Union{Real,AbstractArray})}...)
    if !isa(x, AbstractArray) && isa(A, PooledDataArray)
        # Only perform one pool lookup when assigning a scalar value in
        # a PooledDataArray
        setindex!(A.refs, getpoolidx(A, x), J...)
        return A
    end

    Aextr = daextract(A)
    @ncall N checkbounds A J
    @nexprs N d->(I_d = Base.to_index(J_d))
    stride_1 = 1
    @nexprs N d->(stride_{d+1} = stride_d*size(A,d))
    @nexprs N d->(offset_d = 1)  # really only need offset_$N = 1
    if !isa(x, AbstractArray)
        @nloops N i d->(1:length(I_d)) d->(@inbounds offset_{d-1} = offset_d + (Base.unsafe_getindex(I_d, i_d)-1)*stride_d) begin
            if isa(x, NAtype)
                @inbounds unsafe_setna!(A, Aextr, offset_0)
            else
                @inbounds unsafe_setnotna!(A, Aextr, offset_0)
                @inbounds unsafe_dasetindex!(A, Aextr, x, offset_0)
            end
        end
    else
        X = x
        idxlens = @ncall N index_lengths A I
        @ncall N setindex_shape_check X (d->idxlens[d])
        k = 1
        if isa(A, PooledDataArray) && isa(X, PooledDataArray)
            # When putting one PDA into another, first unify the pools
            # and then translate the references
            poolmap = combine_pools!(A.pool, X.pool)
            Arefs = A.refs
            Xrefs = X.refs
            @nloops N i d->(1:idxlens[d]) d->(@inbounds offset_{d-1} = offset_d + (Base.unsafe_getindex(I_d, i_d)-1)*stride_d) begin
                @inbounds Arefs[offset_0] = Xrefs[k] == 0 ? 0 : poolmap[Xrefs[k]]
                k += 1
            end
        else
            Xextr = daextract(X)
            @nloops N i d->(1:idxlens[d]) d->(@inbounds offset_{d-1} = offset_d + (Base.unsafe_getindex(I_d, i_d)-1)*stride_d) begin
                @inbounds if isa(X, AbstractDataArray) && unsafe_isna(X, Xextr, k)
                    unsafe_setna!(A, Aextr, offset_0)
                else
                    unsafe_setnotna!(A, Aextr, offset_0)
                    unsafe_dasetindex!(A, Aextr, unsafe_getindex_notna(X, Xextr, k), offset_0)
                end
                k += 1
            end
        end
    end
    A
end


text:
s = open("input.txt") do file
    readlines(file);
end

t = open("input2.txt") do file
    readlines(file);
end

function ParseRules(lines)
    rules = Dict();

    for line in lines
        if length(line) < 2
            break;
        end

        parts = split(line, ": ");
        id = parse(Int, parts[1]);
        options = []

        if occursin("|", parts[2])
            ruleBits = split(parts[2]," | ")
        else
            ruleBits = [parts[2]]
        end

        m = match(r"\"([ab])\"", line);
        if m != nothing
            # This is a base rule
            push!(options, String(m.captures[1]));
        else
            println(ruleBits);
            options = map(x->map(y->parse(Int,y), split(x)), ruleBits);
        end

        rules[id] = options;
    end

    return rules;
end

function CheckRule(id, msg, rules)
    r = rules[id];

    if length(r) == 1 && typeof(r[1]) == String
        # base rule of 1 char
        rule = Regex("^"*solved[id]);
        m = match(rule, msg);

        if m == nothing
            return [];
        else
            # println("$(msg) : $(solved[id])");
            return [length(m.match)];
        end
    else
        opts = [];
        for opt in r
            totalLen = [0]
            for rule in opt
                ruleLen = [];

                for tl in totalLen
                    len = CheckRule(rule, msg[(tl+1):end], rules);

                    if length(len) > 0
                        ruleLen = vcat(ruleLen, len.+tl);
                    end
                end

                totalLen = ruleLen;

                if length(ruleLen) == 0
                    break
                end
            end

            if length(totalLen) > 0
                opts = vcat(opts, totalLen);
            end
        end

        return opts;
    end
end

function GetCount(messages, rules, solved)
    count = 0;
    matched = [];
    for msg in messages
        if any(CheckRule(0, msg, rules) .== length(msg))
            count += 1;
            push!(matched, true)
        else
            push!(matched, false)
        end
    end

    println(">>>> $(count)");
    return matched;
end

rules = ParseRules(s)

rules[11] = [[42, 11, 31], [42, 31]];
rules[8] = [[42], [42, 8]];

messages = s[length(rules)+2:end]
uncorrupt = GetCount(messages, rules, solved)

# test
rules_t = ParseRules(t)

rules[11] = [[42, 11, 31], [42, 31]];
rules[8] = [[42], [42, 8]];

messages_t = t[length(rules_t)+2:end]
uncorrupt_t = GetCount(messages_t, rules_t, solved_t)

ext: jl


text:
using LightGraphsIO
using Test

#@testset "LightGraphsIO" begin
    include("internalstests.jl")
    include("DOTtests.jl")
    include("Edgelisttests.jl")
    include("GMLtests.jl")
    include("NETtests.jl")
    include("GEXFtests.jl")
    include("Graph6tests.jl")
    include("GraphMLtests.jl")
    include("CDFtests.jl")
    include("LGCompressedtests.jl")
#end


text:
using StatsBase
using Base.Test

# rle & inverse_rle

z = [1, 1, 2, 2, 2, 3, 1, 2, 2, 3, 3, 3, 3]
(vals, lens) = rle(z)
@test vals == [1, 2, 3, 1, 2, 3]
@test lens == [2, 3, 1, 1, 2, 4]
@test inverse_rle(vals, lens) == z

z = [true, true, false, false, true, false, true, true, true]
vals, lens = rle(z)
@test vals == [true, false, true, false, true]
@test lens == [2, 2, 1, 1, 3]
@test inverse_rle(vals, lens) == z

# levelsmap
a = [1, 1, 2, 2, 2, 3, 1, 2, 2, 3, 3, 3, 3, 2]
b = [true, false, false, true, false, true, true, false]

@test levelsmap(a) == Dict(2=>2, 3=>3, 1=>1)
@test levelsmap(b) == Dict(false=>2, true=>1)

# indicatormat

I = [false true  false false false;
     true  false false false true;
     false false true  true  false]

x = [2, 1, 3, 3, 2]
@test indicatormat(x, 3) == I
@test full(indicatormat(x, 3; sparse=true)) == I

x = ["b", "a", "c", "c", "b"]
@test indicatormat(x) == I
@test full(indicatormat(x; sparse=true)) == I


text:
@testset "Normal form f4 finite fields" begin

    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x", "y", "z"])

    # Regression test, Structural Identifiability
    @test Groebner.normalform([x], R(0)) == R(0)
    @test Groebner.normalform([x], R(1)) == R(1)

    Gs = [
        [x], [x, y], [x, y, z]
    ]

    @test Groebner.normalform(Gs[1], x) == 0
    @test Groebner.normalform(Gs[1], y) == y
    @test Groebner.normalform(Gs[3], x + y^2) == 0
    @test Groebner.normalform(Gs[2], 5x^2 - y^3 - 1) == -1
    @test Groebner.normalform(Gs[2], x + y + 3z) == 3z

    G = [
        x^2 + y*x + 1,
        y^2 - y
    ]

    @test Groebner.normalform(G, x + y + z) == x + y + z
    @test Groebner.normalform(G, x^2 + y^2) == y - y*x - 1
    @test Groebner.normalform(G, y^3) == Groebner.normalform(G, y^4) == y

    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x","y","z"], ordering=:degrevlex)
    G = [
        x^2 + y,
        y^2 + x
    ]
    @test Groebner.normalform(G, x^2 + y^2) == -x - y

    @test_nowarn Groebner.normalform([x, x*y], x)
end

@testset "Normal form -- a lot of variables" begin
    R, xs = PolynomialRing(QQ, ["x$i" for i in 1:63])

    GB = Groebner.groebner(xs)
    @test GB == reverse(xs)
    @test Groebner.normalform(GB, xs[1]) == R(0)

    R, xs = PolynomialRing(QQ, ["x$i" for i in 1:220])

    GB = Groebner.groebner(xs)
    @test GB == reverse(xs)
    @test Groebner.normalform(GB, xs[1]) == R(0)

    #=
    R, xs = PolynomialRing(QQ, ["x$i" for i in 1:1049])

    GB = Groebner.groebner(xs)
    @test GB == reverse(xs)
    @test Groebner.normalform(GB, xs[1]) == R(0)
    =#

end


text:
module WGL

using ..Glimmer
import ..UIVariables
using Parameters
import UUIDs
import Colors;
using OrderedCollections
import JSON

export  WGLColor,
        WGLWireframe,
        WGLMaterial,
        WGLMeshGeometry,
        WGLLinesGeometry,
        WGLPointsGeometry,
        WGLObject,
        WGLTransform,
        WGLSetTransform,
        wglSet!,

        DymmyExport

function typedict(x::AbstractWGLComponent)
    res = OrderedDict()

    for fn in filter(x -> !startswith(string(x), "_"), fieldnames(typeof(x)))
        key = fn
        val = getfield(x, fn)
        if (val === nothing)
            continue
        end
        if (isstructtype(typeof(val)) && !(val isa String))
            if (val isa Dict || val isa OrderedDict)
                # @info "DICT" key, val
                res[key] = val;
            elseif (val isa Vector)
                res[key] = val;
            else
                # @show key
                res[key] = typedict(val)
            end
        else
            res[key] = val;
        end
    end
    # if (isstructtype(typeof(x)))
    # Dict(fn=>getfield(x, fn) for fn ∈ filter(x->!startswith(string(x), "_"), fieldnames(typeof(x))))

    return res;
end

# -----------------------------------------------------------------------

@enum WGLTypes begin
    TMaterial = 1
    TGeometry = 2
    TObject = 3
    TOther = 4
end

@with_kw mutable struct WGLPath <: AbstractWGLComponent
    var::AbstractUIVariable
    path::String
    type::WGLTypes
end

@with_kw mutable struct WGLColor <: AbstractWGLComponent
    r::Float32 = 0.0f
    g::Float32 = 0.0f
    b::Float32 = 0.0f
    a::Float32 = 1.0f
end

WGLColor(c::Tuple{Int64,Int64,Int64}) = WGLColor(
    r=Float32(c[1]) / 255.0f0,
    g=Float32(c[2]) / 255.0f0,
    b=Float32(c[3]) / 255.0f0,
    a=1.0f0)

@with_kw mutable struct WGLWireframe <: AbstractWGLComponent
    show::Bool = false
    lineWidth::Float32 = 1.0f
    color::WGLColor = WGLColor((0, 0, 0))
end

WGLWireframe(lineWidth::Float32, color::WGLColor=WGLColor(0.0f0, 0.0f0, 0.0f0, 1.0f0)) = WGLWireframe(
    show=true,
    lineWidth=lineWidth,
    color=color
)

@with_kw mutable struct WGLTransform <: AbstractWGLComponent
    matrixCoef::Union{Nothing, Vector{Float32}} = nothing

    scale::Union{Nothing, Vector{Float32}} = nothing
    translation::Union{Nothing, Vector{Float32}} = nothing
    rotation::Union{Nothing, Vector{Float32}} = nothing
end

WGLTransform(pos::Tuple{Float32,Float32,Float32}) = WGLTransform(
    translation = [pos[1], pos[2], pos[3]],
    rotation = [0.0f0, 0.0f0, 0.0f0],
    scale = [1.0f0, 1.0f0, 1.0f0],
)

@with_kw mutable struct WGLMaterial <: AbstractWGLMaterial
    name::String
    type::String
    color::WGLColor = WGLColor((0, 0, 0))
    wireframe::Union{Nothing, WGLWireframe} = WGLWireframe(show=false)
end

@with_kw mutable struct WGLMeshGeometry <: AbstractWGLMaterial
    name::String
    vertices::Union{Nothing, Vector{Float32}} = nothing
    triangles::Union{Nothing, Vector{Int64}} = nothing
    uvs::Union{Nothing, Vector{Float32}} = nothing
end

@with_kw mutable struct WGLLinesGeometry <: AbstractWGLMaterial
    name::String
    vertices::Union{Nothing, Vector{Float32}} = nothing
    continuous::Bool = false
    joins::String = "none"                          # none, fill, weld
end

@with_kw mutable struct WGLPointsGeometry <: AbstractWGLMaterial
    name::String
    vertices::Union{Nothing, Vector{Float32}} = nothing
end

@with_kw mutable struct WGLObject <: AbstractWGLMaterial
    name::String
    parentName::String = ""
    transform::WGLTransform = WGLTransform((0f0, 0f0, 0f0))
    geometryName::String
    materialName::String
end

@with_kw mutable struct WGLSetTransform <: AbstractWGLMaterial
    name::String = ""
    transform::WGLTransform = WGLTransform((0f0, 0f0, 0f0))
end

# ---------------------------------------------------------------------------

function wglSendToRenderer!(var::AbstractUIVariable, obj::AbstractWGLComponent, methodName::String; mergeData::Union{Nothing, Dict, OrderedDict} = nothing)
    data = WGL.typedict(obj)

    if (mergeData !== nothing)
        data = merge(data, mergeData)
    end

    json = JSON.json(data)
    js_data = Dict{Symbol, Any}(
        :methodName => methodName,
        :json => json,
    )
    send!(var, js_data);
end

function wglSet!(var::AbstractUIVariable, mat::WGLMaterial)::WGLPath
    wglSendToRenderer!(var, mat, "ProcessMaterial")

    return WGLPath(
        var = var,
        path = mat.name,
        type = TMaterial
    )
end

function wglSet!(var::AbstractUIVariable, geometry::WGLMeshGeometry)::WGLPath
    wglSendToRenderer!(var, geometry, "ProcessTriangleMesh")

    return WGLPath(
        var = var,
        path = geometry.name,
        type = TGeometry
    )
end

function wglSet!(var::AbstractUIVariable, geometry::WGLLinesGeometry)::WGLPath
    wglSendToRenderer!(var, geometry, "ProcessLines")

    return WGLPath(
        var = var,
        path = geometry.name,
        type = TGeometry
    )
end

function wglSet!(var::AbstractUIVariable, geometry::WGLPointsGeometry)::WGLPath
    wglSendToRenderer!(var, geometry, "ProcessPoints")

    return WGLPath(
        var = var,
        path = geometry.name,
        type = TGeometry
    )
end

function wglSet!(var::AbstractUIVariable, object::WGLObject)::WGLPath
    wglSendToRenderer!(var, object, "ProcessObject")

    object_path = object.parentName * "/" * object.name;

    return WGLPath(
        var = var,
        path = object_path,
        type = TObject
    )
end

function wglSet!(path::WGLPath, object::WGLObject)::WGLPath
    wglSendToRenderer!(path.var, object, "ProcessObject"; mergeData = Dict(:parentName => path.path))

    object_path = path.path * "/" * object.name;

    return WGLPath(
        var = path.var,
        path = object_path,
        type = TObject
    )
end

function wglSet!(var::AbstractUIVariable, objectTransform::WGLSetTransform)::WGLPath
    wglSendToRenderer!(var, objectTransform, "ProcessObjectTransform")

    return WGLPath(
        var = var,
        path = objectTransform.name,
        type = TOther
    )
end

function wglSet!(path::WGLPath, objectTransform::WGLSetTransform)::WGLPath
    wglSendToRenderer!(path.var, objectTransform, "ProcessObjectTransform"; mergeData = Dict(:name => path.path))

    return WGLPath(
        var = path.var,
        path = path.path,
        type = TOther
    )
end

end # module WGL


text:
module OCTHaGOn
    using Combinatorics
    using CPLEX
    using DataFrames
    using DocStringExtensions
    using ForwardDiff
    using JuMP
    using LatinHypercubeSampling
    using LinearAlgebra
    using MacroTools
    using Main.IAI
    using MathOptInterface
    using MathOptSetDistances
    using Missings
    using NearestNeighbors
    using Parameters
    using ProgressMeter

    const MOI = MathOptInterface
    const VALID_OPTIMIZERS = [:Cbc, :CPLEX, :Gurobi, :GLPK, :Mosek, :SCIP]
    const SOLVER = CPLEX.Optimizer
    const SOLVER_SILENT = with_optimizer(SOLVER, CPX_PARAM_SCRIND = 0)

    const PROJECT_ROOT = dirname(dirname(@__FILE__))
    const DATA_DIR = PROJECT_ROOT * "\\data\\"
    const BARON_DIR = PROJECT_ROOT * "\\data\\baron\\"
    const GAMS_DIR = PROJECT_ROOT * "\\data\\gams\\"
    const SOL_DIR = PROJECT_ROOT * "\\data\\solutions\\"
    const TREE_DIR = PROJECT_ROOT * "\\data\\trees\\"

    const valid_lowers = ["reg", "lower", "upperlower"]
    const valid_uppers = ["reg", "upper", "upperlower"]
    const valid_singles = ["reg", "upperlower"]
    const valid_pairs =  ["lower" => "upper",
                          "upper" => "lower",
                          "reg" => "upper", # reg with some threshold.second
                          "upper" => "reg"] # reg with some threshold.second

    include("small_scripts.jl")

    include("on_jump.jl")

    include("learners.jl")

    include("exceptions.jl")

    include("settings.jl")

    include("black_box_learners.jl")

    include("root_finding.jl")

    include("iai_wrappers.jl")

    include("global_model.jl")

    include("data_driven.jl")

    include("sample.jl")

    include("constraintify.jl")

    include("fit.jl")

    include("refine.jl")

    include("optimize.jl")

    include("plot.jl")

        # STRUCTS
    export GlobalModel, BlackBoxLearner, BlackBoxClassifier, BlackBoxRegressor,
        LinkedLearner, LinkedClassifier, LinkedRegressor,

    # JuMP.Model extensions to GlobalModel
        set_optimizer, optimize!,

        # PROBLEM CONSTRUCTION

        # Explicit constraints
        add_nonlinear_constraint, add_nonlinear_or_compatible,
        add_linked_constraint,

        # Data driven constraints
        add_variables_from_data!, add_datadriven_constraint, bound_to_data!,

        # CONSTRAINT LEARNING

        # Bound finding
        find_bounds!, find_linear_bounds!,

        # Sampling methods (called on BBLs)
        lh_sample, boundary_sample, last_leaf_sample,
        feasibility_sample, upper_bound_sample,
        uniform_sample_and_eval!, knn_sample,
        build_knn_tree,
        find_knn, classify_patches,

        # Evaluation
        evaluate, evaluate_gradient, update_gradients,
        classify_curvature, eval!, add_data!,

        # Convexity checks
        update_local_convexity, update_vexity, update_leaf_vexity,
        solution, evaluate_accuracy,
        feas_gap, print_feas_gaps,
        nonlinearize!,
        set_param, get_param,

        # Tree generation and learning
        learn_constraint!,
        base_regressor, base_classifier, base_rf_regressor, base_rf_classifier,
        fit_regressor_kwargs, fit_classifier_kwargs,
        regressor_kwargs, classifier_kwargs,

        # Constraint/variable generation/deletion
        set_bigM,
        add_feas_constraints!, add_regr_constraints!,
        add_tree_constraints!, clear_tree_constraints!, update_tree_constraints!,
        clear_lower_constraints!, clear_upper_constraints!,
        add_relaxation_variables!, clear_relaxation_variables!,
        detect_linked_constraints!, merge_linked_constraints!,

        # BlackBoxLearner helpers
        show_trees, all_mi_constraints, active_lower_tree, active_upper_tree,
        active_leaves, find_leaves, get_random_trees,

        # Algorithms
        add_infeasibility_cuts!, boundify, surveysolve,
        relax_objective!, tight_objective!,
        globalsolve!, globalsolve_and_time!, descend!,
        compute_hyperplane_bigM_upper, compute_hyperplane_bigM_lower,
        compute_regression_bigM_lower, compute_regression_bigM_upper,

        # Checks and exceptions
        feasibility, check_accuracy, check_feasibility,
        check_bounds,  check_sampled, check_if_trained,
        OCTHaGOnException, is_feasible, is_sampled,

        # Clean-up
        clear_data!, clear_tree_data!,

        # HELPERS

        # IAI helpers,
        pwl_constraint_data, trust_region_data,

        # Functions on JuMP/MathOptInterface objects
        fetch_variable, get_bounds, get_joint_bounds,
        get_unbounds, flattened_bounds,
        bounded_aux,
        linearize_objective!, classify_constraints,
        bound!, restrict_to_set,
        distance_to_set, get_constant,
        count_constraints, count_variables, count_types,

        # Data manipulation
        data_to_DataFrame, data_to_Dict,

        # Expr manipulation
        functionify, gradientify,

        # Other approximators,
        ridge_regress, u_regress, l_regress, svm, reweight,

        # Functions to import global optimization problems,
        sagemark_to_GlobalModel,
        alphac_to_expr, alphac_to_objexpr, alphac_to_varbound!,

        # Small scripts
        deconstruct, flat, get_datamap, get_var_ranges, get_varmap,
        power, sqr, merge_kwargs, normalized_data, print_details, substitute, vars_from_expr, vars_from_constraint
    end


text:
fn0 = (q, z) -> begin
    m = trace(:m, Normal(q, 1.0))
    q = trace(:q, Normal(z, 5.0))
    m + q
end

fn1 = (q, z, m) -> begin
    t = trace(:z, Normal(5.0, 1.0))
    n = trace(:m, Normal(z, 1.0))
    l = trace(:l, fn0, q, m)
    if l > 5.0
        trace(:q, Normal(1.0, 5.0))
    else
        trace(:q, Normal(5.0, 3.0))
    end
    p = trace(:p, Normal(z, 3.0))
    t = trace(:t, Normal(p, 5.0))
    l + n
end

@testset begin "Dynamic specialization update test 1"
    ret, cl = simulate(fn1, 10.0, 10.0, 1.0)
    prev_score = get_score(cl)

    ret, cl, w, _ = update(Jaynes.Empty(), cl,
                           Δ(5.0, ScalarDiff(-5.0)),
                           Δ(5.0, ScalarDiff(-5.0)),
                           Δ(1.0, NoChange()))
    @test get_score(cl) - w ≈ prev_score
end


text:
using GPMTD
using Distributed

# using BSON: @save, @load
using BSON

using Random
using PDMats
using LinearAlgebra
BLAS.set_num_threads(1)

using Statistics
using Dates

const srnd = parse(Int64, ARGS[1])
const n = parse(Int64, ARGS[2])
const L = parse(Int64, ARGS[3])
const whichdat = ARGS[4]

## Keyword args
# srnd = 427191
# n = 50
# L = 2

## Settings in script
n_burn = Int(50e3)
n_keep = 1000
n_adapt = 2000
n_thin = 20
centerscale = false

## Load data
# whichdat = "Ricker_lag2only"

BSON.@load "data/$(whichdat).bson" y X

y = deepcopy(y[1:n])
X = deepcopy(X[1:n,1:L])

## Setup
mesg = "$(whichdat)_n$(n)_L$(L)_$(srnd)"
Random.seed!(srnd)

mean(y)
std(y)

if centerscale
    X = (X .- mean(y)) ./ std(y)
    y = (y .- mean(y)) ./ std(y)
end

## Prior
prior = Prior_GPMTD(L, InterceptNormal(), MixComponentNormal())
# for j = 1:L # "deterministic"
#     prior.mixcomps[j].σ2.ν = float(n)
#     prior.mixcomps[j].σ2.s0 = 0.0001
# end

## Init
init = State_GPMTD(X, InterceptNormal(), MixComponentNormal())

## log file
logfilename = "postsim_progress/out_prog_$(mesg).txt"

## initialize model
model = Model_GPMTD(y, X, prior, init)

monitor = [:intercept, :lλ, :ζ, :μ, :σ2, :κ, :corParams, :κ_hypers, :corHypers, :fx]
# monitor = [:intercept, :lλ, :ζ :μ, :σ2, :κ, :corParams, :κ_hypers, :corHypers, :fx]

updatevars = [:intercept, :lλ, :ζ, :μ, :σ2, :κ, :corParams, :κ_hypers, :corHypers, :fx]
# updatevars = [:intercept, :lλ, :ζ, :μ, :σ2, :κ, :corParams, :κ_hypers, :corHypers, :fx]

## set up cores for distributed updates of mixcomps
# n_workers = min(length(Sys.cpu_info()) - 1, L)
n_workers = 0
addprocs(n_workers)
n_procs = nprocs()
# n_procs = 1
workers()
@everywhere using Random123, GPMTD

## runs
@time iter, accptr = mcmc!(model, 1000, save=false, n_procs=n_procs,
        report_filename=logfilename, report_freq=100, update=updatevars)
println(model.state.iter)
println(accptr)

adapt!(model, n_iter_collectSS=n_adapt, n_iter_scale=n_adapt,
    adjust_bnds=[0.01, 5.0], maxtries=10;
    report_filename=logfilename, update=updatevars)

## burn-in
timestartburn = Dates.now()
iter, accptr = mcmc!(model, n_burn, thin=1, save=false, n_procs=n_procs,
    report_filename=logfilename, report_freq=10000, update=updatevars, tol_posdef=tol_posdef)
etr(timestartburn, n_burn, n_keep, n_thin, logfilename)

sims, accptr = mcmc!(model, n_keep, save=true, thin=n_thin,
    n_procs=n_procs,
    report_filename=logfilename, report_freq=10000,
    update=updatevars, monitor=monitor)

if n_workers > 0
    rmprocs(workers())
end

using BSON
bson("postsim/postsim_$(mesg).bson",
    model=deepcopy(model), sims=deepcopy(sims),
    init=deepcopy(init), mesg=deepcopy(mesg) )


text:
mutable struct Box
    upper_left_x::Number
    upper_left_y::Number
    lower_right_x::Number
    lower_right_y::Number
end

"""
    Base.:(==)(a::Box, b::Box)::Bool

Tests for equivalence between two Box instances.
"""
function Base.:(==)(a::Box, b::Box)::Bool
    return  (a.upper_left_x == b.upper_left_x) &&
            (a.upper_left_y == b.upper_left_y) &&
            (a.lower_right_x == b.lower_right_x) &&
            (a.lower_right_y == b.lower_right_y)
end

"""
    barycenterbox( x::Number, y::Number, height::Number, width::Number)

Converts a detection box from x,y,h,w coordinates to upper left lower right
coordinates.
"""
function barycenterbox( x::Number, y::Number, height::Number, width::Number)
    wh, hh = width/2, height/2
    return Box( x - wh, y - hh, x + wh, y + hh )
end

"""
    LLURtoULLR( LL_x::Number, LL_y::Number, UR_x::Number, UR_y::Number)

Converts a detection box from lowerleft-upperright coordinates to upper left lower right
coordinates.

Untested...
"""
function LLURtoULLR( LL_x::Number, LL_y::Number, UR_x::Number, UR_y::Number)
    return Box( LL_x, UR_y, UR_x, LL_y )
end

"""
    area( a::Box )

Calculates the area of a Box instance.

"""
area( a::Box ) = ( a.lower_right_x - a.upper_left_x ) * ( a.lower_right_y - a.upper_left_y )

"""
    translate!(a::Box, x::Number,y::Number)::Box

Translates a box by the input x & y coordinates. This is an inplace operation.

"""
function translate!(a::Box, x::Number = 0, y::Number = 0)::Box
    a.upper_left_x  += x
    a.upper_left_y  += y
    a.lower_right_x += x
    a.lower_right_y += y
    return a
end

"""
    translate(a::Box, x::Number,y::Number)::Box

Translates a box by the input x & y coordinates.
"""
translate(a::Box, x::Number = 0, y::Number = 0)::Box = translate!( deepcopy( a ), x, y )

"""
    clamp!(a::Box, bound::Box)::Box

Clamps a `Box` instance `a` to be within a bounded `Box` region `bound`. This is an inplace operation.

"""
function clamp!(a::Box, bound::Box)::Box
    a.upper_left_x  = (a.upper_left_x < bound.upper_left_x) ? bound.upper_left_x : a.upper_left_x
    a.upper_left_y  = (a.upper_left_y < bound.upper_left_y) ? bound.upper_left_y : a.upper_left_y
    a.lower_right_x = (a.lower_right_x > bound.lower_right_x) ? bound.lower_right_x : a.lower_right_x
    a.lower_right_y = (a.lower_right_y > bound.lower_right_y) ? bound.lower_right_y : a.lower_right_y
    return a
end

"""
    clamp(a::Box, bound::Box)::Box

Clamps a `Box` instance `a` to be within a bounded `Box` region `bound`.

"""
clamp(a::Box, bound::Box)::Box = clamp!( deepcopy(a), bound )

"""
    intersection_area( a::Box, b::Box )::Number

Calculates the area of intersection between 2 Box instances.

"""
function intersection_area( a::Box, b::Box )::Number
	ul_x = max( a.upper_left_x,  b.upper_left_x )
	ul_y = max( a.upper_left_y,  b.upper_left_y )
	lr_x = min( a.lower_right_x, b.lower_right_x )
	lr_y = min( a.lower_right_y, b.lower_right_y )
	return max( 0, lr_x - ul_x ) * max( 0, lr_y - ul_y )
end

"""
    intersection_over_union( a::Box, b::Box )::Number

Calculates the intersection over union between 2 Box instances.

"""
function intersection_over_union( a::Box, b::Box )::Number
    intersection = intersection_area( a, b )
    return intersection / ( area( a ) + area( b ) - intersection)
end


text:
abstract type AbstractModel end

struct PoissonLapseUniform <: AbstractModel
    params::Vector{Float64}
    max::Int
end
function init(::Type{<:PoissonLapseUniform}, cm::Dict)
    m = maximum(key.Omissions_plus_one for key in keys(cm))
    return PoissonLapseUniform([10, 2, 0.1], m)
end

Distributions.params(p::PoissonLapseUniform) = p.params
Base.maximum(m::PoissonLapseUniform) = m.max

function Distributions.cdf(p::PoissonLapseUniform, x::Int)
    T, r, ϵ = params(p)
    d = Gamma(T, 1/r)
    # ϵ = 0
    # 1-ϵ e faccio Gamma, ϵ e faccio uniform, sommare e' uguale a "oppure" in probabilita'
    return (1-ϵ)*cdf(d, x) + ϵ*x/p.max
end
##
struct PoissonLapseExponential <: AbstractModel
    params::Vector{Float64}
    max::Int
end
function init(::Type{<:PoissonLapseExponential}, cm::Dict)
    m = maximum(key.Omissions_plus_one for key in keys(cm))
    return PoissonLapseExponential([10, 2, 3], m)
end

Distributions.params(p::PoissonLapseExponential) = p.params
Base.maximum(m::PoissonLapseExponential) = m.max

function Distributions.cdf(p::PoissonLapseExponential, x::Int)
    T, r, ϵ = params(p)
    g = Gamma(T, 1/r)
    l = Exponential(ϵ)
    #ϵ = 0
    # 1-ϵ e faccio Gamma, ϵ e faccio uniform, sommare e' uguale a "oppure" in probabilita'
    #return (1-cdf(l,x))*cdf(d, x) + cdf(l,x)
    return (1-0.01)*cdf(g,x) + 0.01*cdf(l,x)
end


text:
using ConventionalApp

proj = Project(@__DIR__)
setup_env(proj)
@include_startup(proj)

#show(activeproject)

using Revise
using OnDieDecap
#Main.include("src/OnDieDecap.jl")
#Main.OnDieDecap.run_app()


text:
import QuantumGlue

import QuantumOptics
const QO = QuantumOptics

import Yao
using Test

@testset "Convert multi-qubit states" begin
    for f in (Yao.zero_state, Yao.rand_state)
        s = f(3)
        @test s |> QO.Ket |> Yao.register == s
    end
end


text:
using SolidStateDetectors
using IntervalSets
using Plots; pyplot(fmt = :png);

T=Float32

struct CustomChargeDensity{T, N} <: SolidStateDetectors.AbstractChargeDensity{T}
       G1::SolidStateDetectors.AbstractGeometry{T,N}
       G2::SolidStateDetectors.AbstractGeometry{T,N}
       G3::SolidStateDetectors.AbstractGeometry{T,N}
       G4::SolidStateDetectors.AbstractGeometry{T,N}
       G5::SolidStateDetectors.AbstractGeometry{T,N}
       L1::SolidStateDetectors.AbstractChargeDensity{T}
       L2::SolidStateDetectors.AbstractChargeDensity{T}
       L3::SolidStateDetectors.AbstractChargeDensity{T}
       L4::SolidStateDetectors.AbstractChargeDensity{T}
       L5::SolidStateDetectors.AbstractChargeDensity{T}
end

function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::SolidStateDetectors.AbstractCoordinatePoint{T})::T where {T}
	(pt in cdm.G1) * SolidStateDetectors.get_charge_density(cdm.L1, pt) +
	(pt in cdm.G2) * SolidStateDetectors.get_charge_density(cdm.L2, pt) +
	(pt in cdm.G3) * SolidStateDetectors.get_charge_density(cdm.L3, pt) +
	(pt in cdm.G4) * SolidStateDetectors.get_charge_density(cdm.L4, pt) +
	(pt in cdm.G5) * SolidStateDetectors.get_charge_density(cdm.L5, pt)
end

G1 = SolidStateDetectors.Tube{T}(0..0.00725,0..6.28318530718,0..0.002,missing)
G2 = SolidStateDetectors.Tube{T}(0..0.00725,0..6.28318530718,0..0.00001,missing)
G3 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515,0.00001,CartesianPoint{T}(0,0,0.001995),0)
#pstop
G4 = SolidStateDetectors.HexagonalPrism{T}(0.00524, 0.00525,0.00001,CartesianPoint{T}(0,0,0.001995),0)
#adding gradient to p+ side
G5 = SolidStateDetectors.Tube{T}(0..0.00725,0..6.28318530718,0..0.00005,CartesianPoint{T}(0,0,0.000010000001))
CD5 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,-1.9f16),(0.0f0,0.0f0,0.3f21))

CD1 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,5.0f10),(0.0f0,0.0f0,0.0f0))
CD2 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,-1.9f16),(0.0f0,0.0f0,0.0f0))
CD3 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9f16),(0.0f0,0.0f0,0.0f0))
CD4 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,-1.9f16),(0.0f0,0.0f0,0.0f0))

ccdm = CustomChargeDensity(G1,G2,G3,G4,G5,CD1,CD2,CD3,CD4,CD5)

simulation = Simulation{T}("config_files/Si 1Hexagon 10.3-0.1.json")
simulation.detector.semiconductors[1].charge_density_model = ccdm

zs = 0:1e-6:0.002
charge_vec = similar(zs)
for iz in eachindex(zs)
    pt = CartesianPoint{T}(0,0,zs[iz])
    charge_vec[iz]=pt in simulation.detector.semiconductors[1] ?
    		      SolidStateDetectors.get_charge_density(simulation.detector.semiconductors[1], pt) : 0
end

zs2 = 0.0005:1e-6:0.0015
charge_vec2 = similar(zs2)
for iz in eachindex(zs2)
    pt = CartesianPoint{T}(0,0,zs2[iz])
    charge_vec2[iz]=pt in simulation.detector.semiconductors[1] ?
    		      SolidStateDetectors.get_charge_density(simulation.detector.semiconductors[1], pt) : 0
end

calculate_electric_potential!(simulation, max_refinements = 2,depletion_handling=true)

plot(zs, charge_vec, ylabel="impurity density", xlabel = "z")
png("1pixel_density_z")

plot(zs2, charge_vec2, ylabel="impurity density", xlabel = "z")
png("1pixel_density_z2")

plot(simulation.ρ)
png("1pixel_density")

plot(simulation.ρ, z=0.001995)
png("1pixel_density_top")

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
using DrWatson
using DelimitedFiles, Formatting

export read_state

#============== Parsing functions =============#

# rips blank space before new lines in Seba's data
function clean_txt_file(filename)
    open(filename, "r") do f
        s = read(f, String)
    end
    s = replace(s, " \r\n" => "\n")
    open(filename, "w") do f
        write(f, s)
    end
end

function she_p_to_name(gp::GenericParams)
    _p = @strdict N=gp.N dx=gp.Δx eps=gp.p[:ϵ]
    datadir("sims", "she", savename(_p, "values", digits=4))
end

# loads state for given params
function read_state(gp::GenericParams)
    filename = she_p_to_name(gp)
    clean_txt_file(filename)
    u = readdlm(filename, ' ', Float64, '\n')

    State(gp=gp, u=u)
end


text:
struct BlobBit
    data::Blob{UInt64}
    mask::UInt64
end

@inline Base.@propagate_inbounds function Base.getindex(blob::BlobBit)::Bool
    (blob.data[] & blob.mask) != 0
end

@inline Base.@propagate_inbounds function Base.setindex!(blob::BlobBit, v::Bool)::Bool
    c = blob.data[]
    c = ifelse(v, c | blob.mask, c & ~blob.mask)
    blob.data[] = c
    v
end

"A fixed-length bit vector whose data is stored in a Blob."
struct BlobBitVector <: AbstractArray{Bool, 1}
    data::Blob{UInt64}
    length::Int64
end

Base.@propagate_inbounds function get_address(blob::BlobBitVector, i::Int)::BlobBit
    @boundscheck begin
        (i < 1 || i > blob.length) && throw(BoundsError(blob, i))
    end
    i1, i2 = Base.get_chunks_id(i)
    BlobBit(blob.data + (i1-1)*self_size(UInt64), UInt64(1) << i2)
end

# blob interface

function Base.size(blob::Blob{BlobBitVector})
    (blob.length[],)
end

function Base.IndexStyle(_::Type{Blob{BlobBitVector}})
    Base.IndexLinear()
end

@inline Base.@propagate_inbounds function Base.getindex(blob::Blob{BlobBitVector}, i::Int)::BlobBit
    get_address(blob[], i)
end

function unsafe_resize!(blob::BlobBitVector, length::Int64)
    blob.length[] = length
end

# array interface

function Base.size(blob::BlobBitVector)
    (blob.length,)
end

function Base.IndexStyle(_::Type{BlobBitVector})
    Base.IndexLinear()
end

@inline Base.@propagate_inbounds function Base.getindex(blob::BlobBitVector, i::Int)::Bool
    get_address(blob, i)[]
end

@inline Base.@propagate_inbounds function Base.setindex!(blob::BlobBitVector, v::Bool, i::Int)::Bool
    get_address(blob, i)[] = v
end

@inline function Base.findprevnot(blob::BlobBitVector, start::Int)::Union{Nothing,Int}
    start > length(blob) && throw(BoundsError(blob, start))
    # TODO(tjgreen) placeholder slow implementation; should adapt optimized BitVector code
    @inbounds while start > 0 && blob[start]
        start -= 1
    end
    start > 0 ? start : nothing
end

@inline function Base.findnextnot(blob::BlobBitVector, start::Int)::Union{Nothing,Int}
    start > 0 || throw(BoundsError(blob, start))
    # TODO(tjgreen) placeholder slow implementation; should adapt optimized BitVector code
    @inbounds while start <= length(blob) && blob[start]
        start += 1
    end
    start <= length(blob) ? start : nothing
end

@inline function Base.findnext(blob::BlobBitVector, start::Int)::Union{Nothing,Int}
    start > 0 || throw(BoundsError(blob, start))
    # TODO(tjgreen) placeholder slow implementation; should adapt optimized BitVector code
    @inbounds while start <= length(blob) && !blob[start]
        start += 1
    end
    start <= length(blob) ? start : nothing
end


text:
using Base.Threads: sleep

@show ARGS

if isempty(ARGS)
    time = 5
else
    time = parse(Int, ARGS[1])
end

for t in 1:time
    println("Long: $t -> $time")
    println(stderr, "LongError: $t -> $time")
    sleep(rand([0.5, 0.8, 1]))
end
println("Finished after waiting $(time) seconds")

max_stars_repo_name: dmolina/DaemonMode

max_issues_repo_name: dmolina/DaemonMode

max_forks_repo_name: dmolina/DaemonMode


text:
using Documenter, HTsi

makedocs(;
    modules = [HTsi],
    format = Documenter.HTML(),
    pages = [
        "Index" => "index.md",
        "Home"  => "index.md",
        "Geomagnetic ndex"      => [
            "Ap Kp Index" => "ApKp.md",
            ],
        "Flux" => "flux.md",
        hide("Dst.md"),

    ],
    repo = "https://github.com/htyeim/HTsi.jl/blob/{commit}{path}#L{line}",
    sitename = "some space indices",
    authors = "htyeim <htyeim@gmail.com>",
    assets = String[],
)

deploydocs(;
    repo = "github.com/htyeim/HTsi.jl",
)


text:
using Metis
using Random
using Test
using SparseArrays
import LightGraphs

@testset "Metis.permutation" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    perm, iperm = Metis.permutation(S)
    @test isperm(perm)
    @test isperm(iperm)
    @test S == S[perm,perm][iperm,iperm]
end

@testset "Metis.partition" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    T = LightGraphs.smallgraph(:tutte)
    for G in (S, T), alg in (:RECURSIVE, :KWAY), nparts in (3, 4)
        partition = Metis.partition(G, nparts, alg = alg)
        @test extrema(partition) == (1, nparts)
        @test all(x -> findfirst(==(x), partition) !== nothing, 1:nparts)
    end
end

@testset "Metis.separator" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    T = LightGraphs.smallgraph(:tutte)
    for G in (S, T)
        parts = Metis.separator(G)
        @test extrema(parts) == (1, 3)
    end
end

@testset "Metis.MeshToDual" begin
    # Connectivity matrix of 2D mesh of 4 triangles.
    cells = [1 2 2 3;
             2 5 3 6;
             4 4 5 5]
    # Tell Metis to consider a cell connected if it shares two vertices.
    ncommon = 2

    # Get a Graph object. The number of graph verticies is the number of mesh
    # elements.
    G = Metis.graph(cells, ncommon)

    # Partition the graph.
    parts = Metis.partition(G, 2)
    # First two elements should be in the same partition.
    @test parts[1] == parts[2]
    # Last two elements should be in the same partition.
    @test parts[3] == parts[4]
end


text:
include("../src/Snappy.jl")
include("libsnappy.jl")

using Snappy
using BenchmarkTools

# this is rough, sorry!

testfiles = Dict([
    "alice29.txt" => "txt",
    "html" => "html",
    "fireworks.jpeg" => "jpeg",
    "paper-100k.pdf" => "pdf",
    "urls.10K" => "urls",
    "sample-tweet.json" => "json"
])

reports = []
for file in keys(testfiles)
    gc()
    println("\n\n")
    println(" ==== sample file: $file")
    input = read("$(Base.source_dir())/testdata/$(file)")

    # warm up
    ccall_uncompress(ccall_compress(input));
    uncompress(compress(input));

    m1_compress = @benchmark compress($input)
    println(" ======== compress (Julia):")
    display(m1_compress)
    println("")
    m0_compress = @benchmark ccall_compress($input)
    println(" ======== compress (ccall):")
    display(m0_compress)
    println("")
    println("---------------------------")
    println(" ============ compress comparison (median):")
    j = judge(median(m1_compress), median(m0_compress))
    display(j)
    println("\n")
    julia_bps = length(input) / (median(m1_compress).time / 10^9)
    ccall_bps = length(input) / (median(m0_compress).time / 10^9)
    if julia_bps >= 1e9
        julia_pretty = "$(round(julia_bps / 2^30, 3)) GB/s"
    elseif julia_bps >= 1e8
        julia_pretty = "$(Int(floor(julia_bps / 2^20))) MB/s"
    else
        julia_pretty = "$(Int(floor(julia_bps / 2^10))) KB/s"
    end
    if ccall_bps >= 1e9
        ccall_pretty = "$(round(ccall_bps / 2^30, 3)) GB/s"
    elseif ccall_bps >= 1e8
        ccall_pretty = "$(Int(floor(ccall_bps / 2^20))) MB/s"
    else
        ccall_pretty = "$(Int(floor(ccall_bps / 2^10))) KB/s"
    end
    j_pretty = j.ratio.time > 1 ? "+$(round((j.ratio.time-1)*100,2))%" : "-$(round((1-j.ratio.time)*100,2))%"

    outline_a = "|$(testfiles[file])|$(length(input))|$julia_pretty|$ccall_pretty|$j_pretty"
    gc()
    output = compress(input);
    m1_uncompress = @benchmark uncompress($output)
    println(" ======== uncompress (Julia):")
    display(m1_uncompress)
    println("")
    output = ccall_compress(input);
    m0_uncompress = @benchmark ccall_uncompress($output)
    println(" ======== uncompress (ccall):")
    display(m0_uncompress)
    println("")
    println("---------------------------")
    println(" ============ uncompress comparison (median):")
    j = judge(median(m1_uncompress), median(m0_uncompress))
    display(j)
    println("\n")
    julia_bps = length(output) / (median(m1_uncompress).time / 10^9)
    ccall_bps = length(output) / (median(m0_uncompress).time / 10^9)
    if julia_bps >= 1e9
        julia_pretty = "$(round(julia_bps / 2^30, 3)) GB/s"
    elseif julia_bps >= 1e8
        julia_pretty = "$(Int(floor(julia_bps / 2^20))) MB/s"
    else
        julia_pretty = "$(Int(floor(julia_bps / 2^10))) KB/s"
    end
    if ccall_bps >= 1e9
        ccall_pretty = "$(round(ccall_bps / 2^30, 3)) GB/s"
    elseif ccall_bps >= 1e8
        ccall_pretty = "$(Int(floor(ccall_bps / 2^20))) MB/s"
    else
        ccall_pretty = "$(Int(floor(ccall_bps / 2^10))) KB/s"
    end
    j_pretty = j.ratio.time > 1 ? "+$(round((j.ratio.time-1)*100,2))%" : "-$(round((1-j.ratio.time)*100,2))%"

    outline_b = "|$julia_pretty|$ccall_pretty|$j_pretty|"
    println(outline_a * outline_b)
    push!(reports, outline_a*outline_b)
    println("==================================================")
end
println("final report:")
println(join(reports, '\n'))


text:
push!(LOAD_PATH, abspath("src"))
using Yao
using QMPS
using DelimitedFiles, JLD2, FileIO, Pkg

# CUDA switch
const USE_CUDA = haskey(Pkg.installed(), "CuYao")
USE_CUDA && println("Hint: Using CUDA since `CuYao` is detected. Edit file `applications.jl` to modify CUDA settings, like switching computing devices.")
USE_CUDA && include("CuChem.jl")
USE_CUDA && device!(CuDevice(0))

"""Heisenberg model without frustration, with open boundary condition."""
simple_model_heis(size...) = Heisenberg(size...; periodic=false)

"""Heisenberg model with frustration strenth J2 = 0.5, with open boundary condition."""
simple_model_j1j2(size...) = J1J2(size...; periodic=false, J2=0.5)

"""
    simple_ansatz(nbit::Int, symmetry::Symbol, depth::Int; load_params::Bool=false)

Load a predefined MPS inspired ansatz with 4 virtual qubits, batch size 4096.
If `load_params` is `true`, load parameters in training step 500.
"""
function simple_ansatz(nbit::Int, symmetry::Symbol, depth::Int; load_params::Bool=false)
    V = 4     # number of virtual qubits
    batch_size = 4096
    load_step = 500

    ansatz = model(Val(symmetry); nbit=nbit, nlayer=depth, V=V, B=batch_size, pairs=pair_ring(V+1))
    USE_CUDA && (ansatz = ansatz |> cu)

    if load_params
        filename = "data/chem_$(symmetry)_d$(depth)_N$(nbit)_V$(V)_S$(load_step).jld2"
        params = load_training(filename)[3]
        println("loading file $filename")
        println("Number of parameters is ", params |> length)
        dispatch!(ansatz.circuit, params)
    end

    ansatz
end

"""
    save_training(filename, qopt::Adam, loss, params, fidelity)

Save training status.
"""
function save_training(filename, qopt::Adam, loss, params, fidelity)
    save(filename, "qopt", qopt, "loss", loss, "params", params, "fidelity", fidelity)
end

"""
    load_training(filename) -> Tuple

Load training status (qopt, loss, params, fidelity).
"""
function load_training(filename)
    res = load(filename)
    res["qopt"], res["loss"], res["params"], res["fidelity"]
end

"""
    run_train(ansatz, model; SAVE_ID, niter=500, start_point=0, save_step=0)

Run a variational quantum eigensolver to solve a `model` Hamiltonian with MPS inspired `ansatz`

* SAVE_ID, token used for saving data,
* niter, number of training steps,
* start_point, load a save point, default is 0 (random parameters).
* save_step, save the training every `save_step` steps.
"""
function run_train(ansatz, model; SAVE_ID, niter::Int=500, start_point::Int=0, save_step::Int=10)
    nbit = nbit_simulated(ansatz)
    V = ansatz.nbit_virtual
    filename(k::Int) = "data/chem_$(SAVE_ID)_N$(nbit)_V$(V)_S$(k).jld2"
    if start_point==0
        optimizer = Adam(lr=0.1)
        history = Float64[]
        fidelities = Float64[]
    else
        optimizer, history, _params, fidelities = load_training(filename(start_point))
        dispatch!(ansatz.circuit, _params)
    end
    qo = QMPSOptimizer(ansatz, model, optimizer)
    gs = ground_state(model)

    for (k_, p) in enumerate(qo)
        k = k_ + start_point
        curr_loss = energy(ansatz, model)/nbit
        fid = fidelity_exact(ansatz, gs)[]
        push!(history, curr_loss)
        push!(fidelities, fid)
        println("step = $k, energy/site = $curr_loss, fidelity = $(fid)")
        flush(stdout)
        if k%save_step == 0 || k == niter
            save_training(filename(k), qo.optimizer, history, parameters(ansatz.circuit), fidelities)
            k >= niter && break
        end
    end
end

"""
    correlation_matrix(ansatz, op; SAVE_ID)

Calculate and save correlation matrix <Z_i Z_j>.
"""
function correlation_matrix(ansatz; SAVE_ID)
    nbit = nbit_simulated(ansatz)
    V = ansatz.nbit_virtual

    om = zeros(Float64, nbit, nbit)
    for i =1:nbit, j=1:nbit
        if i!=j
            om[i,j] = measure_corr(ansatz, i=>Z, j=>Z) |> real
            println("<σz($i)σz($j)> = $(om[i,j])")
        end
    end

    for (token, var) in [
                         ("om", om),
                        ]
        writedlm("data/chem_$(SAVE_ID)_$(token)$(tag)_N$(nbit)_V$V.dat", var)
    end
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
function validate_available_devices(
    device::DeviceModel{T, U},
    system,
) where {T <: PSY.Device, U <: AbstractDeviceFormulation}
    devices = get_available_components(T, system)
    if isempty(devices)
        @info "The data doesn't include devices of type $(T), consider changing the device models" _group =
            LOG_GROUP_MODELS_VALIDATION
        return false
    end
    PSY.check_components(system, devices)
    return true
end

function validate_service!(
    model::ServiceModel{S, <:AbstractServiceFormulation},
    incompatible_device_types::Set{<:DataType},
    sys::PSY.System,
) where {S <: PSY.Service}
    service = PSY.get_component(S, sys, get_service_name(model))
    if service === nothing
        @info "The data doesn't include services of type $(S) and name $(get_service_name(model)), consider changing the service models" _group =
            LOG_GROUP_MODELS_VALIDATION
        return false
    end

    PSY.check_components(sys, service)
    services_mapping = PSY.get_contributing_device_mapping(sys)

    contributing_devices_ =
        services_mapping[(type=S, name=PSY.get_name(service))].contributing_devices
    contributing_devices = [
        d for d in contributing_devices_ if
        typeof(d) ∉ incompatible_device_types && PSY.get_available(d)
    ]
    if isempty(contributing_devices)
        @warn "The contributing devices for service $(PSY.get_name(service)) is empty, consider removing the service from the system" _group =
            LOG_GROUP_MODELS_VALIDATION
        return false
    end

    PSY.check_components(sys, contributing_devices)
    return true
end

function validate_services!(
    model::ServiceModel{S, <:AbstractServiceFormulation},
    ::Vector{<:DataType},
    sys::PSY.System,
) where {S <: PSY.StaticReserveGroup}
    service = PSY.get_component(S, sys, get_service_name(model))
    if service === nothing
        @info "The data doesn't include services of type $(S) and name $(get_service_name(model)), consider changing the service models" _group =
            LOG_GROUP_MODELS_VALIDATION
        return false
    end

    PSY.check_components(sys, service)
    contributing_services = PSY.get_contributing_services(s)
    if isempty(contributing_services)
        @warn "The contributing services for group service $(PSY.get_name(service)) is empty, consider removing the group service from the system" _group =
            LOG_GROUP_MODELS_VALIDATION
        return false
    end

    PSY.check_components(sys, contributing_devices)
    return true
end


text:
using SimpleSolvers
using Test

n = 2
x = rand(n)
h = zeros(n,n)
T = eltype(x)

function F(x::Vector)
    1 + sum(x.^2)
end

function H!(h::Matrix, x::Vector)
    h .= 0
    for i in eachindex(x)
        h[i,i] = 2
    end
end

H!(h,x)

HPAD = Hessian{T}(F, n; mode = :autodiff)
HPUS = Hessian{T}(H!, n; mode = :user)

@test typeof(HPAD) <: HessianAD
@test typeof(HPUS) <: HessianUser

function test_hessian(h1, h2, atol)
    for i in eachindex(h1,h2)
        @test h1[i] ≈ h2[i] atol=atol
    end
end

had = zero(h)
hus = zero(h)

compute_hessian!(had, x, HPAD)
compute_hessian!(hus, x, HPUS)

test_hessian(had, h, eps())
test_hessian(hus, h, 0)

had1 = zero(h)
hus1 = zero(h)

compute_hessian!(had1, x, F; mode = :autodiff)
compute_hessian!(hus1, x, H!; mode = :user)

test_hessian(had, had1, 0)
test_hessian(hus, hus1, 0)

had2 = zero(h)

compute_hessian_ad!(had2, x, F)

test_hessian(had, had2, 0)


text:
struct Wine <: Tabular end
url(::Wine) = "https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data"
checksum(::Wine) = "6be6b1203f3d51df0b553a70e57b8a723cd405683958204f96d23d7cd6aea659"
prep(::Wine) = path -> preprocess(path, Wine())
target(::Wine) = 1
size(::Wine) = (178, 0, 0)
problem(::Wine) = Classification

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
export SerreGreenNaghdi

"""
    SerreGreenNaghdi(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the Serre-Green-Naghdi model

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`.

## Keywords
- `iterative`: solve the elliptic problem through GMRES if `true`, LU decomposition if `false` (default is `true`);
- `precond`: use a (left) preconditioner for GMRES if `true` (default), choose `precond` as the preconditioner if provided;
- `gtol`: relative tolerance of the GMRES algorithm (default is `1e-14`);
- `restart`: the corresponding option of the GMRES algorithm (default is `100`);
- `maxiter`: the corresponding option of GMRES (default is `nothing`);
- `ktol`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!` and in particular
1. a function `SerreGreenNaghdi.f!` to be called in the time-integration solver;
2. a function `SerreGreenNaghdi.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `SerreGreenNaghdi.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where

    - `η` is the surface deformation;
    - `v` is the derivative of the trace of the velocity potential;
4. additionally, a handy function `SerreGreenNaghdi.mapfrofull` which from data matrix returns the Tuple of real vectors `(η,v,u)`, where

	- `u` corresponds to the layer-averaged velocity.

"""
mutable struct SerreGreenNaghdi <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	mapfrofull	:: Function
	param	:: NamedTuple
	kwargs  :: NamedTuple

    function SerreGreenNaghdi(param::NamedTuple;dealias=0,ktol=0,iterate=true,gtol=1e-14,precond=true,restart=nothing,maxiter=nothing,verbose=true)
		m=WhithamGreenNaghdi(param;SGN=true,dealias=dealias,ktol=ktol,iterate=iterate,gtol=gtol,precond=precond,restart=restart,maxiter=maxiter,verbose=verbose)
        new(m.label, m.f!, m.mapto, m.mapfro, m.mapfrofull, m.param, m.kwargs)
    end
end


text:
#Importing parameters
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport matplotlib.pyplot as plt
@pyimport numpy as np
@pyimport os
@pyimport pickle
@pyimport pandas as pd
@pyimport collections
@pyimport parameters_julia

out_dir = parameters_julia.out_dir
files_ID = parameters_julia.files_ID
month_w = parameters_julia.month_w
year = parameters_julia.year
instances_ = parameters_julia.instances_ID

instance_ = instances_[1];

using JSON

instance1 = instance_

open(out_dir * "instance_comm.txt", "w") do f
    write(f, instance1)
end

week_day_Apr_list = parameters_julia.week_day_list
#day = week_day_Apr_list[1]

using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport utils_julia
@pyimport GLS_julia

using JuMP, Ipopt
cnt = 0
poaDict = Dict{}()
social = Dict{}()
user = Dict{}()
flow_user, link_vector = GLS_julia.GLS_juliaf(instance_);

#load OD pair-route incidence
odPairRoute = readstring(out_dir * "od_pair_route_incidence_" * instance1 *  files_ID * ".json");
odPairRoute = JSON.parse(odPairRoute);

for route in keys(odPairRoute)
    if odPairRoute[route]>0
        odPairRoute[route] = 1
    end
end

#load link-route incidence
linkRoute = readstring(out_dir * "link_route_incidence_" * instance1 *  files_ID * ".json");
linkRoute = JSON.parse(linkRoute);

#load OD pair labels
odPairLabel = readstring(out_dir *  "od_pair_label_dict_refined.json");
odPairLabel = JSON.parse(odPairLabel);

odPairLabel_ = readstring(out_dir  *"od_pair_label_dict__refined.json");
odPairLabel_ = JSON.parse(odPairLabel_);

#load link labels
linkLabel = readstring(out_dir * "link_label_dict.json");
linkLabel = JSON.parse(linkLabel);

linkLabel_ = readstring(out_dir * "link_label_dict_.json");
linkLabel_ = JSON.parse(linkLabel_);

#load node-link incidence
nodeLink = readstring(out_dir * "node_link_incidence.json");
nodeLink = JSON.parse(nodeLink);

include("Julia_files/load_network_uni_class.jl");

for day in week_day_Apr_list

	#load OD pair-route incidence
	#odPairRoute = readstring(out_dir * "demandsDict/jacobi" * string(day) * "_" * month_w * "_" *instance1 * ".json");
	#odPairRoute = JSON.parse(odPairRoute);
	#odPairRoute = odPairRoute[:,:]
	#println(odPairRoute)
	#println(length(odPairRoute))
	#println(length(odPairRoute[1]))
	#odPairRoute = reshape(odPairRoute, length(odPairRoute), length(odPairRoute[1]))

	ta_data = load_ta_network_(out_dir, files_ID, month_w, string(day) , instance1);

	capacity = ta_data.capacity;
	free_flow_time = ta_data.free_flow_time;
	#=
	############
	#Read in the demand file
	file = open(out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_trips_" * month_w * "_" * string(day) * "_" * instance1 * ".txt")
	demands = Dict()
	s = 0
	for line in eachline(file)
	    if contains(line, "Origin")
	        s = Int(parse(Float64,(split(line)[2])))
	    elseif contains(line, ";")
	        pairs = split(line, ";")
	        for pair in pairs[1:end-1]
	            if !contains(pair, "\n")
	                pair_vals = split(pair, ":")
	                t, demand = Int(parse(Float64,pair_vals[1])), parse(Float64,pair_vals[2])
	                demands[(s,t)] = demand
	            end
	        end
	    end
	end
	close(file);
	print(demands)
	=#
	demands = readstring(out_dir * "demandsDict/demandsDictFixed" * string(day) * "_" * month_w * "_"  * instance1 * ".json")
	demands = JSON.parse(demands)

	#println(demands)

	demandsVec = zeros(length(odPairLabel_))

	println(demandsVec)

	for i = 1:length(demandsVec)
		#println( "(" * string(odPairLabel_["$i"][1]) * ", " * string(odPairLabel_["$i"][2]) * ")" )
	    demandsVec[i] = demands["(" * string(odPairLabel_["$i"][1]) * ", " * string(odPairLabel_["$i"][2]) * ")"]
	end
	println(demandsVec)
	for key=keys(odPairRoute)
	    if contains(key, "56-")
	        println(key)
	    end
	end

	coeffs_dict_Apr_PM_ = readstring(out_dir * "coeffs_dict_" * month_w *  "_" * instance1 * ".json")
	coeffs_dict_Apr_PM_ = JSON.parse(coeffs_dict_Apr_PM_)

	best_key = readstring(out_dir * "cross_validation_best_key/cross_validation_best_key_" * month_w * "_" * string(day) * "_" * instance1 * ".json")
	best_key = JSON.parse(best_key)
	#best_key = "(7, 0.5, 1000.0, 1)"
	fcoeffs = coeffs_dict_Apr_PM_[best_key]
	polyDeg = length(fcoeffs)

	# m = Model(solver=GurobiSolver(OutputFlag=false))
	m = Model(solver=IpoptSolver());

	numLinks = length(linkLabel_)
	numRoute = length(linkRoute)
	#numRoute = 400

	numOD = length(demandsVec)

	@variable(m, linkFlow[1:numLinks])
	@variable(m, pathFlow[1:numRoute])

	pathFlowSum = Dict()

	for i=1:numOD
	    pathFlowSum[i] = 0
	    for j=1:numRoute
	        if "$(i)-$(j)" in keys(odPairRoute)
	            pathFlowSum[i] += pathFlow[j]
	        end
	    end
	    @constraint(m, pathFlowSum[i] == demandsVec[i])
	end

	pathFlowLinkSum = Dict()

	for a=1:numLinks
	    pathFlowLinkSum[a] = 0
	    for j=1:numRoute
	        if "$(a)-$(j)" in keys(linkRoute)
	            pathFlowLinkSum[a] += pathFlow[j];
	        end
	    end
	    @constraint(m, pathFlowLinkSum[a] == linkFlow[a]);
	end

	for j=1:numRoute
	    @constraint(m, pathFlow[j] >= 0);
	end

	#@expression(m, f, sum{free_flow_time[a]*linkFlow[a] + .03*free_flow_time[a]*((linkFlow[a])^5)/((capacity[a])^4), a = 1:numLinks} )

	@NLexpression(m, f, sum{ free_flow_time[a] * fcoeffs[i]  *linkFlow[a]^i / capacity[a]^(i-1) , i = 1:polyDeg , a = 1:numLinks }) ;
#=
	@NLexpression(m, f, sum{free_flow_time[a] * fcoeffs[1] * linkFlow[a] +
	        free_flow_time[a] * fcoeffs[2] * linkFlow[a]^2 / capacity[a] +
	        free_flow_time[a] * fcoeffs[3] * linkFlow[a]^3 / capacity[a]^2 +
	        free_flow_time[a] * fcoeffs[4] * linkFlow[a]^4 / capacity[a]^3 +
	        free_flow_time[a] * fcoeffs[5] * linkFlow[a]^5 / capacity[a]^4 +
	        free_flow_time[a] * fcoeffs[6] * linkFlow[a]^6 / capacity[a]^5 +
			free_flow_time[a] * fcoeffs[7] * linkFlow[a]^7 / capacity[a]^6 , a = 1:numLinks})
=#
	@NLobjective(m, Min, f);
	#print(m)

	solve(m);

	flows = Dict();

	for i = 1:length(ta_data.start_node)
	    key = (ta_data.start_node[i], ta_data.end_node[i]);
	    flows[key] = getvalue(linkFlow)[i];
	end

	function socialObj(linkFlowVec)
	    objVal =  sum(sum(free_flow_time[a] * fcoeffs[i] * linkFlowVec[a]^i / capacity[a]^(i-1) for i=1:polyDeg)
	        for a = 1:numLinks)
	    return objVal
	end

	PoA_dict = Dict();
	tapSocialFlowDicDict = Dict();
	tapSocialFlowVecDict = Dict();

	#cnt = 0
	#for i = 1:length(week_day_Apr_list)
	#poaDictAprPM =  socialObj(flow_user[:,cnt+1])/getobjectivevalue(m)
	#println(socialObj(flow_user[:, cnt]))
	cnt += 1

	user_flow_tap_MSA = readstring(out_dir * "demandsDict/tapFlowDicDict" * string(day) * "_" * month_w * "_"  * instance1 * ".json")
	user_flow_tap_MSA = JSON.parse(user_flow_tap_MSA)

	k = maximum(map(k->parse(Int,k),keys(user_flow_tap_MSA)))
	user_flow_tap_MSA = user_flow_tap_MSA[string(k)]

	user_flow = Dict();
	for a = 1:length(linkLabel)
		#println( "(" * string(odPairLabel_["$i"][1]) * ", " * string(odPairLabel_["$i"][2]) * ")" )
	    user_flow[a] = user_flow_tap_MSA[ "(" * string(linkLabel[string(a-1)][1]) * ", " * string(linkLabel[string(a-1)][3]) * ")"]
	    #demandsVec[i] = demands["(" * string(odPairLabel_["$i"][1]) * ", " * string(odPairLabel_["$i"][2]) * ")"]
	end

	social[day] = getobjectivevalue(m)
	user[day] = socialObj(user_flow)
	#user[day] = socialObj(flow_user[:, cnt])
	poaDict[day] = (user[day]/social[day])

end

outfile =  open(out_dir * "cong_" * month_w * "_" * instance1 * ".json", "w")
JSON.print(outfile, user)
close(outfile)

outfile =  open(out_dir * "PoA_dict_" * month_w * "_" * instance1 * ".json", "w")
JSON.print(outfile, poaDict)
close(outfile)

user
social

poaDict

#=

=#

max_stars_repo_path: 05. TAP POA.jl

max_issues_repo_path: 05. TAP POA.jl

max_forks_repo_path: 05. TAP POA.jl


text:
"""
    STD_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)

Performs standard probing using linear underestimator forms via
McCormick relaxations, sparse calculations, and contracts `X::Vector{Interval{Float64}}` in place.
The `opt` object has the fields `.numVar`, `.numConstr`, `.f`, `.g`, `.gL`, `.gU`,
`.gL_Loc`, `.gU_Loc` and `solver.LP_solver`. `f(x)` is the objective, `g(x)` is
the constraint function, `numVar` is the number of decision variables, `numConstr`
is the number of constraints, `gL` is the lower bound, `gU` is the upper bound,
`gL_Loc` is an index to indicate the lower bound is finite, and `gU_Loc` is an index
that indicates the upper bound is finite. The upper bound is `UBD::Float64`.
"""
function STD_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)

    mu_temp::Int64 = copy(EAGOSmoothMcCormickGrad.MC_param.mu)
    set_diff_relax(0)

    # constructs LP relaxation
    Xlo::Vector{Float64} = [X[i].lo for i=1:opt[1].numVar]
    Xhi::Vector{Float64} = [X[i].hi for i=1:opt[1].numVar]
    x0::Vector{Float64} = (Xlo + Xhi)/2.0
    x_SMC::Vector{SMCg{opt[1].numVar,Float64}} = [SMCg{opt[1].numVar,Float64}(x0[i],
                                                                              x0[i],
                                                                              seed_g(Float64,i,opt[1].numVar),
                                                                              seed_g(Float64,i,opt[1].numVar),
                                                                              X[i],
                                                                              false,
                                                                              X,
                                                                              x0) for i=1:opt[1].numVar]
    # probes upper bound
    f_mc::SMCg{opt[1].numVar,Float64} = opt[1].f(x_SMC)
    f_cv::Float64 = f_mc.cv
    if opt[1].numConstr>0
        c::Vector{SMCg{opt[1].numVar,Float64}} = opt[1].g(x_SMC)
        dcdx::SparseMatrixCSC{Float64,Int64} = spzeros(length(opt[1].gL_loc)+length(opt[1].gU_loc)+1,opt[1].numVar)
    else
        dcdx = spzeros(1,opt[1].numVar)
    end

    dcdx[1,:] = f_mc.cv_grad
    if opt[1].numConstr>0
        cx_ind1::Int64 = 2
        for i in opt[1].gL_loc
            for j=1:opt[1].numVar
                if (c[i].cv_grad[j] != 0.0)
                    dcdx[cx_ind1,j] = c[i].cv_grad[j]
                end
            end
            cx_ind1 += 1
        end
        for i in opt[1].gU_loc
            for j=1:opt[1].numVar
                if (c[i].cc_grad[j] != 0.0)
                    dcdx[cx_ind1,j] = -c[i].cc_grad[j]
                end
            end
            cx_ind1 += 1
        end
    end

    # forms rhs of linear relaxations
    if opt[1].numConstr>0
        rhs::Vector{Float64} = zeros(Float64,length(opt[1].gL_loc)+length(opt[1].gU_loc)+1)
    else
        rhs = zeros(Float64,1)
    end
    rhs[1] = UBD + f_mc.cv - sum(x0.*f_mc.cv_grad)

    if opt[1].numConstr>0
        cx_ind2::Int64 = 2
        for i in opt[1].gU_loc
            rhs[cx_ind2] = sum(x0[:].*c[i].cv_grad[:])+opt[1].gU[i]-c[i].cv
            cx_ind2 += 1
        end
        for i in opt[1].gL_loc
            rhs[cx_ind2] = sum(-x0[:].*c[i].cc_grad[:])-opt[1].gL[i]+c[i].cc
            cx_ind2 += 1
        end
    end

    if (opt[1].numConstr>0)
        temp_model = buildlp([f_mc.cv_grad[i] for i=1:opt[1].numVar], dcdx, '<', rhs, Xlo, Xhi, opt[1].solver.LP_solver)
    else
        temp_model = buildlp([f_mc.cv_grad[i] for i=1:opt[1].numVar], zeros(opt[1].numVar,opt[1].numVar), '<', zeros(opt[1].numVar), Xlo, Xhi, opt[1].solver.LP_solver)
    end

    for i=1:opt[1].numVar

        # probes upper bound
        MathProgBase.setvarLB!(temp_model,[i == j ? Xhi[i] : Xlo[i] for j=1:opt[1].numVar])
        MathProgBase.setvarUB!(temp_model,Xhi)
        result = solvelp(temp_model)
        if (result.status == :Optimal)
            val::Float64 = result.objval + f_cv - sum([x0[i]*f_mc.cv_grad[i] for i=1:opt[1].numVar])
            pnt::Vector{Float64} = result.sol
            mult::Vector{Float64} = result.attrs[:redcost]
            mult_lo::Vector{Float64} = [tol_eq(X[i].lo,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            mult_hi::Vector{Float64} = [tol_eq(X[i].hi,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            Variable_DR!(X,mult_lo,mult_hi,val,UBD)
        end
        # probes lower bound
        MathProgBase.setvarLB!(temp_model,Xlo)
        MathProgBase.setvarUB!(temp_model,[i == j ? Xhi[i] : Xlo[i] for j=1:opt[1].numVar])
        result = solvelp(temp_model)
        if (result.status == :Optimal)
            val = result.objval + f_cv - sum([x0[i]*f_mc.cv_grad[i] for i=1:opt[1].numVar])
            pnt = result.sol
            mult = result.attrs[:redcost]
            mult_lo = [tol_eq(X[i].lo,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            mult_hi = [tol_eq(X[i].hi,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            Variable_DR!(X,mult_lo,mult_hi,val,UBD)
        end
    end

    # reset McCormick relaxation to original value
    set_diff_relax(mu_temp)

    # stores outputs
    X[:] = X
    return true
end

"""
    Imp_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)

Performs standard probing using linear underestimator forms via
McCormick relaxations, sparse calculations, and contracts `X::Vector{Interval{Float64}}`
in place only tightening the following intervals in the vector `X[(nx+1):(nx+np)]`.
The `opt` object has the fields `.numVar`, `.numConstr`, `.f`, `.g`, `.gL`, `.gU`,
`.gL_Loc`, `.gU_Loc` and `solver.LP_solver`. `f(x)` is the objective, `g(x)` is
the constraint function, `numVar` is the number of decision variables, `numConstr`
is the number of constraints, `gL` is the lower bound, `gU` is the upper bound,
`gL_Loc` is an index to indicate the lower bound is finite, and `gU_Loc` is an index
that indicates the upper bound is finite. The upper bound is `UBD::Float64`.
"""
function Imp_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)

    # constructs LP relaxation
    Xlo::Vector{Float64} = [X[i].lo for i=1:opt[1].numVar]
    Xhi::Vector{Float64} = [X[i].hi for i=1:opt[1].numVar]
    x0::Vector{Float64} = (Xlo + Xhi)/2.0
    x_SMC::Vector{SMCg{opt[1].numVar,Float64}} = [SMCg{opt[1].numVar,Float64}(x0[i],
                                                                              x0[i],
                                                                              seed_g(Float64,i,opt[1].numVar),
                                                                              seed_g(Float64,i,opt[1].numVar),
                                                                              X[i],
                                                                              false,
                                                                              X,
                                                                              x0) for i=1:opt[1].numVar]
    # probes upper bound
    f_mc::SMCg{opt[1].numVar,Float64} = opt[1].f(x_SMC)
    f_cv::Float64 = f_mc.cv
    c::Vector{SMCg{opt[1].numVar,Float64}} = opt[1].g(x_SMC)
    # forms coefficients of linear relaxations
    dcdx::SparseMatrixCSC{Float64,Int64} = spzeros(length(opt[1].gL_loc)+length(opt[1].gU_loc)+1,opt[1].numVar)
    dcdx[1,:] = f_mc.cv_grad
    cx_ind1::Int64 = 2
    for i in opt[1].gL_loc
        for j=1:opt[1].numVar
            if (c[i].cv_grad[j] != 0.0)
                dcdx[cx_ind1,j] = c[i].cv_grad[j]
            end
        end
        cx_ind1 += 1
    end
    for i in opt[1].gU_loc
        for j=1:opt[1].numVar
            if (c[i].cc_grad[j] != 0.0)
                dcdx[cx_ind1,j] = -c[i].cc_grad[j]
            end
        end
        cx_ind1 += 1
    end

    # forms rhs of linear relaxations
    rhs::Vector{Float64} = zeros(Float64,length(opt[1].gL_loc)+length(opt[1].gU_loc)+1)
    rhs[1] = UBD + f_mc.cv - sum(x0.*f_mc.cv_grad)
    cx_ind2::Int64 = 2
    for i in opt[1].gU_loc
        rhs[cx_ind2] = sum(x0[:].*c[i].cv_grad[:])+opt[1].gU[i]-c[i].cv
        cx_ind2 += 1
    end
    for i in opt[1].gL_loc
        rhs[cx_ind2] = sum(-x0[:].*c[i].cc_grad[:])-opt[1].gL[i]+c[i].cc
        cx_ind2 += 1
    end

    if (opt[1].numConstr>0)
        temp_model = buildlp([f_mc.cv_grad[i] for i=1:opt[1].numVar], dcdx, '<', rhs, Xlo, Xhi, opt[1].solver.LP_solver)
    else
        temp_model = buildlp([f_mc.cv_grad[i] for i=1:opt[1].numVar], zeros(opt[1].numVar,opt[1].numVar), '<', zeros(opt[1].numVar), Xlo, Xhi, opt[1].solver.LP_solver)
    end

    for i=1:opt[1].numVar

        # probes upper bound
        MathProgBase.setvarLB!(temp_model,[i == j ? Xhi[i] : Xlo[i] for j=1:opt[1].numVar])
        MathProgBase.setvarUB!(temp_model,Xhi)
        result = solvelp(temp_model)
        if (result.status == :Optimal)
            val::Float64 = result.objval + f_cv - sum([x0[i]*f_mc.cv_grad[i] for i=1:opt[1].numVar])
            pnt::Vector{Float64} = result.sol
            mult::Vector{Float64} = result.attrs[:redcost]
            mult_lo::Vector{Float64} = [tol_eq(X[i].lo,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            mult_hi::Vector{Float64} = [tol_eq(X[i].hi,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            Variable_DR!(X,mult_lo,mult_hi,val,UBD)
        end
        # probes lower bound
        MathProgBase.setvarLB!(temp_model,Xlo)
        MathProgBase.setvarUB!(temp_model,[i == j ? Xhi[i] : Xlo[i] for j=1:opt[1].numVar])
        result = solvelp(temp_model)
        if (result.status == :Optimal)
            val = result.objval + f_cv - sum([x0[i]*f_mc.cv_grad[i] for i=1:opt[1].numVar])
            pnt = result.sol
            mult = result.attrs[:redcost]
            mult_lo = [tol_eq(X[i].lo,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            mult_hi = [tol_eq(X[i].hi,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
            Variable_DR!(X,mult_lo,mult_hi,val,UBD)
        end
    end
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using DelimitedFiles
using JLD2
Random.seed!(233)
ADCME.options.mpi.printlevel = 0

mutable struct MPIConfig
    n::Int64
    N::Int64
    h::Float64
    function MPIConfig(n::Int64)
        r = mpi_size()
        m = Int64(round(sqrt(r)))
        if m^2 != r
            error("$r must be a squared number")
        end
        h = 1/(n*m+1)
        new(n, m, h)
    end
end

function getIJ(mc::MPIConfig)
    r = mpi_rank()
    N = mc.N
    I = r÷N + 1
    J = r%N + 1
    return I, J
end

function get_xy(mc::MPIConfig)
    n, N, h = mc.n, mc.N, mc.h
    x = Float64[]
    y = Float64[]
    I, J = getIJ(mc)
    x0 = (I-1)*n*h
    y0 = (J-1)*n *h
    for i = 1:n
        for j = 1:n
            push!(x, i*h+x0)
            push!(y, j*h+y0)
        end
    end
    X, Y = Array(reshape(x, (mc.n, mc.n))'), Array(reshape(y, (mc.n, mc.n))')
end

"""
Given an extended `κ_ext` (size = `(n+2)×(n+2)`), solves the forward computation problem
"""
function poisson_solver(κ_local, f_local, mc::MPIConfig)
    @assert size(κ_local) == (mc.n, mc.n)
    @assert size(f_local) == (mc.n, mc.n)
    global κ_ext = mpi_halo_exchange(κ_local, mc.N, mc.N; fill_value = 1.0)
    h = mc.h
    f_local = constant(f_local)
    rhs = 2*h^2*reshape(f_local, (-1,))

    poisson_linear_solver(κ_ext, rhs, mc)
end

function poisson_linear_solver(κ_ext, rhs, mc)
    function forward(κ_ext)
        A = get_poisson_matrix(κ_ext, mc)
        A\rhs
    end

    function backward(du, u, κ_ext)
        A = get_poisson_matrix(κ_ext, mc)
        B = A'
        x = -(B\du)
        x = reshape(x, (mc.n, mc.n))
        u = reshape(u, (mc.n, mc.n))
        uext = mpi_halo_exchange(u, mc.N, mc.N; fill_value = 0.0)
        out = get_poisson_grad(x, uext, mc.N)
        set_shape(out, (mc.n+2, mc.n+2))
    end

    register(forward, backward)(κ_ext)
end

function get_poisson_grad(x,uext,N)
    get_poisson_grad_ = load_op_and_grad("$(@__DIR__)/build/libGetPoissonMatrix","get_poisson_grad")
    x,uext,cn = convert_to_tensor(Any[x,uext,N], [Float64,Float64,Int64])
    get_poisson_grad_(x,uext,cn)
end

function dofmap(mc::MPIConfig)
    n, N = mc.n, mc.N
    global_to_local = zeros(Int64, (n*N)^2)
    local_to_global = zeros(Int64, (n*N)^2)
    for I = 1:N
        for J = 1:N
            for i = 1:n
                for j = 1:n
                    ii = (I-1)*n + i
                    jj = (J-1)*n + j
                    global_to_local[ (ii-1)*n*N + jj ] = ((I-1)*N+J-1)*n^2 + (i-1)*n+j
                    local_to_global[ ((I-1)*N+J-1)*n^2 + (i-1)*n+j ] = (ii-1)*n*N + jj
                end
            end
        end
    end
    return global_to_local, local_to_global
end

function get_colext(mc::MPIConfig)
    I, J = getIJ(mc)
    n, N = mc.n, mc.N
    NDOF = n*N
    DMAP, _ = dofmap(mc)

    colext = zeros(Int64, n+2, n+2)
    for i = 1:n+2
        for j = 1:n+2
            ii = (I-1)*n+(i-1)
            jj = (J-1)*n+(j-1)
            if 1<=ii<=NDOF && 1<=jj<=NDOF
                idx = (ii-1)*NDOF + jj
                colext[i,j] = DMAP[idx]
            end
        end
    end
    colext = colext .- 1

    ncolsize = Int64[]
    for i = 2:n+1
        for j = 2:n+1
            ns = (colext[i-1,j]>=0) + (colext[i+1,j]>=0) +
                    (colext[i,j-1]>=0) + (colext[i,j+1]>=0) + (colext[i,j]>=0)
            push!(ncolsize, ns)
        end
    end

    ilower = ((I-1)*N + J-1)*n^2
    iupper = ((I-1)*N + J)*n^2-1

    rows = Int32.(Array(0:n^2 .- 1))
    rows = rows .+ ilower

    return colext, sum(ncolsize), Int32.(rows), Int32.(ncolsize), ilower, iupper
end

function get_poisson_matrix(kext,mc; deps = missing)
    colext, colssize, rows, ncols, ilower, iupper = get_colext(mc)
    if ismissing(deps)
        deps = kext[1,1]
    end
    get_poisson_matrix_ = load_op_and_grad("$(@__DIR__)/build/libGetPoissonMatrix","get_poisson_matrix", multiple=true)
    kext,colext,colssize,deps = convert_to_tensor(Any[kext,colext,colssize,deps], [Float64,Int64,Int64,Float64])
    global cols, vals = get_poisson_matrix_(kext,colext,colssize,deps)
    mpi_SparseTensor(rows, ncols, cols, vals, ilower, iupper, (mc.n*mc.N)^2)
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module StirlingNumbers

using Nemo, IterTools, Triangles

export ModuleStirlingNumbers
export I132393, L132393, V132393, M132393
export I048993, L048993, V048993, M048993
export I271703, L271703, V271703, M271703
export I094587, L094587, V094587, M094587
export I008279, L008279, V008279, M008279
# export Lah, StirlingSet, StirlingCycle

"""
* I132393, L132393, V132393, M132393, I048993, L048993, V048993, M048993, I271703, L271703, V271703, M271703, I094587, L094587, V094587, M094587, I008279, L008279, V008279, M008279
"""
const ModuleStirlingNumbers = ""

"""
Recurrence for A132393, StirlingCycle numbers.
"""
function R132393(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    (n - 1)*prevrow(k) + prevrow(k - 1)
end
"""
Recurrence for A048993, StirlingSet numbers.
"""
function R048993(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    k*prevrow(k) + prevrow(k - 1)
end
"""
Recurrence for A271703, Lah numbers.
"""
function R271703(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    (k - 1 + n)*prevrow(k) + prevrow(k - 1)
end
"""
Recurrence for A094587.
"""
function R094587(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    (n - k)*prevrow(k) + prevrow(k - 1)
end
"""
Recurrence for A008279. Number of permutations of n things k at a time.
"""
function R008279(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    prevrow(k) + k*prevrow(k - 1)
end
"""
Iterates over the first n rows of `A132393`.
"""
I132393(n) = RecTriangle(n, R132393)
"""
Iterates over the first n rows of `A048993`.
"""
I048993(n) = RecTriangle(n, R048993)
"""
Iterates over the first n rows of `A271703`.
"""
I271703(n) = RecTriangle(n, R271703)
"""
Iterates over the first n rows of `A094587`.
"""
I094587(n) = RecTriangle(n, R094587)
"""
Iterates over the first n rows of `A008279`.
"""
I008279(n) = RecTriangle(n, R008279)
"""
Lists the first n rows of `A132393` by concatenating.
"""
L132393(n) = vcat(I132393(n)...)
"""
Lists the first n rows of `A048993` by concatenating.
"""
L048993(n) = vcat(I048993(n)...)
"""
Lists the first n rows of `A271703` by concatenating.
"""
L271703(n) = vcat(I271703(n)...)
"""
Lists the first n rows of `A094587` by concatenating.
"""
L094587(n) = vcat(I094587(n)...)
"""
Lists the first n rows of `A008279` by concatenating.
"""
L008279(n) = vcat(I008279(n)...)
"""
Return the triangular array as a square matrix with dim rows.
"""
M132393(dim) = fromΔ(I132393(dim))
"""
Return the triangular array as a square matrix with dim rows.
"""
M048993(dim) = fromΔ(I048993(dim))
"""
Return the triangular array as a square matrix with dim rows.
"""
M271703(dim) = fromΔ(I271703(dim))
"""
Return the triangular array as a square matrix with dim rows.
"""
M094587(dim) = fromΔ(I094587(dim))
"""
Return the triangular array as a square matrix with dim rows.
"""
M008279(dim) = fromΔ(I008279(dim))
"""
Return row n of A132393 based on the iteration `I132393`(n).
"""
V132393(n) = nth(I132393(n+1), n+1)
"""
Return row n of A048993 based on the iteration `I048993`(n).
"""
V048993(n) = nth(I048993(n+1), n+1)
"""
Return row n of A271703 based on the iteration `I271703`(n).
"""
V271703(n) = nth(I271703(n+1), n+1)
"""
Return row n of A094587 based on the iteration `I094587`(n).
"""
V094587(n) = nth(I094587(n+1), n+1)
"""
Return row n of A008279 based on the iteration `I008279`(n).
"""
V008279(n) = nth(I008279(n+1), n+1)

#START-TEST-########################################################

using Test, Combinatorics

function test()

    @testset "Stirling1" begin

        for (n, S1) in enumerate(I132393(10))
            SC = [stirlings1(n-1, k) for k in 0:n-1]
            @test S1 == SC
        end
    end

    @testset "Stirling2" begin

        for (n, S2) in enumerate(I048993(10))
            SC = [stirlings2(n-1, k) for k in 0:n-1]
            @test S2 == SC
        end
    end
end # test

function demo()
    for row in I132393(9) println(row) end
    println()

    for n in 0:8 println(n, ": ", V132393(n)) end
    println()

    ShowAsΔ(I132393(9), "\t")

    println(L132393(9))
    println()

    ShowAsMatrix(I132393(9))
    println()

    println(M132393(9))
    println()

    println(V132393(8))
    println()

    # ===================================

    for row in I048993(9) println(row) end
    println()

    for n in 0:8 println(n, ": ", V048993(n)) end
    println()

    ShowAsΔ(I048993(9), "\t")

    println(L048993(9))
    println()

    ShowAsMatrix(I048993(9))
    println()

    println(M048993(9))
    println()

    println(V048993(8))
    println()

    # ===================================

    for row in I008279(9) println(row) end
    println()

    for n in 0:8 println(n, ": ", V008279(n)) end
    println()

    ShowAsΔ(I008279(9), "\t")

    println(L008279(9))
    println()

    ShowAsMatrix(I008279(9))
    println()

    println(M008279(9))
    println()

    println(V008279(8))
    println()
end

"""
(for row in I132393(500) end)
    0.246288 seconds (379.25 k allocations: 8.748 MiB)
"""
function perf()
    GC.gc()
    @time (for row in I132393(500) end)
    GC.gc()
    @time (for row in I048993(500) end)
    GC.gc()
    @time (for row in I271703(500) end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
1;
0,    1;
0,    1,     1;
0,    2,     3,     1;
0,    6,    11,     6,    1;
0,   24,    50,    35,   10,    1;
0,  120,   274,   225,   85,   15,   1;
0,  720,  1764,  1624,  735,  175,  21,  1;
0, 5040, 13068, 13132, 6769, 1960, 322, 28, 1;

1
0  1
0  1   1
0  1   3     1
0  1   7     6     1
0  1   15    25    10     1
0  1   31    90    65     15     1
0  1   63    301   350    140    21    1
0  1   127   966   1701   1050   266   28   1

1
0  1
0  2       1
0  6       6       1
0  24      36      12      1
0  120     240     120     20      1
0  720     1800    1200    300     30      1
0  5040    15120   12600   4200    630     42    1
0  40320   141120  141120  58800   11760   1176  56   1

1
1       1
2       2       1
6       6       3       1
24      24      12      4       1
120     120     60      20      5      1
720     720     360     120     30     6     1
5040    5040    2520    840     210    42    7    1
40320   40320   20160   6720    1680   336   56   8    1

1
1  1
1  2  2
1  3  6   6
1  4  12  24   24
1  5  20  60   120  120
1  6  30  120  360  720   720
1  7  42  210  840  2520  5040  5040
1  8  56  336  1680 6720  20160 40320  40320
=#