text:
init_var(s::Symbol) = Symbol(string(s, "00"))

struct RecSystem
    vars::Vector{Symbol}
    params::Vector{SymOrNum}
    init::Matrix{<:Number}
    body::Matrix{<:Number}
end

Base.size(s::RecSystem) = length(s.vars)

function value(l::RecSystem, k::Int)
    @assert k >= 0
    iszero(k) && return tuple(l.init...)
    tuple((l.body^k * l.init)...)
end

value(l::RecSystem, r::UnitRange{Int}) = [value(l, k) for k in r]

function sequentialize(M::Matrix, v::Vector)
    LinearAlgebra.istriu(M) && return M, v
    h = size(M, 2)
    V = [mkvar("_" * string(x)) for x in v]
    V = [V; v]
    tl = zeros(Int, h, h)
    tr = UniformScaling(1)
    bl = LowerTriangular(M) - Diagonal(M)
    br = M - bl
    B = vcat(hcat(tl, tr), hcat(bl, br))
    zcols = findall(iszero, collect(eachcol(B)))
    B[setdiff(1:end, zcols), setdiff(1:end, zcols)], V[setdiff(1:end, zcols)]
end

function loop(l::RecSystem)
    body, vars = sequentialize(l.body, map(mkvar, l.vars))
    lhss = (Meta.parse ∘ string).(body * vars)
    lhss = [replace(x, CONST_ONE_SYM, one(Int)) for x in lhss]
    rhss = [replace(x, CONST_ONE_SYM, one(Int)) for x in l.vars]
    pars = map(x->(x isa Symbol ? mkpoly(init_var(x)) : mkpoly(x)), l.params)
    init = [:($rhs = $lhs) for (rhs,lhs) in zip(rhss, l.init*pars) if rhs != one(Int)]
    assign = [:($rhs = $lhs) for (rhs,lhs) in zip(rhss, lhss) if rhs != lhs]
    striplines(quote
        $(init...)
        while true
            $(assign...)
        end
    end)
end

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

const CONST_ONE_SYM = :_c
const ROOT_SYM = :_w

struct RecurrenceTemplate
    vars::Vector{Symbol}
    params::Vector{SymOrNum}
    body::Matrix{<:Poly}
    init::Matrix{<:Poly}
end

RecurrenceTemplate(r::RecTemplate, vars::Vector{Symbol}) =
    RecurrenceTemplate(vars, r.params, r.body, r.init)

Base.size(rt::RecurrenceTemplate) = length(rt.vars)

variables_init(r::RecurrenceTemplate) =
    [Symbol(string(v)) for v in r.init if !isconstant(v)]

variables_body(r::RecurrenceTemplate) =
    [Symbol(string(v)) for v in r.body if !isconstant(v)]

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

struct ClosedFormTemplate
    vars::Vector{Symbol}
    params::Vector{SymOrNum}
    multiplicities::Vector{Int}

    function ClosedFormTemplate(v::Vector{Symbol}, p::Vector{SymOrNum}, m::Vector{Int})
        @assert length(v) == sum(m) "The sum of multiplicites has to match the number of variables. Got $(sum(m)), need $(length(v))"
        new(v, p, m)
    end
end

ClosedFormTemplate(rt::RecurrenceTemplate, ms::Vector{Int}) = ClosedFormTemplate(rt.vars, rt.params, ms)

function closed_form(ct::ClosedFormTemplate, lc::Symbol, var::Symbol, arg)
    ms = ct.multiplicities
    rs = symroot(length(ms))
    idx = findfirst(x->x==var, ct.vars)
    nroots = length(rs)
    parg = mkpoly(arg)
    pairs = [replacement_pair(NExp{lc}(rs[i], parg)) for i in 1:nroots]
    params = map(mkpoly, ct.params)
    poly = sum(coeff_var(i, j, idx, params=params) * first(pairs[i]) * parg^(j-1) for i in 1:nroots for j in 1:ms[i])
    CFiniteExpr{lc}(poly, Dict(pairs))
end

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

coeff_vec(i::Int, j::Int, rows::Int; params::Vector{<:Poly}) =
    [coeff_var(i, j, k, params=params) for k in 1:rows]

coeff_var(i::Int, j::Int, varidx::Int; params::Vector{<:Poly}) =
    sum(mkvar("c$i$j$(varidx)$l") * p for (l,p) in enumerate(params))

symroot(n::Int) = [mkvar(string(ROOT_SYM, i)) for i in 1:n]

cforms(varcnt::Int, rs::Vector{<:Var}, ms::Vector{Int}; lc::Union{Int,Var}, exp::Int, params::Vector{<:Poly}) =
    sum(coeff_vec(i, j, varcnt, params=params) * rs[i]^exp * lc^(j-1) for i in 1:length(rs) for j in 1:ms[i])

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

struct SynthesisProblem
    inv::Invariant
    rt::RecurrenceTemplate
    ct::ClosedFormTemplate

    function SynthesisProblem(inv::Invariant, rt::RecurrenceTemplate, ct::ClosedFormTemplate)
        @assert issubset(program_variables(inv), rt.vars) "$(program_variables(inv)) is not a subset of $(rt.vars)"
        @assert rt.vars == ct.vars && rt.params == ct.params
        new(inv, rt, ct)
    end
end

vars(s::SynthesisProblem) = s.rt.vars
params(s::SynthesisProblem) = s.rt.params
params(::Type{Poly}, s::SynthesisProblem) = map(mkpoly, params(s))
params(::Type{Symbol}, s::SynthesisProblem) = filter(x->x isa Symbol, params(s))
params(::Type{Var}, s::SynthesisProblem) = map(mkvar, params(Symbol, s))

body(s::SynthesisProblem) = s.rt.body
init(s::SynthesisProblem) = s.rt.init

multiplicities(s::SynthesisProblem) = s.ct.multiplicities
roots(s::SynthesisProblem) = symroot(length(multiplicities(s)))

function constraints(sp::SynthesisProblem; progress::Bool=true)
    cscforms = cstr_cforms(sp)
    csinit   = cstr_init(sp)
    csroots  = cstr_roots(sp)
    csrel    = cstr_algrel(sp)
    @debug "Constraints" cscforms csinit csroots csrel
    pcp = csroots & cscforms & csinit & csrel
    if progress
        pcp &= cstr_progress(sp)
    end
    pcp
end

function create_solver(sp::SynthesisProblem, T::Type{<:NLSolver}; progress::Bool=true)
    pcp = constraints(sp; progress=progress)
    vars = NLSat.variables(pcp)
    varmap = convert(Dict{Symbol,Type}, Dict(v=>AlgebraicNumber for v in vars))
    @debug "Variables" varmap

    vs = variables_body(sp.rt)
    soft = ClauseSet([Clause(Constraint{EQ}(v)) for v in vs])

    solver = T()
    NLSat.add_vars!(solver, varmap)
    NLSat.add!(solver, pcp)
    # NLSat.add_soft!(solver, soft)
    solver
end

function parse_model(sp::SynthesisProblem, model::Dict{Symbol,T}) where {T<:Number}
    _A, _B = init(sp), body(sp)
    A = T[get(model, Symbol(string(b)), b) for b in _A]
    B = T[get(model, Symbol(string(b)), b) for b in _B]
    RecSystem(vars(sp), params(sp), A, B)
end

function solve(sp::SynthesisProblem, solver::S; timeout::Int) where {S<:NLSolver}
    status, elapsed, model = NLSat.solve(solver, timeout = timeout)
    @debug "Solver result" status model
    if status == NLSat.sat
        return model, SynthesisResult(status, parse_model(sp, model), elapsed, sp)
    end
    nothing, SynthesisResult(status, nothing, elapsed, sp)
end

function solve(sp::SynthesisProblem; solver::Type{<:NLSolver}=Z3Solver, progress::Bool=true, timeout::Int=10)
    _solver = create_solver(sp, solver; progress=progress)
    last(solve(sp, _solver, timeout=timeout))
end

"Generate constraints ensuring that p is an algebraic relation."
function cstr_algrel(sp::SynthesisProblem)
    res = constraint_walk(sp.inv) do poly
        expr = function_walk(poly) do func, args
            @assert length(args) == 1 "Invariant not properly preprocessed"
            :($(closed_form(sp.ct, sp.inv.lc, func, args[1])))
        end
        expr = geometric_walk(expr) do base, exp
            @assert exp == sp.inv.lc
            parg = mkpoly(exp)
            pair = replacement_pair(NExp{exp}(mkpoly(base), parg))
            CFiniteExpr{exp}(pair[1], Dict(pair))
        end
        cfin = eval(expr)
        splitvars = Symbol[params(Symbol, sp); sp.inv.lc]
        cstr = cfinite_constraints(cfin; split_vars=splitvars)
        :($cstr)
    end
    eval(res)
end

"Generate constraints ensuring that the root symbols are roots of the characteristic polynomial of B."
function cstr_roots(sp::SynthesisProblem)
    B, rs, ms = body(sp), roots(sp), multiplicities(sp)
    λ = mkvar(gensym_unhashed(:x))
    BB = -B
    for i in diagind(B)
        BB[i] = λ + BB[i]
    end
    cpoly = det(BB)
    factors = prod((λ - r)^m for (r, m) in zip(rs,ms))
    ps = destructpoly(cpoly - factors, λ)
    cs1 = ClauseSet(map(Clause ∘ Constraint{EQ}, ps))

    ps = [r1-r2 for (r1,r2) in combinations(rs, 2)]
    cs2 = ClauseSet(map(Clause ∘ Constraint{NEQ}, ps))

    cs3 = ClauseSet([Clause(Constraint{NEQ}(r)) for r in rs])

    cs1 & cs2 & cs3
end

"Generate constraints defining the initial values."
function cstr_init(sp::SynthesisProblem)
    B, rs, ms = body(sp), roots(sp), multiplicities(sp)
    pars = params(Poly, sp)
    s = size(B, 1)
    d = sum(ms)
    A = init(sp) * pars
    cstr = Poly[]
    for i in 0:d-1
        M = cforms(s, rs, ms, lc=i, exp=i, params=pars)
        X = iszero(i) ? A : B^i*A
        append!(cstr, X - M)
    end
    ps = destructpoly(cstr, params(Var, sp))
    ClauseSet(map(Clause ∘ Constraint{EQ}, ps))
end

"Generate constraints describing the relationship between entries of B and the closed forms."
function cstr_cforms(sp::SynthesisProblem)
    B, rs, ms = body(sp), roots(sp), multiplicities(sp)
    pars = params(Poly, sp)
    t = length(ms)
    rows = size(B, 1)
    Ds = [sum(binomial(k-1, j-1) * coeff_vec(i, k, rows, params=pars) * rs[i] for k in j:ms[i]) - B * coeff_vec(i, j, rows, params=pars) for i in 1:t for j in 1:ms[i]]
    ps = destructpoly(collect(Iterators.flatten(Ds)), params(Var, sp))
    ClauseSet(map(Clause ∘ Constraint{EQ}, ps))
end

"Generate constraints ensuring that sequence is not constant, i.e. B*B*x0 != B*x0."
function cstr_progress(sp::SynthesisProblem)
    pars = params(Poly, sp)
    A, B = init(sp)*pars, body(sp)
    cs1 = B * B * A - B * A
    cs2 = B * B * B * A - B * A
    _cstr_progress(sp, cs1, pars) & _cstr_progress(sp, cs2, pars)
end

function _cstr_progress(sp::SynthesisProblem, cs, pars)
    # do not consider variables which only occurr as initial variable in polys
    invvars = program_variables(sp.inv)
    # filter!(!isinitvar, vars)
    nonloopvars = setdiff(vars(sp), invvars)
    inds = [i for (i, v) in enumerate(vars(sp)) if !(v in invvars)]
    deleteat!(cs, inds)
    res = Clause()
    for c in cs
        ps = destructpoly([c], params(Var, sp))
        res |= Clause(map(Constraint{NEQ}, ps))
    end
    ClauseSet(res)
end

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

struct SynthesisResult
    status::NLStatus
    recsystem::Union{RecSystem,Nothing}
    elapsed::TimePeriod
    problem::SynthesisProblem
end

issat(s::SynthesisResult) = s.status == NLSat.sat

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

function Base.summary(io::IO, r::Union{RecurrenceTemplate,RecSystem})
    compact = get(io, :compact, false)
    if compact
        print(io, "$(typeof(r)) ($(size(r)))")
    else
        print(io, "$(typeof(r)) of size $(size(r))")
    end
end

function Base.show(io::IO, r::Union{RecurrenceTemplate,RecSystem})
    summary(io, r)
    compact = get(io, :compact, false)
    if !compact
        println(io, ":")
        _show(io, r)
    end
end

function _show(io::IO, r::Union{RecurrenceTemplate,RecSystem})
    init = r.init * map(mkpoly, r.params)
    vars = [v == CONST_ONE_SYM ? "1" : "$v" for v in r.vars]
    _print_recsystem(io, vars, r.body, init)
end

function Base.show(io::IO, s::SynthesisResult)
    compact = get(io, :compact, false)
    if compact
        print(io, "SynthesisResult($(s.status), $(s.elapsed))")
    else
        print(io, "SynthesisResult ($(s.status)) in $(s.elapsed)")
        if !isnothing(s.recsystem)
            println(io)
            _show(io, s.recsystem)
        end
    end
end


text:
rootdir ="/mnt/e/xxx号机4月4日~5日数据/c"

## ###############所有指定时间段内的数据
startTime= "20190304" * "200109"
endTime=  "20190304" * "231544"
## ###############所有时间全部分析
startTime =endTime =""
##########

println()
# fileRawDataChannel= Channel(300)
return (rootdir,startTime,endTime)

# function  writeTitle(CumSumDataobj)
# end

ext: jl

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# Similar to JuMP/test/solvers.jl

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

mos = try_import(:Mosek)
scs = try_import(:SCS)
csdp = try_import(:CSDP)

isscs(solver) = contains(string(typeof(solver)),"SCSSolver")

# Semidefinite solvers
sdp_solvers = Any[]
# Need at least Mosek 8 for sosdemo3 to pass; see:
# https://github.com/JuliaOpt/Mosek.jl/issues/92
# and at least Mosek 8.0.0.41 for sosdemo6 to pass; see:
# https://github.com/JuliaOpt/Mosek.jl/issues/98
mos && push!(sdp_solvers, Mosek.MosekSolver(LOG=0))
# Need 54000 iterations for sosdemo3 to pass on Linux 64 bits
# With 55000, sosdemo3 passes for every platform except Windows 64 bits on AppVeyor
scs && push!(sdp_solvers, SCS.SCSSolver(eps=1e-6, max_iters=60000, verbose=0))

csdp && push!(sdp_solvers, CSDP.CSDPSolver(printlevel=1, maxiter=60000, axtol=1.0e-8, atytol=1.0e-8))

max_stars_repo_name: polislizarralde/jsos

max_issues_repo_name: polislizarralde/jsos

max_forks_repo_name: polislizarralde/jsos


text:
using ForwardDiff

println("Testing Partials...")
tic()
include("PartialsTest.jl")
println("done (took $(toq()) seconds).")

println("Testing Dual...")
tic()
include("DualTest.jl")
println("done (took $(toq()) seconds).")

println("Testing derivative functionality...")
tic()
include("DerivativeTest.jl")
println("done (took $(toq()) seconds).")

println("Testing gradient functionality...")
tic()
include("GradientTest.jl")
println("done (took $(toq()) seconds).")

println("Testing jacobian functionality...")
tic()
include("JacobianTest.jl")
println("done (took $(toq()) seconds).")

println("Testing hessian functionality...")
tic()
include("HessianTest.jl")
println("done (took $(toq()) seconds).")

println("Testing perturbation confusion functionality...")
tic()
include("ConfusionTest.jl")
println("done (took $(toq()) seconds).")

println("Testing miscellaneous functionality...")
tic()
include("MiscTest.jl")
println("done (took $(toq()) seconds).")

if Base.JLOptions().opt_level >= 3
    println("Testing SIMD vectorization...")
    tic()
    include("SIMDTest.jl")
    println("done (took $(toq()) seconds).")
end

println("Testing deprecations...")
tic()
include("DeprecatedTest.jl")
println("done (took $(toq()) seconds).")


text:
export ccat,
	squeeze,
    all_αs,
    all_ws_sansbn,
    flip_batch!,
    shift_batch!,
	cutout_batch!,
	norm_batch!,
    TrainCuIterator,
    EvalCuIterator,
    TestCuIterator,
    CosineAnnealing,
	apply!,
	depthwiseconv!,
	cudnnDepthwiseConvolutionDescriptor
using Adapt
using CUDA
using Flux
import Flux.Optimise.apply!

ccat(X...) = cat(X...; dims=3)

function squeeze(A::AbstractArray) #generalize this?
    if ndims(A) == 3
        if size(A, 3) > 1
            return dropdims(A; dims = (1))
        elseif size(A, 3) == 1
            return dropdims(A; dims = (1, 3))
        end
    elseif ndims(A) == 4
        if size(A, 4) > 1
            return dropdims(A; dims = (1, 2))
        elseif size(A, 4) == 1
            return dropdims(A; dims = (1, 2, 4))
        end
    end
    return A
end

function all_params(submodels)
    ps = Params()
    for submodel in submodels
        Flux.params!(ps, submodel)
    end
    return ps
end

all_αs(model) = Flux.params([model.normal_αs, model.reduce_αs])
all_ws(model) = Flux.params([model.stem, model.cells..., model.global_pooling, model.classifier])

function all_ws_sansbn(model) #without batchnorm params
    all_w = Flux.params([model.stem, model.cells..., model.global_pooling, model.classifier])
    for (i,cell) in enumerate(model.cells)
        for (j,mixedop) in enumerate(cell.mixedops)
            for (k,op) in enumerate(mixedop.ops)
                for (l,layer) in enumerate(op.op)
                    if typeof(layer) <: Flux.BatchNorm
                        delete!(all_w, layer.γ)
                        delete!(all_w, layer.β)
                    end
                end
            end
        end
    end
    for (l,layer) in enumerate(model.stem.layers)
        if typeof(layer) <: Flux.BatchNorm
            delete!(all_w, layer.γ)
            delete!(all_w, layer.β)
        end
    end
    all_w
end

CIFAR_MEAN = [0.49139968f0, 0.48215827f0, 0.44653124f0]
CIFAR_STD = [0.24703233f0, 0.24348505f0, 0.26158768f0]

function flip_batch!(batch::Array{Float32,4})
	flips = falses(size(batch,4))
	for image in 1:size(batch,4)
		orig = copy(batch[:,:,:,image])
		flip = rand(Bool)
		if flip
			flipped = reverse(orig, dims=2)
			batch[:,:,:,image] = flipped
		end
		flips[image] = flip
	end
	flips
end

function shift_batch!(batch::Array{Float32,4})
	shifts = Array{Int64}(undef,size(batch,4),2)
	for image in 1:size(batch,4)
		orig = copy(batch[:,:,:,image])
		shiftx = rand(-4:4)
		shifty = rand(-4:4)
		if shiftx > 0
			batch[1:size(batch,1)-shiftx,:,:,image] = orig[shiftx+1:size(batch,1),:,:]
			batch[size(batch,1)-shiftx+1:size(batch,1),:,:,image] .= 0f0
		elseif shiftx < 0
			batch[1:-shiftx,:,:,image] .= 0f0
			batch[1-shiftx:size(batch,1),:,:,image] = orig[1:size(batch,1)+shiftx,:,:]
		end
		orig = copy(batch[:,:,:,image])
		if shifty > 0
			batch[:,1:size(batch,2)-shifty,:,image] = orig[:,shifty+1:size(batch,2),:]
			batch[:,size(batch,2)-shifty+1:size(batch,2),:,image] .= 0f0
		elseif shifty < 0
			batch[:,1:-shifty,:,image] .= 0f0
			batch[:,1-shifty:size(batch,2),:,image] = orig[:,1:size(batch,2)+shifty,:]
		end
		shifts[image,:] = [shiftx;shifty]
	end
	shifts
end

function cutout_batch!(batch::Array{Float32,4}, cutout::Int = -1)
	cutouts = Array{Int64}(undef,size(batch,4),2)
	if cutout > 0
		for image in 1:size(batch,4)
			cutx = rand(1:size(batch,1))
			cuty = rand(1:size(batch,2))
			minx = maximum([cutx-cutout÷2,1])
			maxx = minimum([cutx+cutout÷2-1,size(batch,1)])
			miny = maximum([cuty-cutout÷2,1])
			maxy = minimum([cuty+cutout÷2-1,size(batch,2)])
			batch[minx:maxx,miny:maxy,:,image] .= 0f0
			cutouts[image,:] = [cutx;cuty]
		end
	end
end

function norm_batch!(batch::Array{Float32,4})
	mean_im = repeat(reshape(CIFAR_MEAN, (1,1,3,1)), outer = [32,32,1,size(batch,4)])
	std_im = repeat(reshape(CIFAR_STD, (1,1,3,1)), outer = [32,32,1,size(batch,4)])
	batch = (batch.-mean_im)./std_im
end

mutable struct TrainCuIterator{B}
    batches::B
    previous::Any
    TrainCuIterator(batches) = new{typeof(batches)}(batches)
end
function Base.iterate(c::TrainCuIterator, state...)
    item = iterate(c.batches, state...)
    isdefined(c, :previous) && foreach(CUDA.unsafe_free!, c.previous)
    item === nothing && return nothing
    batch, next_state = item
	flip_batch!(batch[1])
	shift_batch!(batch[1])
	norm_batch!(batch[1])
    cubatch = map(x -> adapt(CuArray, x), batch)
    c.previous = cubatch
	return cubatch, next_state
end

mutable struct EvalCuIterator{B}
    batches::B
    previous::Any
    EvalCuIterator(batches) = new{typeof(batches)}(batches)
end
function Base.iterate(c::EvalCuIterator, state...)
    item = iterate(c.batches, state...)
    isdefined(c, :previous) && foreach(CUDA.unsafe_free!, c.previous)
    item === nothing && return nothing
    batch, next_state = item
	flip_batch!(batch[1])
	shift_batch!(batch[1])
	norm_batch!(batch[1])
	cutout_batch!(batch[1], 16)
    cubatch = map(x -> adapt(CuArray, x), batch)
    c.previous = cubatch
    return cubatch, next_state
end

mutable struct TestCuIterator{B}
    batches::B
    previous::Any
    TestCuIterator(batches) = new{typeof(batches)}(batches)
end
function Base.iterate(c::TestCuIterator, state...)
    item = iterate(c.batches, state...)
    isdefined(c, :previous) && foreach(CUDA.unsafe_free!, c.previous)
    item === nothing && return nothing
    batch, next_state = item
	norm_batch!(batch[1])
    cubatch = map(x -> adapt(CuArray, x), batch)
    c.previous = cubatch
    return cubatch, next_state
end

mutable struct CosineAnnealing
  tmax::Int64
  t::Int64
end

CosineAnnealing(tmax::Int64 = 1) = CosineAnnealing(tmax, 0)

function Flux.Optimise.apply!(o::CosineAnnealing, x, Δ)
  tmax = o.tmax
  t = o.t
  Δ .*= (1 + cos(t/tmax*pi))/2
  return Δ
end

import NNlib: depthwiseconv!, ∇depthwiseconv_data!, ∇depthwiseconv_filter!

using CUDA.APIUtils: @workspace

using CUDA.CUDNN:
    cudnnConvolutionForward,
    cudnnConvolutionForward!,
    cudnnConvolutionBackwardFilter,
    cudnnConvolutionBackwardData,
    cudnnGetConvolutionNdForwardOutputDim,
    cudnnSetConvolutionMathType,
    cudnnSetConvolutionReorderType,
    cudnnSetConvolutionGroupCount,
    cudnnFindConvolutionForwardAlgorithmEx,
        cudnnConvolutionFwdAlgoPerf_t,
    cudnnFindConvolutionBackwardFilterAlgorithmEx,
        cudnnConvolutionBwdFilterAlgoPerf_t,
    cudnnFindConvolutionBackwardDataAlgorithmEx,
        cudnnConvolutionBwdDataAlgoPerf_t,
    cudnnConvolutionDescriptor,
        cudnnConvolutionDescriptor_t,
        cudnnCreateConvolutionDescriptor,
        cudnnSetConvolutionNdDescriptor,
        cudnnDestroyConvolutionDescriptor,
    cudnnConvolutionMode_t,
        CUDNN_CONVOLUTION,       # 0
        CUDNN_CROSS_CORRELATION, # 1
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,      # 0
        CUDNN_ACTIVATION_RELU,         # 1
        CUDNN_ACTIVATION_TANH,         # 2
        CUDNN_ACTIVATION_CLIPPED_RELU, # 3
        CUDNN_ACTIVATION_ELU,          # 4
        CUDNN_ACTIVATION_IDENTITY,     # 5
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnMathType_t,
        CUDNN_DEFAULT_MATH,                    # 0
        CUDNN_TENSOR_OP_MATH,                  # 1
        CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION, # 2
        CUDNN_FMA_MATH,                        # 3
    cudnnReorderType_t,
        CUDNN_DEFAULT_REORDER, # 0
        CUDNN_NO_REORDER,      # 1
    cudnnConvolutionFwdAlgo_t,
        CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM,         # 0
        CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM, # 1
        CUDNN_CONVOLUTION_FWD_ALGO_GEMM,                  # 2
        CUDNN_CONVOLUTION_FWD_ALGO_DIRECT,                # 3
        CUDNN_CONVOLUTION_FWD_ALGO_FFT,                   # 4
        CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING,            # 5
        CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD,              # 6
        CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD_NONFUSED,     # 7
        CUDNN_CONVOLUTION_FWD_ALGO_COUNT,                 # 8
    cudnnConvolutionBwdFilterAlgo_t,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0,                 # 0, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1,                 # 1,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT,               # 2,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3,                 # 3, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD,          # 4, /* not implemented */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD_NONFUSED, # 5,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT_TILING,        # 6,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_COUNT,             # 7
    cudnnConvolutionBwdDataAlgo_t,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_0,                 # 0, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_1,                 # 1,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT,               # 2,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING,        # 3,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD,          # 4,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD_NONFUSED, # 5,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_COUNT,             # 6
    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,        # 0, /* row major (wStride = 1, hStride = w) */
        CUDNN_TENSOR_NHWC,        # 1, /* feature maps interleaved ( cStride = 1 )*/
        CUDNN_TENSOR_NCHW_VECT_C, # 2, /* each image point is vector of element of C, vector length in data type */
    cudnnDataType,
    convdims,
    math_mode,
    handle,
	scalingParameter,
	cudnnTensorDescriptor,
	cudnnFilterDescriptor,
	cudnnConvolutionBwdDataAlgoPerf,
	cudnnConvolutionBwdFilterAlgoPerf,
	cudnnConvolutionBackwardData,
	cudnnConvolutionBackwardFilter

function nnlibPadding(dims)
    pd = NNlib.padding(dims)
    if !all(pd[1:2:end] .== pd[2:2:end])
        @warn "cuDNN does not support asymmetric padding; defaulting to symmetric choice" maxlog=1
    end
    return pd[1:2:end]
end
const CUDNNFloat = Union{Float16,Float32,Float64}

function cudnnDepthwiseConvolutionDescriptor(cdims::DepthwiseConvDims, x::DenseCuArray{T}) where T
    cudnnConvolutionDescriptor(convdims(nnlibPadding(cdims),size(x)), convdims(NNlib.stride(cdims),size(x)), convdims(NNlib.dilation(cdims),size(x)), CUDNN_CONVOLUTION, cudnnDataType(T), CUDNN_TENSOR_OP_MATH, CUDNN_DEFAULT_REORDER, Cint(size(x,3)))
end

function depthwiseconv!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T}, cdims::DepthwiseConvDims;
               alpha=1, beta=0, algo=-1) where T<:Union{Float16,Float32,Float64}
    d = cudnnDepthwiseConvolutionDescriptor(cdims, x)
    cudnnConvolutionForward!(y, w, x, d; alpha, beta, z=y)
end

function ∇depthwiseconv_data!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, w::DenseCuArray{T},
                     cdims::DepthwiseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(dx), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(w)
    convDesc = cudnnDepthwiseConvolutionDescriptor(cdims, dx)
    p = cudnnConvolutionBwdDataAlgoPerf(wDesc, w, yDesc, dy, convDesc, xDesc, dx)
    @workspace size=p.memory workspace->cudnnConvolutionBackwardData(handle(), alpha, wDesc, w, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, xDesc, dx)
    return dx
end

function ∇depthwiseconv_filter!(dw::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                       cdims::DepthwiseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(x), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(dw)
    convDesc = cudnnDepthwiseConvolutionDescriptor(cdims, x)
    p = cudnnConvolutionBwdFilterAlgoPerf(xDesc, x, yDesc, dy, convDesc, wDesc, dw);
    @workspace size=p.memory workspace->cudnnConvolutionBackwardFilter(handle(), alpha, xDesc, x, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, wDesc, dw);
    return dw
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
@testset "attributes" begin
    l = 10
    I = [1, missing, 3, missing, 5]

    for i in I
        mhv = MaybeHotVector(i, l)
        @test size(mhv) == (l,)
        @test length(mhv) == l
    end

    mhm = MaybeHotMatrix(I, l)
    @test size(mhm) == (l, length(I))
    @test length(mhm) == l * length(I)
end

@testset "type construction" begin
    @test MaybeHotVector(1, 10) isa AbstractVector{Bool}
    @test MaybeHotVector(missing, 10) isa AbstractVector{Missing}
    @test MaybeHotMatrix([1, 2], 10) isa AbstractMatrix{Bool}
    @test MaybeHotMatrix([missing, missing], 10) isa AbstractMatrix{Missing}
    @test MaybeHotMatrix([1, missing], 10) isa AbstractMatrix{Union{Bool, Missing}}
end

@testset "comparisons" begin
    @test MaybeHotVector(2, 10) == MaybeHotVector(2, 10)
    @test isequal(MaybeHotVector(1, 10), MaybeHotVector(1, 10))

    @test isequal(MaybeHotVector(missing, 10), MaybeHotVector(missing, 10))

    @test MaybeHotMatrix([7,2], 10) == MaybeHotMatrix([7,2], 10)
    @test isequal(MaybeHotMatrix([1,3,9], 10), MaybeHotMatrix([1,3,9], 10))

    @test isequal(MaybeHotMatrix([1,2,missing], 10), MaybeHotMatrix([1,2,missing], 10))
    @test !isequal(MaybeHotMatrix([1,2,missing,3], 10) == MaybeHotMatrix([1,2,missing,4], 10))
end

@testset "hcat" begin
    l = 10
    I = [1, missing, 3, missing, 5]
    mhm = MaybeHotMatrix(I, l)
    mhm2 = MaybeHotMatrix([1,5,3], l)
    mhm3 = MaybeHotMatrix(fill(missing, 4), l)
    mhvs = [MaybeHotVector.(I, l)...]

    @test all(mhv -> areequal(hcat(mhv), reduce(hcat, [mhv]), MaybeHotMatrix(mhv)),
              mhvs)

    @test all(mhm -> areequal(hcat(mhm), reduce(hcat, [mhm]), mhm), [mhm, mhm2, mhm3])

    @test areequal(hcat(mhvs...), reduce(hcat, mhvs), mhm)

    @test areequal(hcat(mhm, mhm),
                   reduce(hcat, [mhm, mhm]),
                   hcat(mhvs..., mhvs...),
                   reduce(hcat, [mhvs..., mhvs...]),
                   MaybeHotMatrix(vcat(I, I), l))

    @test areequal(hcat(mhm, mhm2), reduce(hcat, [mhm, mhm2]), MaybeHotMatrix(vcat(mhm.I, mhm2.I), l))
    @test areequal(hcat(mhm, mhm3), reduce(hcat, [mhm, mhm3]), MaybeHotMatrix(vcat(mhm.I, mhm3.I), l))
    @test areequal(hcat(mhm3, mhm2), reduce(hcat, [mhm3, mhm2]), MaybeHotMatrix(vcat(mhm3.I, mhm2.I), l))

    @inferred hcat(mhm)
    @inferred reduce(hcat, [mhm])
    @inferred hcat(mhm2)
    @inferred reduce(hcat, [mhm2])
    @inferred hcat(mhm3)
    @inferred reduce(hcat, [mhm3])
    @inferred hcat(mhm, mhm2)
    @inferred reduce(hcat, [mhm, mhm2])
    @inferred hcat(mhm, mhm3)
    @inferred reduce(hcat, [mhm, mhm3])
    @inferred hcat(mhm, mhm2, mhm3)
    @inferred reduce(hcat, [mhm, mhm2, mhm3])
    @inferred hcat(mhvs...)
    @inferred reduce(hcat, mhvs)

    @test_throws DimensionMismatch hcat(MaybeHotVector.([1, 2], [l, l+1])...)
    @test_throws DimensionMismatch hcat(MaybeHotMatrix.([[1], [2, 3]], [l, l+1])...)
    @test_throws DimensionMismatch reduce(hcat, MaybeHotVector.([1, 2], [l, l+1]))
    @test_throws DimensionMismatch reduce(hcat, MaybeHotMatrix.([[1], [2, 3]], [l, l+1]))
end

@testset "indexing" begin
    l = 10
    I = [1, missing, 3, missing, 5]

    for i in I
        mhv = MaybeHotVector(i, l)
        if ismissing(i)
            @test all(isequal.(missing, mhv))
        else
            @test Vector(mhv) == onehot(i, 1:l)
        end
        @test_throws BoundsError mhv[0]
        @test_throws BoundsError mhv[l+1]
        @test all(isequal.(mhv[:], mhv))
    end

    mhm = MaybeHotMatrix(I, l)
    m = Matrix(mhm)
    for (k,i) in I |> enumerate
        if ismissing(i)
            @test all(isequal.(missing, m[:, k]))
            @test all(isequal.(missing, mhm[:, k]))
        else
            @test mhm[:, k] == m[:, k] == onehot(i, 1:l)
        end
    end
    @test isequal(mhm[[1,2,7], 3], m[[1,2,7], 3])
    @test isequal(mhm[CartesianIndex(2, 4)], m[2,4])
    for k in eachindex(I)
        @test isequal(mhm[:, k], MaybeHotVector(I[k], l))
    end
    @test isequal(mhm, mhm[:, eachindex(I)])
    @test isequal(mhm, mhm[:, eachindex(I) |> collect])
    @test isequal(mhm, mhm[:, :])
    @test isequal(mhm, hcat(MaybeHotVector.(I, l)...))
    @test isequal(mhm[:, [1,2,5]], hcat(MaybeHotVector.(I[[1,2,5]], l)...))

    @test_throws BoundsError mhm[0, 1]
    @test_throws BoundsError mhm[2, -1]
    @test_throws BoundsError mhm[CartesianIndex()]
    @test_throws BoundsError mhm[CartesianIndex(1)]
end

@testset "multiplication" begin
    W = rand(10, 10)
    x1 = MaybeHotVector(1, 10)
    x2 = MaybeHotVector(8, 10)
    x3 = MaybeHotVector(missing, 10)
    X1 = MaybeHotMatrix([7, 10], 10)
    X2 = MaybeHotMatrix([missing, missing], 10)
    X3 = MaybeHotMatrix([3, missing, 1, missing], 10)

    @test isequal(W * x1, W * Vector(x1))
    @test isequal(W * x2, W * Vector(x2))
    @test isequal(W * x3, W * Vector(x3))
    @test isequal(W * X1, W * Matrix(X1))
    @test isequal(W * X2, W * Matrix(X2))
    @test isequal(W * X3, W * Matrix(X3))

    @test eltype(W * x1) === eltype(W * x2) === eltype(W * X1) === eltype(W)
    @test eltype(W * x3) === eltype(W * X2) === Missing
    @test eltype(W * X3) === Union{Missing, eltype(W)}

    @inferred W * x1
    @inferred W * x2
    @inferred W * x3
    @inferred W * X1
    @inferred W * X2
    @inferred W * X3

    @test_throws DimensionMismatch W * MaybeHotVector(1, 5)
    @test_throws DimensionMismatch W * MaybeHotVector(missing, 3)
    @test_throws DimensionMismatch W * MaybeHotMatrix([1, 2], 9)
    @test_throws DimensionMismatch W * MaybeHotMatrix([1, missing, 2], 9)
    @test_throws DimensionMismatch W * MaybeHotMatrix([missing, missing], 9)
end

@testset "equality" begin
    mhv1 = MaybeHotVector(1, 10)
    mhv2 = MaybeHotVector(1, 10)
    mhv3 = MaybeHotVector(1, 11)
    mhv4 = MaybeHotVector(missing, 11)
    mhv5 = MaybeHotVector(2, 10)
    @test mhv1 == mhv2
    @test mhv1 != mhv3
    @test !isequal(mhv1, mhv4)
    @test mhv1 != mhv5

    mhm1 = MaybeHotMatrix([1,2], 10)
    mhm2 = MaybeHotMatrix([1,2], 10)
    mhm3 = MaybeHotMatrix([1], 10)
    mhm4 = MaybeHotMatrix([missing], 10)
    mhm5 = MaybeHotMatrix([1,2], 11)
    @test mhm1 == mhm2
    @test mhm1 != mhm3
    @test !isequal(mhv1, mhv4)
    @test mhm1 != mhm5
end

@testset "onehot and onehotbatch" begin
    i = 1
    b = MaybeHotVector(i, 10)
    @test onehot(b) == onehot(i, 1:length(b))
    b = MaybeHotVector(missing, 10)
    @test_throws MethodError onehot(b)
    I = [3, 1, 2]
    B = MaybeHotMatrix(I, 10)
    @test onehotbatch(B) == onehotbatch(I, 1:size(B, 1))
    B = MaybeHotMatrix([3, missing, 2], 10)
    @test_throws MethodError onehotbatch(B)
end

@testset "AbstractMatrix * MaybeHotVector{<:Integer} gradtest" begin
    # for MaybeHot types with missing elements, it doesn't make sense to compute gradient
    for (m, n) in product(fill((1, 5, 10, 20), 2)...), i in [rand(1:n) for _ in 1:3]
        A = randn(m, n)
        b = MaybeHotVector(i, n)

        @test gradtest(A -> A * b, A)

        dA, db = gradient(sum ∘ *, A, b)
        @test dA ≈ gradient(A -> sum(A * onehot(b)), A) |> only
        @test isnothing(db)
    end
end

@testset "AbstractMatrix * MaybeHotMatrix{<:Integer} gradtest" begin
    # for MaybeHot types with missing elements, it doesn't make sense to compute gradient
    for (m, n, k) in product(fill((1, 5, 10, 20), 3)...), I in [rand(1:n, k) for _ in 1:3]
        A = randn(m, n)
        B = MaybeHotMatrix(I, n)

        @test gradtest(A -> A * B, A)

        dA, dB = gradient(sum ∘ *, A, B)
        @test dA ≈ gradient(A -> sum(A * onehotbatch(B)), A) |> only
        @test isnothing(dB)
    end
end

@testset "maybehot" begin
    @test_throws ArgumentError maybehot(0, 1:3)
    @test_throws ArgumentError maybehot(4, 1:3)
    @test_throws ArgumentError maybehot("a", [1, 2])

    @test ismissing(maybehot(missing, 1:3).i)
    @test maybehot(2, 1:3).i == 2

    @test maybehot(missing, 1:3) isa AbstractVector{Missing}
    @test maybehot(2, 1:3) isa AbstractVector{Bool}

    @test maybehot(missing, 1:3) isa MaybeHotVector{Missing}
    @test maybehot(2, 1:3) isa MaybeHotVector{UInt32}
end

@testset "maybehotbatch" begin
    @test_throws ArgumentError maybehotbatch([1, 2, 3, 0], 1:3)
    @test_throws ArgumentError maybehotbatch([4], 1:3)
    @test_throws ArgumentError maybehotbatch([2, "a"], [1, 2])

    mhm = maybehotbatch([missing, missing], 1:3)

    @test all(isequal.(mhm.I, [missing, missing]))
    @test mhm isa AbstractMatrix{Missing}
    @test mhm isa MaybeHotMatrix{Missing}

    mhm = maybehotbatch([3, 1], 1:3)

    @test mhm.I == [3, 1]
    @test mhm isa AbstractMatrix{Bool}
    @test mhm isa MaybeHotMatrix{UInt32}

    mhm = maybehotbatch([1, missing], 1:3)

    @test all(isequal.(mhm.I, [1, missing]))
    @test mhm isa AbstractMatrix{Union{Bool, Missing}}
    @test mhm isa MaybeHotMatrix{Union{UInt32, Missing}}
end


text:
using FFTResampling, Documenter

 # set seed fixed for documentation
DocMeta.setdocmeta!(FFTResampling, :DocTestSetup, :(using FFTResampling); recursive=true)
makedocs(modules = [FFTResampling],
         sitename = "FFTResampling.jl",
         pages = Any[
            "FFTResampling.jl" => "index.md",
         ]
        )

deploydocs(repo = "github.com/roflmaostc/FFTResampling.jl.git", devbranch="main")


text:
# ---
# title: 1124. Longest Well-Performing Interval
# id: problem1124
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/longest-well-performing-interval/description/>
# hidden: true
# ---
#
# We are given `hours`, a list of the number of hours worked per day for a given
# employee.
#
# A day is considered to be a _tiring day_ if and only if the number of hours
# worked is (strictly) greater than `8`.
#
# A _well-performing interval_ is an interval of days for which the number of
# tiring days is strictly larger than the number of non-tiring days.
#
# Return the length of the longest well-performing interval.
#
#
#
# **Example 1:**
#
#
#
#     Input: hours = [9,9,6,0,6,6,9]
#     Output: 3
#     Explanation: The longest well-performing interval is [9,9,6].
#
#
#
#
# **Constraints:**
#
#   * `1 <= hours.length <= 10000`
#   * `0 <= hours[i] <= 16`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module TestChildren

include("../src/GraphIdx.jl")

using Test
import .GraphIdx.Tree: ChildrenIndex, node_degrees, root_node, find_root

@testset "generate(width=5): low level   " begin
    # root, parent = TreeInstance.generate(5)[3:4]
    root = 1
    parent = Int[1, 1, 2, 3, 4, 1, 8, 9, 4, 9]
    cidx = ChildrenIndex(parent, root)
    value_, idx_ = cidx.value, cidx.idx
    @test value_ == [1, 2, 6, 3, 4, 5, 9, 7, 8, 10]
    @test idx_ == [2, 4, 5, 6, 8, 8, 8, 8, 9, 11, 11]
    @test root_node(cidx) == root
    @test length(cidx) == length(parent)
end

@testset "generate(width=5): index access" begin
    # root, parent = TreeInstance.generate(5)[3:4]
    root = 1
    parent = Int[1, 1, 2, 3, 4, 1, 8, 9, 4, 9]
    cidx = ChildrenIndex(parent)
    @test cidx[1] == [2,6]
    @test cidx[2] == [3]
    @test cidx[9] == [8, 10]
    @test cidx[4] == [5, 9]
    @test cidx[10] == []
    @test cidx[7] == []
    @test root_node(cidx) == root
    @test length(cidx) == length(parent)
end

@testset "generate(width=5): index access" begin
    parent0 = Int[1, 1, 2, 3, 4, 1, 8, 9, 4, 9]
    root0 = 1
    @test find_root(parent0) == root0
    n = length(parent0)

    # import Random
    # perm = Random.randperm(Random.MersenneTwister(42), 10)
    perm = [9, 3, 1, 5, 10, 7, 2, 6, 4, 8]
    iperm = invperm(perm)
    @test perm[iperm] == 1:n
    @test iperm[perm] == 1:n
    # perm: old  -> new
    # iperm: new -> old
    # parent0: old i -> old parent
    # parent1: new i -> old i -> old parent -> new parent
    parent1 = perm[parent0[iperm]]
    root1 = find_root(parent1)
    @test root1 == perm[root0]

    cidx = ChildrenIndex(parent1)
    @test root_node(cidx) == root1
    @test length(cidx) == length(parent1)

    @test Set(cidx[perm[1]]) == Set(perm[[2,6]])
    @test Set(cidx[perm[2]]) == Set(perm[[3]])
    @test Set(cidx[perm[9]]) == Set(perm[[8, 10]])
    @test Set(cidx[perm[4]]) == Set(perm[[5, 9]])
    @test Set(cidx[perm[10]]) == Set(perm[[]])
    @test Set(cidx[perm[7]]) == Set(perm[[]])
end

@testset "generate(width=21): some tests " begin
    # pi = TreeInstance.generate(21)[4]
    pi = [1, 1, 2, 3, 4, 7, 8, 29, 8, 9, 32, 33, 34, 35, 36, 37, 16,
    19, 40, 21, 42, 1, 22, 25, 4, 25, 26, 27, 28, 29, 10, 31, 34, 55,
    56, 57, 36, 17, 18, 41, 20, 63, 64, 23, 44, 67, 48, 69, 28, 29,
    72, 51, 52, 33, 76, 55, 56, 57, 38, 59, 60, 63, 84, 65, 86, 65,
    66, 47, 70, 49, 70, 73, 74, 95, 74, 75, 98, 99, 78, 59, 80, 81,
    104, 105, 86, 87, 108, 67, 68, 111, 92, 71, 92, 93, 94, 95, 76,
    97, 98, 79, 122, 103, 82, 103, 104, 127, 108, 109, 110, 89, 110,
    91, 114, 93, 114, 115, 138, 117, 98, 119, 120, 121, 144, 103, 146,
    105, 128, 107, 130, 131, 110, 131, 134, 155, 134, 137, 116, 137,
    138, 141, 142, 121, 164, 145, 124, 145, 126, 127, 128, 149, 152,
    153, 132, 153, 154, 155, 156, 137, 138, 159, 140, 163, 142, 163,
    166, 145, 146, 147]

    cidx = ChildrenIndex(pi)
    @test cidx[162] == []
    @test Set(cidx[137]) == Set([138, 136, 158])
    @test Set(cidx[1]) == Set([2, 22])
    @test cidx[89] == [110]
    @test root_node(cidx) == 1
    @test length(cidx) == length(pi)
end

@testset "generate(width=50)             " begin
    # pi =  TreeInstance.generate(50)[4]

    pi = [1, 1, 2, 3, 4, 7, 57, 7, 59, 9, 61, 13, 63, 64, 16, 17, 67,
    68, 18, 21, 71, 21, 73, 23, 24, 25, 28, 78, 28, 31, 81, 82, 32,
    84, 85, 35, 36, 88, 40, 90, 91, 43, 93, 45, 95, 45, 97, 98, 99,
    100, 101, 51, 3, 55, 5, 55, 56, 108, 60, 110, 111, 61, 62, 63, 64,
    116, 117, 67, 68, 69, 70, 71, 72, 73, 76, 26, 127, 128, 78, 79,
    82, 132, 82, 83, 135, 36, 137, 89, 139, 89, 90, 91, 92, 144, 94,
    95, 98, 99, 149, 99, 102, 103, 53, 54, 55, 107, 157, 107, 59, 160,
    110, 62, 112, 115, 165, 117, 167, 168, 118, 70, 71, 72, 73, 74,
    75, 125, 126, 127, 128, 129, 130, 182, 134, 184, 185, 86, 136,
    137, 138, 90, 140, 141, 193, 145, 195, 96, 146, 198, 199, 149,
    152, 153, 103, 104, 105, 155, 207, 159, 209, 159, 162, 112, 162,
    214, 164, 165, 217, 167, 119, 120, 172, 122, 123, 175, 176, 126,
    176, 177, 178, 179, 180, 181, 182, 183, 184, 187, 137, 138, 239,
    140, 141, 191, 192, 144, 245, 246, 198, 248, 198, 199, 151, 203,
    153, 254, 155, 156, 206, 258, 259, 211, 261, 211, 212, 213, 165,
    215, 267, 268, 169, 219, 271, 272, 224, 274, 175, 176, 228, 278,
    179, 180, 232, 182, 283, 233, 236, 186, 187, 288, 238, 241, 242,
    243, 193, 243, 244, 247, 297, 247, 248, 300, 201, 302, 252, 255,
    205, 257, 207, 259, 309, 310, 260, 212, 262, 263, 315, 216, 266,
    269, 219, 271, 321, 273, 223, 324, 276, 326, 327, 328, 329, 230,
    231, 281, 282, 334, 284, 236, 286, 287, 288, 291, 241, 291, 243,
    344, 245, 246, 298, 348, 349, 299, 302, 352, 304, 254, 304, 256,
    257, 309, 359, 309, 261, 311, 312, 313, 316, 317, 367, 368, 269,
    319, 320, 323, 373, 323, 324, 376, 326, 378, 330, 331, 332, 333,
    283, 333, 336, 337, 338, 288, 289, 341, 291, 341, 293, 343, 344,
    345, 346, 347, 348, 400, 401, 353, 403, 304, 405, 306, 307, 357,
    358, 359, 311, 412, 313, 363, 364, 367, 368, 369, 319, 420, 321,
    371, 372, 324, 374, 377, 427, 377, 380, 330, 380, 381, 433, 385,
    386, 436, 437, 338, 390, 391, 441, 442, 392, 344, 345, 395, 347,
    348, 398, 399, 451, 452, 402, 454, 455, 356, 457, 358, 359, 460,
    410, 462, 414, 464, 365, 417, 367, 368, 418, 419, 371, 372, 424,
    374, 424, 427, 428, 478, 430, 380, 432, 382, 432, 384, 485, 486,
    487, 388, 438, 390, 491, 441, 442, 494, 444, 445, 448, 398, 499,
    500, 452, 502, 503, 453, 505, 455, 456, 459, 460, 510, 411, 463,
    464, 514, 415, 467, 417, 467, 419, 520, 472, 422, 423, 424, 425,
    526, 427, 479, 429, 530, 431, 532, 533, 483, 535, 536, 488, 438,
    539, 489, 541, 493, 443, 493, 445, 446, 498, 548, 498, 499, 502,
    503, 504, 505, 555, 507, 508, 509, 510, 511, 561, 462, 463, 564,
    514, 517, 467, 568, 518, 570, 522, 472, 573, 574, 575, 527, 477,
    529, 479, 580, 581, 582, 583, 484, 585, 586, 536, 488, 538, 490,
    540, 592, 544, 494, 546, 547, 548, 598, 548, 600, 601, 551, 603,
    604, 556, 506, 507, 557, 509, 561, 611, 561, 562, 563, 564, 516,
    568, 569, 570, 571, 521, 571, 572, 575, 625, 575, 527, 528, 529,
    579, 580, 632, 582, 534, 584, 587, 588, 589, 539, 589, 590, 591,
    592, 593, 645, 546, 598, 599, 600, 650, 602, 652, 604, 605, 606,
    556, 557, 558, 610, 611, 661, 611, 663, 664, 616, 566, 567, 617,
    618, 621, 622, 572, 622, 625, 626, 627, 677, 578, 628, 631, 581,
    631, 632, 633, 685, 686, 587, 588, 640, 590, 640, 692, 693, 645,
    646, 696, 646, 698, 648, 700, 652, 653, 654, 604, 656, 706, 607,
    608, 609, 659, 662, 712, 664, 665, 615, 716, 666, 718, 619, 669,
    721, 671, 723, 673, 625, 675, 678, 628, 678, 630, 631, 683, 633,
    683, 686, 687, 637, 638, 688, 640, 641, 693, 694, 644, 694, 746,
    647, 697, 698, 750, 702, 652, 753, 654, 706, 707, 657, 658, 708,
    660, 710, 713, 663, 713, 714, 766, 667, 719, 669, 670, 722, 772,
    722, 725, 726, 776, 677, 678, 779, 729, 730, 731, 734, 735, 736,
    737, 687, 788, 738, 690, 740, 743, 744, 694, 744, 796, 748, 749,
    699, 749, 752, 753, 754, 704, 754, 806, 758, 708, 758, 761, 811,
    712, 764, 714, 815, 816, 766, 818, 719, 769, 821, 773, 823, 824,
    825, 777, 727, 777, 780, 830, 731, 781, 833, 783, 735, 787, 737,
    787, 788, 789, 741, 791, 792, 795, 845, 846, 847, 799, 749, 850,
    851, 803, 753, 803, 804, 856, 808, 758, 759, 760, 861, 811, 812,
    764, 814, 815, 867, 868, 820, 821, 822, 872, 822, 825, 875, 825,
    826, 778, 879, 880, 881, 782, 832, 884, 885, 786, 836, 888, 838,
    890, 840, 841, 844, 845, 895, 896, 846, 798, 899, 849, 901, 802,
    803, 853, 854, 855, 856, 857, 809, 861, 862, 863, 864, 865, 866,
    816, 866, 867, 868, 869, 870, 871, 872, 924, 876, 877, 927, 828,
    878, 879, 882, 883, 933, 883, 886, 836, 886, 938, 939, 889, 890,
    842, 892, 893, 945, 895, 896, 848, 898, 850, 902, 903, 904, 854,
    855, 856, 906, 907, 959, 911, 861, 913, 863, 913, 965, 866, 916,
    968, 869, 870, 871, 972, 873, 923, 924, 976, 977, 927, 928, 929,
    930, 931, 932, 884, 936, 937, 887, 937, 989, 939, 940, 892, 893,
    943, 995, 947, 897, 947, 950, 900, 952, 953, 903, 955, 956, 906,
    907, 957, 960, 961, 911, 961, 964, 914, 966, 916, 966, 967, 968,
    971, 972, 973, 974, 975, 925, 975, 976, 977, 929, 930, 931, 932,
    933, 983, 935, 987, 988, 938, 988, 940, 941, 993, 994, 944, 994,
    997, 947, 948, 998, 999]

    n = length(pi)
    @test all(pi .>= 1)
    @test all(pi .<= n)
    deg = node_degrees(pi)
    @test maximum(deg) == 3
    cidx = ChildrenIndex(pi)
    @test root_node(cidx) == 1
    @test length(cidx) == length(pi)
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
using Base: @propagate_inbounds

# AbstractArray interface

Base.IndexStyle(a::ExtensionArray) = IndexCartesian()
Base.parent(a::ExtensionArray) = a.parent
Base.eltype(a::ExtensionArray) = eltype(parent(a))
Base.axes(a::ExtensionArray) = axes(parent(a))
Base.size(a::ExtensionArray) = size(parent(a))
Base.length(a::ExtensionArray) = length(parent(a))
Base.iterate(a::ExtensionArray, state...) = iterate(parent(a), state...)
@propagate_inbounds Base.getindex(a::ExtensionArray, i...) =
    getindex(parent(a), i...)
@propagate_inbounds Base.setindex!(a::ExtensionArray, x, i...) =
    setindex!(parent(a), x, i...)

# ExtensionArray interface

extend(a::AbstractArray, extension::Extension) = ExtensionArray(a, extension)

extension(a::AbstractArray) = ExtensionNone()
extension(a::ExtensionArray) = a.extension

eltype_extension(a::ExtensionArray{<:Any,<:Any,<:Any,ExtensionConstant{S}}) where S = S

@propagate_inbounds getindex_extension(a::ExtensionArray, i) =
    getindex_extension(a, i, a.extension)
@propagate_inbounds setindex_extension!(a::ExtensionArray, x, i) =
    setindex_extension!(a, x, i, a.extension)


text:
# What is Meta? A place to store extra information. All of
# the extra information that I could find is here
mutable struct Meta
    rule
    line::Int64
    column::Int64
    end_line::Int64
    end_column::Int64
    start_pos::Int64
    end_pos::Int64
    empty::Bool
    Meta() = new()
end

# Note that a Python slotted tree is just a tree with
# restricted possible fields, like a Julia struct.
"""
The main tree class.

    Creates a new tree, and stores "data" and "children" in attributes of the same name.
    Trees can be hashed and compared.

    Parameters:
        data: The name of the rule or alias
        children: List of matched sub-rules and terminals
        meta: Line & Column numbers (if ``propagate_positions`` is enabled).
            meta attributes: line, column, start_pos, end_line, end_column, end_pos
"""
mutable struct Tree
    data
    children
    _meta::Meta
end

Tree(data,children) = Tree(data,children,Meta())
Tree(data,children,n::Nothing) = Tree(data,children)

meta(t::Tree) = t._meta

_pretty_label(t::Tree) = t.data

pretty(t::Tree,level,indent_str::String) = begin
    if length(t.children) == 1 && !(typeof(t.children[1]) == Tree)
        return [repeat(indent_str,level), _pretty_label(t), "\t", "$(t.children[1])", "\n"]
    end

    l = [ repeat(indent_str,level), _pretty_label(t), "\n" ]
    for n in t.children
        if typeof(n) == Tree
            append!(l, pretty(n,level+1, indent_str))
        else
            append!(l, [ repeat(indent_str,level+1), "$n", "\n" ])
        end
    end
    return l
end

pretty(t::Tree; indent_str = "  ") = join(pretty(t,0, indent_str),"")

#Base.show(io::IO,t::Tree) = println(io,pretty(t))
Base.show(io::IO,t::Tree) = println(io,"Tree($(t.data), $(t.children))")
"""
Expand children at the given indices. Rewritten for Julia as we
can't expand inline.
"""
expand_kids_by_index!(t::Tree, indices) = begin
    #println("Expanding $t and $indices")
    new_children = []
    for (i,c) in enumerate(t.children)
        if i in indices
            append!(new_children,c.children)
        else
            push!(new_children,c)
        end
    end
    t.children = new_children
    #println("Returning: $t")
end

Base.:(==)(t1::Tree,t2::Tree) = begin
    return t1.data == t2.data && t1.children == t2.children
end

# When comparing a tree with something else, we compare only
# on data (as a string) and children. This is mainly for testing.
#==Base.:(==)(t1::Tree,t2::Array) = begin
    if length(t2) != 2
        println("Not equal: length != 2")
        return false
    end
    if t2[1] != t1.data
        println("Not equal: $(t2[1]) != $(t1.data)")
        return false
    end
    if length(t2[2]) != length(t1.children)
        println("Not equal: child length differs")
        return false
    end
    for (f,s) in zip(t1.children,t2[2])
        if f != s
            println("Not equal: $f != $s")
            return false
        end
    end
    return true
end
==#

Base.hash(t1::Tree) = begin
    return hash(t1.data,hash(t1.children))
end

find_pred(t::Tree,pred) = begin
    Iterators.filter(pred, t)
end

find_data(t::Tree,data) = begin
    find_pred(t, x -> x.data == data)
end

# Test only non-tree values for `pred`
scan_values(t::Tree,pred) = begin
    result = Any[]
    for c in t.children
        if typeof(c) == Tree
            append!(result,scan_values(c,pred))
        else
            if pred(c)
                push!(result,c)
            end
        end
    end
    return result
end

#==
scan_values(t::Tree,pred) = Channel() do val_chan
    begin
        for c in t.children
            if typeof(c) == Tree
                for t in scan_values(c,pred)
                    put!(val_chan,t)
                end
            else
                if pred(c)
                    put!(val_chan,c)
                end
            end
        end
    end
end
==#

Base.iterate(t::Tree) = begin
    visited = Set()
    q = [t]
    l = []
    while !isempty(q)
        subtree = pop!(q)
        push!(l,subtree)
        if objectid(subtree) in visited
            continue
        end
        push!(visited,objectid(subtree))
        append!(q,[c for c in subtree.children if typeof(c) == Tree])
    end
    seen = Set()
    idx = length(l)
    if idx == 0 return nothing end
    push!(seen,objectid(l[idx]))
    return l[idx],(seen,l,idx)
end

# `s` is a tuple `(seen,l,idx)` where `idx` is the most
# recently returned index.
Base.iterate(t::Tree,s) = begin
    seen,l,idx = s
    idx = idx -1
    while idx > 0 && objectid(l[idx]) in seen
        idx = idx - 1
    end
    if idx == 0 return nothing end
    push!(seen,objectid(l[idx]))
    return l[idx],(seen,l,idx)
end

tree_set(t::Tree,data,children) = begin
    t.data = data
    t.children = children
end


text:
#=  After Azevedo and Leshno (2016, §IV.A), an example of how in some situations
	schools may be disincentivized to improve their quality.   =#

using DeferredAcceptance

function argsort(vec::AbstractArray{<:Real, 1})::AbstractArray{<:Real, 1}
    return invperm(sortperm(vec))
end

# Might add these functions to the module later.
function quality(assn, capacities, scores)
	(m, ) = size(capacities)
	res = zeros(m)
	for (s, c) in enumerate(assn)
		if c ≤ m
			res[c] += scores[s, c]
		end
	end
	res ./= capacities
	return res
end

function cutoffs(assn, capacities, scores)
	(m, ) = size(capacities)
	res = zeros(m)
	for c in 1:m
		res[c] = minimum(scores[:, c][assn .== c])
	end
	return res
end

students = [2 1 1 1 1;
			1 2 2 2 2]

# School 2 improves, winning over one applicant
students_delta = [2 2 1 1 1;
				  1 1 2 2 2]

# For illustration purposes
scores = [ 0.95  0.96  0.89  0.83  0.61;
 		   0.97  0.84  0.82  0.9   0.68 ]

schools = mapslices(argsort, 1 .- scores, dims=2)

@assert schools == [2 1 3 4 5;
		            1 3 4 2 5]

capacities = [2, 3]

out_normal = DA(students, schools', capacities)
out_delta = DA(students_delta, schools', capacities)

println("Normal assn:  ", out_normal[1])
println("    Quality:  ", quality(out_normal[1], capacities, scores'))
println("    Cutoffs:  ", cutoffs(out_normal[1], capacities, scores'))
println("Delta:        ", out_delta[1])
println("    Quality:  ", quality(out_delta[1], capacities, scores'))
println("    Cutoffs:  ", cutoffs(out_delta[1], capacities, scores'))
# Normal assn:  [2, 1, 1, 2, 2]
#     Quality:  [0.925, 0.8500000000000001]
#     Cutoffs:  [0.89, 0.68]
# Delta:        [2, 2, 1, 1, 2]
#     Quality:  [0.86, 0.8300000000000001]
#     Cutoffs:  [0.83, 0.68]

#=  School 2 ends up losing a great student (4), because persuading student 2
	away from school 1 meant that school 1 lowered its cutoff and 4 got in. =#

# Shows that the stable match is unique.
println("Normal-rev:   ",DA(students, schools', capacities, rev=true)[1])
println("Delta-rev:    ",DA(students_delta, schools', capacities, rev=true)[1])
# Normal-rev:   [2, 1, 1, 2, 2]
# Delta-rev:    [2, 2, 1, 1, 2]

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
function reflection_coefficient(ω::T, m_wave::MatchWave{T}, medium::Medium{T}, specie::Specie{T};
            kws...) where T<:AbstractFloat

    R_eff = reflection_coefficient(ω, m_wave.effective_waves, medium, [specie]; x = m_wave.x_match[end], kws...)
    R_discrete = reflection_coefficient(ω, m_wave.average_wave, medium, specie; kws...)

    return  R_eff + R_discrete
end


text:
# Autogenerated wrapper script for UDUNITS_jll for x86_64-linux-gnu
export libudunits2, udunits2

using Expat_jll
JLLWrappers.@generate_wrapper_header("UDUNITS")
JLLWrappers.@declare_library_product(libudunits2, "libudunits2.so")
JLLWrappers.@declare_executable_product(udunits2)
function __init__()
    JLLWrappers.@generate_init_header(Expat_jll)
    JLLWrappers.@init_library_product(
        libudunits2,
        "lib/libudunits2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        udunits2,
        "bin/udunits2",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using SparseArrays
using ViewerGL
GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#Test spatial_arrangement on octahemioctahedron DEBUG
octahemioctahedron = Lar.octahemioctahedron()
V,EV,FV,CV = Lar.struct2lar(octahemioctahedron)
GL.VIEW([ GL.GLGrid(V,FV, GL.Point4d(1,1,1,0.1))]);

cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');#cobordo 1

V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(
W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);


text:
# syntax: proto3
using Compat
using ProtoBuf
import ProtoBuf.meta

struct __enum_ConnectionRequest_Type <: ProtoEnum
    RPC::Int32
    STREAM::Int32
    __enum_ConnectionRequest_Type() = new(0,1)
end #struct __enum_ConnectionRequest_Type
const ConnectionRequest_Type = __enum_ConnectionRequest_Type()

mutable struct ConnectionRequest <: ProtoType
    _type::Int32
    client_name::AbstractString
    client_identifier::Array{UInt8,1}
    ConnectionRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ConnectionRequest

struct __enum_ConnectionResponse_Status <: ProtoEnum
    OK::Int32
    MALFORMED_MESSAGE::Int32
    TIMEOUT::Int32
    WRONG_TYPE::Int32
    __enum_ConnectionResponse_Status() = new(0,1,2,3)
end #struct __enum_ConnectionResponse_Status
const ConnectionResponse_Status = __enum_ConnectionResponse_Status()

mutable struct ConnectionResponse <: ProtoType
    status::Int32
    message::AbstractString
    client_identifier::Array{UInt8,1}
    ConnectionResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ConnectionResponse

mutable struct Argument <: ProtoType
    position::UInt32
    value::Array{UInt8,1}
    Argument(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Argument

mutable struct ProcedureCall <: ProtoType
    service::AbstractString
    procedure::AbstractString
    service_id::UInt32
    procedure_id::UInt32
    arguments::Base.Vector{Argument}
    ProcedureCall(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ProcedureCall
const __fnum_ProcedureCall = Int[1,2,4,5,3]
meta(t::Type{ProcedureCall}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ProcedureCall, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct Request <: ProtoType
    calls::Base.Vector{ProcedureCall}
    Request(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Request

mutable struct Error <: ProtoType
    service::AbstractString
    name::AbstractString
    description::AbstractString
    stack_trace::AbstractString
    Error(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Error

mutable struct ProcedureResult <: ProtoType
    error::Error
    value::Array{UInt8,1}
    ProcedureResult(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ProcedureResult

mutable struct Response <: ProtoType
    error::Error
    results::Base.Vector{ProcedureResult}
    Response(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Response

mutable struct StreamResult <: ProtoType
    id::UInt64
    result::ProcedureResult
    StreamResult(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StreamResult

mutable struct StreamUpdate <: ProtoType
    results::Base.Vector{StreamResult}
    StreamUpdate(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StreamUpdate

mutable struct Class <: ProtoType
    name::AbstractString
    documentation::AbstractString
    Class(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Class

mutable struct EnumerationValue <: ProtoType
    name::AbstractString
    value::Int32
    documentation::AbstractString
    EnumerationValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EnumerationValue

mutable struct Enumeration <: ProtoType
    name::AbstractString
    values::Base.Vector{EnumerationValue}
    documentation::AbstractString
    Enumeration(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Enumeration

mutable struct Exception <: ProtoType
    name::AbstractString
    documentation::AbstractString
    Exception(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Exception

struct __enum_Type_TypeCode <: ProtoEnum
    NONE::Int32
    DOUBLE::Int32
    FLOAT::Int32
    SINT32::Int32
    SINT64::Int32
    UINT32::Int32
    UINT64::Int32
    BOOL::Int32
    STRING::Int32
    BYTES::Int32
    CLASS::Int32
    ENUMERATION::Int32
    EVENT::Int32
    PROCEDURE_CALL::Int32
    STREAM::Int32
    STATUS::Int32
    SERVICES::Int32
    TUPLE::Int32
    LIST::Int32
    SET::Int32
    DICTIONARY::Int32
    __enum_Type_TypeCode() = new(0,1,2,3,4,5,6,7,8,9,100,101,200,201,202,203,204,300,301,302,303)
end #struct __enum_Type_TypeCode
const Type_TypeCode = __enum_Type_TypeCode()

mutable struct _Type <: ProtoType
    code::Int32
    service::AbstractString
    name::AbstractString
    types::Base.Vector{_Type}
    _Type(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct _Type

mutable struct Parameter <: ProtoType
    name::AbstractString
    _type::_Type
    default_value::Array{UInt8,1}
    Parameter(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Parameter

mutable struct Procedure <: ProtoType
    name::AbstractString
    parameters::Base.Vector{Parameter}
    return_type::_Type
    return_is_nullable::Bool
    documentation::AbstractString
    Procedure(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Procedure

mutable struct Service <: ProtoType
    name::AbstractString
    procedures::Base.Vector{Procedure}
    classes::Base.Vector{Class}
    enumerations::Base.Vector{Enumeration}
    exceptions::Base.Vector{Exception}
    documentation::AbstractString
    Service(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Service

mutable struct Services <: ProtoType
    services::Base.Vector{Service}
    Services(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Services

mutable struct Tuple_ <: ProtoType
    items::Base.Vector{Array{UInt8,1}}
    Tuple_(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Tuple

mutable struct List <: ProtoType
    items::Base.Vector{Array{UInt8,1}}
    List(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct List

mutable struct Set <: ProtoType
    items::Base.Vector{Array{UInt8,1}}
    Set(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Set

mutable struct DictionaryEntry <: ProtoType
    key::Array{UInt8,1}
    value::Array{UInt8,1}
    DictionaryEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DictionaryEntry

mutable struct Dictionary <: ProtoType
    entries::Base.Vector{DictionaryEntry}
    Dictionary(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Dictionary

mutable struct Stream <: ProtoType
    id::UInt64
    Stream(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Stream

mutable struct Event <: ProtoType
    stream::Stream
    Event(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Event

mutable struct Status <: ProtoType
    version::AbstractString
    bytes_read::UInt64
    bytes_written::UInt64
    bytes_read_rate::Float32
    bytes_written_rate::Float32
    rpcs_executed::UInt64
    rpc_rate::Float32
    one_rpc_per_update::Bool
    max_time_per_update::UInt32
    adaptive_rate_control::Bool
    blocking_recv::Bool
    recv_timeout::UInt32
    time_per_rpc_update::Float32
    poll_time_per_rpc_update::Float32
    exec_time_per_rpc_update::Float32
    stream_rpcs::UInt32
    stream_rpcs_executed::UInt64
    stream_rpc_rate::Float32
    time_per_stream_update::Float32
    Status(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Status

mutable struct MultiplexedRequest <: ProtoType
    connection_request::ConnectionRequest
    request::Request
    MultiplexedRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MultiplexedRequest

mutable struct MultiplexedResponse <: ProtoType
    response::Response
    stream_update::StreamUpdate
    MultiplexedResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MultiplexedResponse

export ConnectionRequest_Type, ConnectionRequest, ConnectionResponse_Status, ConnectionResponse, Request, ProcedureCall, Argument, Response, ProcedureResult, Error, StreamUpdate, StreamResult, Services, Service, Procedure, Parameter, Class, Enumeration, EnumerationValue, Exception, Type_TypeCode, _Type, Tuple, List, Set, Dictionary, DictionaryEntry, Stream, Event, Status, MultiplexedRequest, MultiplexedResponse


text:
@testset "cartesian" begin

# test of cartesian
sc = SparkContext(master="local")

nums1 = parallelize(sc, 1:3)
nums2 = parallelize(sc, 1:2)
nums3 = map(nums2, it -> it + 10)
rdd = cartesian(nums1, nums3)
rdd2 = map(rdd, it-> it[1] + it[2])
@test reduce(rdd2, +) == 81

close(sc)

end


text:
using ParameterizedFunctions, MATLABDiffEq, OrdinaryDiffEq, ODEInterface,
      ODEInterfaceDiffEq, Plots, Sundials, SciPyDiffEq, deSolveDiffEq
using DiffEqDevTools
using LinearAlgebra

f = @ode_def_bare LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
p = [1.5,1,3,1]
tspan = (0.0,10.0)
u0 = [1.0,1.0]
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern7())
          Dict(:alg=>MATLABDiffEq.ode45())
          Dict(:alg=>MATLABDiffEq.ode113())
          Dict(:alg=>SciPyDiffEq.RK45())
          Dict(:alg=>SciPyDiffEq.LSODA())
          Dict(:alg=>SciPyDiffEq.odeint())
          Dict(:alg=>deSolveDiffEq.lsoda())
          Dict(:alg=>deSolveDiffEq.ode45())
          Dict(:alg=>CVODE_Adams())
  ]

names = [
  "Julia: DP5"
  "Hairer: dopri5"
  "Julia: Tsit5"
  "Julia: Vern7"
  "MATLAB: ode45"
  "MATLAB: ode113"
  "SciPy: RK45"
  "SciPy: LSODA"
  "SciPy: odeint"
  "deSolve: lsoda"
  "deSolve: ode45"
  "Sundials: Adams"
  ]

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10)
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names = names,
                      appxsol=test_sol,dense=false,
                      save_everystep=false,numruns=100,maxiters=10000000,
                      timeseries_errors=false,verbose=false)
plot(wp,title="Non-stiff 1: Lotka-Volterra")

f = @ode_def_bare RigidBodyBench begin
  dy1  = -2*y2*y3
  dy2  = 1.25*y1*y3
  dy3  = -0.5*y1*y2 + 0.25*sin(t)^2
end
prob = ODEProblem(f,[1.0;0.0;0.9],(0.0,100.0))
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern7())
          Dict(:alg=>MATLABDiffEq.ode45())
          Dict(:alg=>MATLABDiffEq.ode113())
          Dict(:alg=>SciPyDiffEq.RK45())
          Dict(:alg=>SciPyDiffEq.LSODA())
          Dict(:alg=>SciPyDiffEq.odeint())
          Dict(:alg=>deSolveDiffEq.lsoda())
          Dict(:alg=>deSolveDiffEq.ode45())
          Dict(:alg=>CVODE_Adams())
  ]

names = [
  "Julia: DP5"
  "Hairer: dopri5"
  "Julia: Tsit5"
  "Julia: Vern7"
  "MATLAB: ode45"
  "MATLAB: ode113"
  "SciPy: RK45"
  "SciPy: LSODA"
  "SciPy: odeint"
  "deSolve: lsoda"
  "deSolve: ode45"
  "Sundials: Adams"
  ]

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10)
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names = names,
                      appxsol=test_sol,dense=false,
                      save_everystep=false,numruns=100,maxiters=10000000,
                      timeseries_errors=false,verbose=false)
plot(wp,title="Non-stiff 2: Rigid-Body")

rober = @ode_def begin
  dy₁ = -k₁*y₁+k₃*y₂*y₃
  dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  dy₃ =  k₂*y₂^2
end k₁ k₂ k₃
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])
sol = solve(prob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (7:8)
reltols = 1.0 ./ 10.0 .^ (3:4);

setups = [Dict(:alg=>Rosenbrock23())
          Dict(:alg=>TRBDF2())
          Dict(:alg=>RadauIIA5())
          Dict(:alg=>rodas())
          Dict(:alg=>radau())
          Dict(:alg=>MATLABDiffEq.ode23s())
          Dict(:alg=>MATLABDiffEq.ode15s())
          Dict(:alg=>SciPyDiffEq.LSODA())
          Dict(:alg=>SciPyDiffEq.BDF())
          Dict(:alg=>SciPyDiffEq.odeint())
          Dict(:alg=>deSolveDiffEq.lsoda())
          Dict(:alg=>CVODE_BDF())
          ]

names = [
  "Julia: Rosenbrock23"
  "Julia: TRBDF2"
  "Julia: radau"
  "Hairer: rodas"
  "Hairer: radau"
  "MATLAB: ode23s"
  "MATLAB: ode15s"
  "SciPy: LSODA"
  "SciPy: BDF"
  "SciPy: odeint"
  "deSolve: lsoda"
  "Sundials: CVODE"
  ]

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names = names,print_names = true,
                      dense=false,verbose = false,
                      save_everystep=false,appxsol=test_sol,
                      maxiters=Int(1e5))
plot(wp,title="Stiff 1: ROBER", legend=:topleft)

f = @ode_def Hires begin
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
           0.43*y6 + 0.69*y7
  dy7 = 280.0*y6*y8 - 1.81*y7
  dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(f,u0,(0.0,321.8122))

sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);

setups = [Dict(:alg=>Rosenbrock23())
          Dict(:alg=>TRBDF2())
          Dict(:alg=>RadauIIA5())
          Dict(:alg=>rodas())
          Dict(:alg=>radau())
          Dict(:alg=>MATLABDiffEq.ode23s())
          Dict(:alg=>MATLABDiffEq.ode15s())
          Dict(:alg=>SciPyDiffEq.LSODA())
          Dict(:alg=>SciPyDiffEq.BDF())
          Dict(:alg=>SciPyDiffEq.odeint())
          Dict(:alg=>deSolveDiffEq.lsoda())
          Dict(:alg=>CVODE_BDF())
          ]

names = [
  "Julia: Rosenbrock23"
  "Julia: TRBDF2"
  "Julia: radau"
  "Hairer: rodas"
  "Hairer: radau"
  "MATLAB: ode23s"
  "MATLAB: ode15s"
  "SciPy: LSODA"
  "SciPy: BDF"
  "SciPy: odeint"
  "deSolve: lsoda"
  "Sundials: CVODE"
  ]

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names = names,print_names = true,
                      save_everystep=false,appxsol=test_sol,
                      maxiters=Int(1e5),numruns=100)
plot(wp,title="Stiff 2: Hires",legend=:topleft)


text:
function consumer_amqp_config!(chan, name)
    declare_queue!(chan, name)
end


text:
using ITensors,
      Test
import ITensors.SmallString

@testset "QN" begin

  @testset "QNVal Basics" begin
    qv = QNVal()
    @test !isactive(qv)

    qv = QNVal("Sz",0)
    @test name(qv) == SmallString("Sz")
    @test val(qv) == 0
    @test modulus(qv) == 1
    @test isactive(qv)

    qv = QNVal("A",1,2)
    @test name(qv) == SmallString("A")
    @test val(qv) == 1
    @test modulus(qv) == 2
    @test !isfermionic(qv)

    qv = QNVal("Nf",1,-1)
    @test name(qv) == SmallString("Nf")
    @test val(qv) == 1
    @test modulus(qv) == -1
    @test isfermionic(qv)
  end

  @testset "QN Basics" begin
    q = QN()
    @test length(sprint(show,q)) > 1

    q = QN(("Sz",1))
    @test length(sprint(show,q)) > 1
    @test isactive(q[1])
    @test val(q,"Sz") == 1

    q = QN("Sz",1)
    @test length(sprint(show,q)) > 1
    @test isactive(q[1])
    @test val(q,"Sz") == 1

    q = QN("P",1,2)
    @test length(sprint(show,q)) > 1
    @test isactive(q[1])
    @test val(q,"P") == 1
    @test modulus(q,"P") == 2

    q = QN(("A",1),("B",2))
    @test isactive(q[1])
    @test isactive(q[2])
    @test val(q,"A") == 1
    @test val(q,"B") == 2
    @test modulus(q,"A") == 1
    @test modulus(q,"B") == 1

    q = QN(("B",2),("A",1))
    @test val(q,"A") == 1
    @test val(q,"B") == 2
  end

  @testset "Comparison" begin
    @test QN() == QN()
    @test QN("A",1) == QN("A",1)
    @test QN(("A",1),("B",3)) == QN(("A",1),("B",3))
    @test QN(("A",1),("B",3)) == QN(("B",3),("A",1))

    # Zero value and missing sector treated the same:
    @test QN(("A",0),("B",3)) == QN("B",3)
    @test QN(("B",3),("A",0)) == QN("B",3)
  end

  @testset "Arithmetic" begin
    @test QN("Sz",1) + QN() == QN("Sz",1)
    @test QN("Sz",1) + QN("Sz",2) == QN("Sz",3)
    @test QN("Sz",1) + QN("Sz",-2) == QN("Sz",-1)

    @test QN(("A",1),("Sz",0)) + QN(("A",0),("Sz",1)) == QN(("A",1),("Sz",1))

    @test QN("P",0,2) + QN("P",1,2) == QN("P",1,2)
    @test QN("P",1,2) + QN("P",1,2) == QN("P",0,2)
  end

  @testset "Ordering" begin
    z = QN()
    qa = QN(("Sz",1),("Nf",1))
    qb = QN(("Sz",0),("Nf",2))
    qc = QN(("Sz",1),("Nf",2))
    qd = QN(("Sz",1),("Nf",2))
    qe = QN(("Sz",-1),("Nf",-2))

    @test !(z < z)
    @test !(qa < z)
    @test (z < qa)
    @test (z < qb)
    @test !(qb < z)
    @test (z < qc)
    @test !(qc < z)
    @test (z < qd)
    @test !(qd < z)
    @test !(z < qe)
    @test (qe < z)

    @test (qa > qb)
    @test !(qb > qa)
    @test !(qb == qa)
    @test (qb < qc)
    @test !(qc < qb)
    @test !(qc == qb)
    @test (qc == qd)
    @test !(qc < qd)
    @test !(qd < qc)
  end

end


text:
# This file is a part of StruckVMEDevices.jl, licensed under the MIT License (MIT).

function read_data(input::IO; nbuffers = typemax(Int))
    daqevtno = Vector{Int32}()
    bufferno = Vector{Int32}()
    channel = Vector{Int32}()
    timestamp = Vector{Int64}()
    energy = Vector{Int32}()
    trigmax = Vector{Int32}()
    samples = VectorOfVectors{UInt32}()

    tmpbuffer = Vector{Int32}()
    tmpevtdata = Vector{UInt8}()

    evtno = 0
    nbufread = 0

    while !eof(input) && nbufread < nbuffers
        bufinfo = read(input, BankChannelHeaderInfo)
        nbufread += 1
        bufno = bufinfo.bufferno

        resize!(tmpevtdata, sizeof(UInt32) * bufinfo.nevents * bufinfo.nwords_per_event)
        read!(input, tmpevtdata)
        evtdatabuf = IOBuffer(tmpevtdata)

        for i in 1:bufinfo.nevents
            evtno += 1
            evt = read(evtdatabuf, RawChEvent, bufinfo.nmawvalues, bufinfo.firmware_type, tmpbuffer)

            push!(daqevtno, evtno)
            push!(channel, evt.chid + 1)
            push!(bufferno, bufno)
            push!(timestamp, evt.timestamp)
            push!(energy, isnothing(evt.energy) ? Int32(0) : evt.energy.maximum)
            push!(trigmax, isnothing(evt.trig_maw) ? Int32(0) : evt.trig_maw.maximum)
            push!(samples, evt.samples)
        end
        # info("Read buffer $bufno, channel $chno with $(length(buffer.events)) events")
    end
    # info("Read $nbufread buffers * channels")

    Table(
        daqevtno = daqevtno,
        bufferno = bufferno,
        channel = channel,
        timestamp = timestamp,
        energy = energy,
        trigmax = trigmax,
        samples = samples,
    )
end


text:
@testset "spr" begin
    binary_tree = ParseNewick("((raccoon:19.19959,bear:6.80041):0.84600,((sea_lion:11.99700,seal:12.00300):7.52973,((monkey:100.85930,cat:47.14069):20.59201,weasel:18.87953):2.09460):3.87382);")

    error_tree = ParseNewick("(raccoon:19.19959):0.84600;")
    error_tree_not_binary = ParseNewick("((raccoon:19.19959,bear:6.80041):0.84600,((sea_lion:11.99700,seal:12.00300,lizard:5.03):7.52973,((monkey:100.85930,cat:47.14069):20.59201,weasel:18.87953):2.09460):3.87382);")
    tree_binary_two = ParseNewick("((raccoon:19.19959,bear:6.80041)50:0.84600,((sea_lion:11.99700,seal:12.00300)100:7.52973,((monkey:100.85930,cat:47.14069)80:20.59201,weasel:18.87953)75:2.09460)50:3.87382,dog:25.46154);")

    @test false == check_binary(error_tree)
    @test false == check_binary(error_tree_not_binary)
    @test_throws ErrorException MCPhyloTree.SPR(error_tree)
    @test true == check_binary(tree_binary_two)
    spr_binary = MCPhyloTree.SPR(binary_tree)

    @test length(Set([n.num for n in post_order(spr_binary)])) == length(Set([n.num for n in post_order(binary_tree)]))

    @test round(tree_length(binary_tree);digits=3) == round(tree_length(spr_binary);digits=3)

    @test length(post_order(spr_binary)) == length(post_order(binary_tree))

    spr_binary_two = MCPhyloTree.SPR(tree_binary_two)

    @test length(Set([n.num for n in post_order(spr_binary_two)])) == length(Set([n.num for n in post_order(tree_binary_two)]))

    @test round(tree_length(tree_binary_two);digits=3) == round(tree_length(spr_binary_two);digits=3)

    @test length(post_order(spr_binary_two)) == length(post_order(tree_binary_two))

    spr_binary_risky = MCPhyloTree.risky_SPR(binary_tree)

    @test length(Set([n.num for n in post_order(spr_binary_risky)])) == length(Set([n.num for n in post_order(binary_tree)]))
    @test round(tree_length(spr_binary_risky);digits=3) == round(tree_length(spr_binary);digits=3)
    @test length(post_order(spr_binary_risky)) == length(post_order(binary_tree))
end


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

@testset "Reproject a Geometry" begin
    println("Method 1")
    source = AG.unsafe_importEPSG(2927); target = AG.unsafe_importEPSG(4326)
        transform = AG.unsafe_createcoordtrans(source, target)
            point = AG.unsafe_fromWKT("POINT (1120351.57 741921.42)")
                println("Before: $(AG.toWKT(point))")
                AG.transform!(point, transform)
                println("After: $(AG.toWKT(point))")
    AG.destroy(point); AG.destroy(transform);
    AG.destroy(target); AG.destroy(source)

    println("Method 2")
    AG.importEPSG(2927) do source; AG.importEPSG(4326) do target
        AG.createcoordtrans(source, target) do transform
            AG.fromWKT("POINT (1120351.57 741921.42)") do point
                println("Before: $(AG.toWKT(point))")
                AG.transform!(point, transform)
                println("After: $(AG.toWKT(point))")
    end end end end
end

@testset "Get Projection" begin
    AG.registerdrivers() do
        AG.read("ospy/data1/sites.shp") do dataset
            layer = AG.getlayer(dataset, 0)
            spatialref = AG.getspatialref(layer)
            println(AG.toWKT(spatialref))
            println(AG.toWKT(spatialref, false))
            println(AG.toWKT(spatialref, true))
            println(AG.toPROJ4(spatialref))
            println(AG.toXML(spatialref))
            println(AG.toMICoordSys(spatialref))
            AG.nextfeature(layer) do feature
                feature |>
                    AG.getgeom |>
                    AG.getspatialref |>
                    AG.toWKT |>
                    println
            end
        end
    end

    AG.importEPSG(26912) do spatialref
        println("before (Proj4): $(AG.toPROJ4(spatialref))")
        println("before (WKT): $(AG.toWKT(spatialref))")
        AG.morphtoESRI!(spatialref)
        println("after (Proj4): $(AG.toPROJ4(spatialref))")
        println("after: $(AG.toWKT(spatialref))")
    end
end


text:
"""
    AbstractLinearSolver

Abstract type for solvers of interval linear systems.
"""
abstract type AbstractLinearSolver end

"""
    AbstractDirectSolver <: AbstractLinearSolver

Abstract type for direct solvers of interval linear systems, such as Gaussian elimination
and Hansen-Bliek-Rohn.
"""
abstract type AbstractDirectSolver <: AbstractLinearSolver end

"""
    AbstractIterativeSolver <: AbstractLinearSolver

Abstract type for iterative solvers of interval linear systems, such as Jacobi or
Gauss-Seidel.
"""
abstract type AbstractIterativeSolver <: AbstractLinearSolver end

"""
    HansenBliekRohn <: AbstractDirectSolver

Type for the `HansenBliekRohn` solver of the square interval linear system ``Ax=b``.
For more details see section 5.6.2 of [[HOR19]](@ref)

### Notes

- Hansen-Bliek-Rohn works with H-matrices without precondition and with strongly regular
  matrices using [`InverseMidpoint`](@ref) precondition
- If the midpoint of ``A`` is a diagonal matrix, then the algorithm returns the exact hull.
- An object of type Hansen-Bliek-Rohn is a callable function with method

        (hbr::HansenBliekRohn)(A::AbstractMatrix{T},
                               b::AbstractVector{T}) where {T<:Interval}

### Examples

```jldoctest
julia> A = [2..4 -1..1;-1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> b = [-2..2, -1..1]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-1, 1]

julia> hbr = HansenBliekRohn()
HansenBliekRohn linear solver

julia> hbr(A, b)
2-element Vector{Interval{Float64}}:
 [-1.66667, 1.66667]
 [-1.33334, 1.33334]
```
"""
struct HansenBliekRohn <: AbstractDirectSolver end

function (hbr::HansenBliekRohn)(A::AbstractMatrix{T},
                                b::AbstractVector{T}) where {T<:Interval}
    n = length(b)
    compA = comparison_matrix(A)
    compA_inv, cert = epsilon_inflation(compA, Diagonal(ones(n)))

    cert || @warn "Could not find a verified enclosure of ⟨A⟩⁻¹"

    all(sum(compA_inv; dims=2) .> 0) || throw(ArgumentError("applying Hanben-Bliek-Rohn to a non-H-matrix."))
    u = compA_inv * mag.(b)
    d = diag(compA_inv)

    _α = sup.(diag(compA) .- 1 ./ d)
    α = Interval.(-_α, _α)

    _β = @. sup(u/d - mag(b))
    β = Interval.(-_β, _β)

    return (b .+ β) ./ (diag(A) .+ α)

end

"""
    GaussianElimination <: AbstractDirectSolver

Type for the Gaussian elimination solver of the square interval linear system ``Ax=b``.
For more details see section 5.6.1 of [[HOR19]](@ref)

### Notes

- An object of type `GaussianElimination` is a callable function with method

        (ge::GaussianElimination)(A::AbstractMatrix{T},
                                  b::AbstractVector{T}) where {T<:Interval}

### Examples

```jldoctest
julia> A = [2..4 -1..1;-1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> b = [-2..2, -1..1]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-1, 1]

julia> ge = GaussianElimination()
GaussianElimination linear solver

julia> ge(A, b)
2-element Vector{Interval{Float64}}:
 [-1.66667, 1.66667]
 [-1.33334, 1.33334]
```
"""
struct GaussianElimination <: AbstractDirectSolver end

function (ge::GaussianElimination)(A::AbstractMatrix{T},
                                   b::AbstractVector{T}) where {T<:Interval}
    n = length(b)
    Abrref = rref([A b])

    # backsubstitution
    x = similar(b)
    x[end] = Abrref[n, n+1]/Abrref[n, n]
    @inbounds for i = n-1:-1:1
        x[i] = (Abrref[i, n+1] - sum(Abrref[i, j]*x[j] for j in i+1:n))/Abrref[i, i]
    end
    return x
end

## JACOBI
"""
    Jacobi <: AbstractIterativeSolver

Type for the Jacobi solver of the interval linear system ``Ax=b``.
For details see Section 5.7.4 of [[HOR19]](@ref)

### Fields

- `max_iterations` -- maximum number of iterations (default 20)
- `atol`           -- absolute tolerance (default 0), if at some point ``|xₖ - xₖ₊₁| < atol``
                      (elementwise), then stop and return ``xₖ₊₁``.
                      If `atol=0`, then `min(diam(A))*1e-5` is used.

### Notes

- An object of type `Jacobi` is a function with method

        (jac::Jacobi)(A::AbstractMatrix{T},
                      b::AbstractVector{T},
                      [x]::AbstractVector{T}=enclose(A, b)) where {T<:Interval}

    #### Input
    - `A`   -- N×N interval matrix
    - `b`   -- interval vector of length N
    - `x`   -- (optional) initial enclosure for the solution of ``Ax = b``. If not given,
               it is automatically computed using [`enclose`](@ref enclose)

### Examples

```jldoctest
julia> A = [2..4 -1..1;-1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> b = [-2..2, -1..1]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-1, 1]

julia> jac = Jacobi()
Jacobi linear solver
max_iterations = 20
atol = 0.0

julia> jac(A, b)
2-element Vector{Interval{Float64}}:
 [-1.66668, 1.66668]
 [-1.33335, 1.33335]
```
"""
struct Jacobi <: AbstractIterativeSolver
    max_iterations::Int
    atol::Float64
end

Jacobi() = Jacobi(20, 0.0)

function (jac::Jacobi)(A::AbstractMatrix{T},
                       b::AbstractVector{T},
                       x::AbstractVector{T}=enclose(A, b)) where {T<:Interval}

    n = length(b)
    atol = iszero(jac.atol) ? minimum(diam.(A))*1e-5 : jac.atol

    for _ in 1:jac.max_iterations
        xold = copy(x)
        @inbounds @simd for i in 1:n
            x[i] = b[i]
            for j in 1:n
                (i == j) || (x[i] -= A[i, j] * xold[j])
            end
            x[i] = (x[i]/A[i, i]) ∩ xold[i]
        end
        all(interval_isapprox.(x, xold; atol=atol)) && break
    end
    return x
end

## GAUSS SEIDEL
"""
    GaussSeidel <: AbstractIterativeSolver

Type for the Gauss-Seidel solver of the interval linear system ``Ax=b``.
For details see Section 5.7.4 of [[HOR19]](@ref)

### Fields

- `max_iterations` -- maximum number of iterations (default 20)
- `atol`           -- absolute tolerance (default 0), if at some point ``|xₖ - xₖ₊₁| < atol``
                      (elementwise), then stop and return ``xₖ₊₁``.
                      If `atol=0`, then `min(diam(A))*1e-5` is used.

### Notes

- An object of type `GaussSeidel` is a function with method

        (gs::GaussSeidel)(A::AbstractMatrix{T},
                          b::AbstractVector{T},
                          [x]::AbstractVector{T}=enclose(A, b)) where {T<:Interval}

    #### Input
    - `A`   -- N×N interval matrix
    - `b`   -- interval vector of length N
    - `x`   -- (optional) initial enclosure for the solution of ``Ax = b``. If not given,
               it is automatically computed using [`enclose`](@ref enclose)

### Examples

```jldoctest
julia> A = [2..4 -1..1;-1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> b = [-2..2, -1..1]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-1, 1]

julia> gs = GaussSeidel()
GaussSeidel linear solver
max_iterations = 20
atol = 0.0

julia> gs(A, b)
2-element Vector{Interval{Float64}}:
 [-1.66668, 1.66668]
 [-1.33334, 1.33334]
```
"""
struct GaussSeidel <: AbstractIterativeSolver
    max_iterations::Int
    atol::Float64
end

GaussSeidel() = GaussSeidel(20, 0.0)

function (gs::GaussSeidel)(A::AbstractMatrix{T},
                           b::AbstractVector{T},
                           x::AbstractVector{T}=enclose(A, b)) where {T<:Interval}
    n = length(b)

    atol = iszero(gs.atol) ? minimum(diam.(A))*1e-5 : gs.atol
    @inbounds for _ in 1:gs.max_iterations
        xold = copy(x)
        @inbounds for i in 1:n
            x[i] = b[i]
            @inbounds for j in 1:n
                (i == j) || (x[i] -= A[i, j] * x[j])
            end
            x[i] = (x[i]/A[i, i]) .∩ xold[i]
        end
        all(interval_isapprox.(x, xold; atol=atol)) && break
    end
    return x
end

## KRAWCZYK
"""
    LinearKrawczyk <: AbstractIterativeSolver

Type for the Krawczyk solver of the interval linear system ``Ax=b``.
For details see Section 5.7.3 of [[HOR19]](@ref)

### Fields

- `max_iterations` -- maximum number of iterations (default 20)
- `atol`           -- absolute tolerance (default 0), if at some point ``|xₖ - xₖ₊₁| < atol``
                      (elementwise), then stop and return ``xₖ₊₁``.
                      If `atol=0`, then `min(diam(A))*1e-5` is used.

### Notes

- An object of type `LinearKrawczyk` is a function with method

        (kra::LinearKrawczyk)(A::AbstractMatrix{T},
                              b::AbstractVector{T},
                              [x]::AbstractVector{T}=enclose(A, b)) where {T<:Interval}

    #### Input
    - `A`   -- N×N interval matrix
    - `b`   -- interval vector of length N
    - `x`   -- (optional) initial enclosure for the solution of ``Ax = b``. If not given,
               it is automatically computed using [`enclose`](@ref enclose)

### Examples

```jldoctest
julia> A = [2..4 -1..1;-1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> b = [-2..2, -1..1]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-1, 1]

julia> kra = LinearKrawczyk()
LinearKrawczyk linear solver
max_iterations = 20
atol = 0.0

julia> kra(A, b)
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-2, 2]
```
"""
struct LinearKrawczyk <: AbstractIterativeSolver
    max_iterations::Int
    atol::Float64
end

LinearKrawczyk() = LinearKrawczyk(20, 0.0)

function (kra::LinearKrawczyk)(A::AbstractMatrix{T},
                               b::AbstractVector{T},
                               x::AbstractVector{T}=enclose(A, b)) where {T<:Interval}

    atol = iszero(kra.atol) ? minimum(diam.(A))*1e-5 : kra.atol

    C = inv(mid.(A))
    for i = 1:kra.max_iterations
        xnew  = (C*b  - C*(A*x) + x) .∩ x
        all(interval_isapprox.(x, xnew; atol=atol)) && return xnew
        x = xnew
    end
    return x
end

# custom printing for solvers
function Base.string(s::AbstractLinearSolver)

    str="""$(typeof(s)) linear solver
    """

    fields = fieldnames(typeof(s))
    for field in fields
        str *= """$field = $(getfield(s, field))
        """
    end
    return str
end

Base.show(io::IO, s::AbstractLinearSolver) = print(io, string(s))


text:
# using POMDPModelTools
const NEIGHBOR_DIST = 1 # how big a neighborhood, in terms of Euclidean distance

# # not used now
# const WIND_DIR_DICT = Dict(:north=>1,
#                         :north_east=>2,
#                         :east=>3,
#                         :south_east=>4,
#                         :south=>5,
#                         :south_west=>6,
#                         :west=>7,
#                         :north_west=>8)

function POMDPs.transition(pomdp::FireWorld, state::FireState, action::Array{Int64,1})
    total_size = pomdp.grid_size * pomdp.grid_size
    # 1. update transition of state
    sp_trans_dist = transition_of_states(pomdp, state, action)

    neighbors = FireState[]
    probabilities = Float64[]

    for sp_trans in sp_trans_dist.vals
        # 2. update fuels - other cells that were burning but no action was taking to put it out
        prob = pdf(sp_trans_dist, sp_trans)

        sp_burn = sp_trans.burning
        sp_burn_probs = sp_trans.burn_probs
        sp_fuels = sp_trans.fuels

        # get indices of cells where no put out action was applied to
        all_cells = Array[1:total_size]
        no_action_cells = all_cells[1][minus(action, all_cells[1])]

        # update fuel levels of cells that were burning yet no action was applied to
        sp_fuels_new = deepcopy(sp_fuels)
        for cell in no_action_cells
            if sp_burn[cell]
                new_fuel = sp_fuels[cell]-1
                if new_fuel < 1
                    sp_fuels_new[cell] = 0
                    sp_burn[cell] = false
                else
                    sp_fuels_new[cell] = new_fuel
                end
            end
        end
        # want to keep the burning indicator even if fuel of a cell is 0 for step 3 below
        # as that cell may spread the fire to other cells despite itself being burned to fuel exhaustion
        sp_new = FireState(sp_burn, sp_burn_probs, sp_fuels_new)

        # 3. update fire spread
        P_xy = fire_spread(pomdp, sp_new)
        sp = update_burn_probs(sp_new, P_xy)

        push!(neighbors, sp)
        push!(probabilities, prob)
    end
    return SparseCat(neighbors, normalize!(probabilities, 1))
end

# Transition helper
function transition_of_states(pomdp::FireWorld, state::FireState, action::Array{Int64,1})
    tprob = pomdp.tprob
    burning = state.burning
    burn_probs = state.burn_probs
    fuels = state.fuels
    neighbors = FireState[]
    probabilities = Array{Float64,1}(undef,0)

    # action is an array of indices that we will apply fire fighting actions to
    # those cells may or may not be burning
    # 1. find cell indices of the action set where cell is burning now
    act_burn = burning[action] .* action
    act_on_burning_cells = act_burn[findall(x->x>0, act_burn)]
    num_act_n_burning = length(act_on_burning_cells)

    # check if any cell is burning now
    if num_act_n_burning > 0
        # 1.(a) decrement fuels in those cells
        fuels_new = [i in act_on_burning_cells ? max(0, fuels[i]-1) : fuels[i] for i in 1:length(fuels) ]
        # 1.(b) account for varying fire fighting outcomes
        # some cells may be put out successfully and some may not be
        # for i in space([false, true], num_act_n_burning) # tuple is faster
        for i in space((false, true), num_act_n_burning) # different permutations of fire fighting outcomes
            burn = collect(i)
            burning_new = update_fullburn(burning, burn, act_on_burning_cells)
            prob = (1-tprob)^sum(burn)*(tprob)^(num_act_n_burning - sum(burn)) # basically binomial
            push!(neighbors, FireState(burning_new, burn_probs, fuels_new))
            push!(probabilities, prob)
        end
    else # none of the cells we apply action to is burning to begin with
        push!(neighbors, FireState(burning, burn_probs, fuels))
        push!(probabilities, 1.0)
    end
    return SparseCat(neighbors, probabilities)
end

# update burning or not for cells that were burning
function update_fullburn(burning::BitArray{1}, burn_perm::Array{Bool,1}, burn_indices::Array{Int64,1})
    burns = deepcopy(burning)
    for i in 1:length(burn_perm)
        burn_new = burn_perm[i]
        idx = burn_indices[i]
        burns[idx] = burn_new
    end
    return burns
end

# FIRE SPREAD
# returns new FireState with different probability burn map
function update_burn_probs(s::FireState, P_xy::Array{Float64,2})
    burn = s.burning
    burn_probs = s.burn_probs
    fuels = s.fuels

    total_size = length(s.burning)
    burn_probs_update = zeros(total_size)
    fuels_left = findall(x->x>0, fuels)
    for i in fuels_left
        burn_probs_update[i] = burn_probs[i]
    end

    burn_probs_result = zeros(total_size)
    for i in 1:total_size
        not_spread_i = 1
        for j in 1:total_size
            P_xy[i,j] , burn_probs_update[j]
            not_spread_j = 1 - P_xy[i,j] * burn_probs_update[j]
            not_spread_i *= not_spread_j
            burn_probs_result[i] = 1 - not_spread_i
        end
    end

    all_cells = collect(1:total_size)
    no_fuels = all_cells[minus(fuels_left, all_cells)]
    for i in no_fuels
        burn_probs_result[i] = 0
    end
    return FireState(burn, burn_probs_result, fuels)
end

function fire_spread(pomdp::FireWorld, s::FireState)
    total_size = pomdp.grid_size * pomdp.grid_size
# indices conversion
    cartesian = CartesianIndices((1:pomdp.grid_size, 1:pomdp.grid_size))
    linear = LinearIndices((1:pomdp.grid_size, 1:pomdp.grid_size))
#     "Lambda won't be input;
#     this needs to take input of distance, wind, and constant lambda_b (terrain-specifc)
#     to update lambda
#     Distance: can calculate from grid; the further, the harder to spread
#    # 2020/05/18 update: restrict to neighboring cells only
#     Wind: Rating of how strong wind is - want to get at direction and speed (relative to direction of two cells)
#     i.e. need 8 directions (for cell and wind)
#     Lambda_b: say, fuel level at the cell"
    wind = pomdp.wind
    wind_strength = wind[1]
    wind_acc = wind[2]
    wind_dir = wind[3]
    longest_dist = euclidean([1,1], [GRID_SIZE, GRID_SIZE])
    to_norm = wind_strength * wind_acc * DEFAULT_FUEL
    lambdas = zeros((total_size, total_size))
    P_xy = zeros((total_size, total_size))
    burning = s.burning
    fuels = s.fuels
    for i in 1:total_size
        fuel_level = fuels[i]
        cart_i = cartesian[i]
        for j in 1:total_size
            cart_j = cartesian[j]
            rel_pos = relative_direction(cart_i, cart_j)
            wind_factor = find_wind_dir_factor(wind_dir, rel_pos)
            distance_ij = euclidean([cart_i[1],cart_i[2]], [cart_j[1], cart_j[2]])
            if distance_ij > NEIGHBOR_DIST
                distance_ij = 0
            end
            rel_dist = distance_ij/longest_dist
            lambda_ij = wind_strength*(wind_acc*rel_dist)*wind_factor*fuel_level/to_norm
            lambdas[i,j] = lambda_ij
            P_xy[i,j]= 1 - exp(-lambda_ij)
        end
    end
    return P_xy
end

function find_wind_dir_factor(wind_dir::Int, cell_dir::Int)
    alignment = abs(wind_dir - cell_dir)
    if alignment == 0
        return 1
    elseif alignment == 1 || alignment == 7
        return 0.5
    else
        return 0 # the rest are not considered
    end
end

function relative_direction(cart_i::CartesianIndex{2}, cart_j::CartesianIndex{2})
    row_i = cart_i[1]
    col_i = cart_i[2]
    row_j = cart_j[1]
    col_j = cart_j[2]
    if row_i == row_j # same row
        if col_j > col_i # j is east of i
            return 3
        else
            return 7
        end
    elseif row_j > row_i # j is south of i
        if col_i == col_j
            return 4
        else
            return 6
        end
    elseif row_j < row_i # j is west of i
        if col_i == col_j
            return 2
        else
            return 8
        end
    else # same column
        if row_j > row_i # j is south of i
            return 5
        else
            return 1
        end
    end
end

# to get cells no actions applied to
minus(indx, x) = setdiff(1:length(x), indx)

# get permutations
space(x, n) = vec(collect(Iterators.product([x for i = 1:n]...)))


text:
@testset "Utilities" begin
    deg  = 1
    aux  = true
    red  = false
    case = "case5_spm.m"

    file  = joinpath("../test/data/matpower", case)

    # solve problem
    result = _SPM.run_sopf_iv(file, _PM.IVRPowerModel, ipopt_solver, aux=aux, deg=deg, red=red)
    @test result["termination_status"] == LOCALLY_SOLVED

    pg_coeff = pce_coeff(result, "gen", 1, "pg")
    @test isapprox(pg_coeff, [1.5000000149998045;    0.000459749295155349;   0.00010187050082378998])

    vr_coeff = pce_coeff(result, "bus", 2, "vr")
    @test isapprox(vr_coeff, [0.9902405059792688;  -0.0015349621281835933;  -0.0002292009456770499  ])

    pg_sample = sample(result, "gen", 1, "pg"; sample_size=10)
    @test length(pg_sample) == 10

    pg_density = density(result, "gen", 1, "pg"; sample_size=10)
    # @test TODO: what to test here? do we re-export KernelDensity.UnivariateKDE?

end