text:
if basename(pwd()) == "aoc"
    cd("2017/21")
end

struct RuleDict
    rules2to3::Dict{BitMatrix,BitMatrix}
    rules3to4::Dict{BitMatrix,BitMatrix}
end
RuleDict() = RuleDict(Dict{BitMatrix,BitMatrix}(), Dict{BitMatrix,BitMatrix}())

function Base.getindex(rules::RuleDict, index::BitMatrix)
    dims = size(index)
    if dims == (2, 2)
        rules.rules2to3[index]
    elseif dims == (3, 3)
        rules.rules3to4[index]
    else
        throw(KeyError(index))
    end
end

Base.isempty(rd::RuleDict) = isempty(rd.rules2to3) && isempty(rd.rules3to4)
Base.length(rd::RuleDict) = length(rd.rules2to3) + length(rd.rules3to4)

function Base.setindex!(rules::RuleDict, value::BitMatrix, index::BitMatrix)
    push!(
        if size(index) == (2, 2) && size(value) == (3, 3)
            rules.rules2to3
        elseif size(index) == (3, 3) && size(value) == (4, 4)
            rules.rules3to4
        else
            throw(ArgumentError("Can't set a pair of $(summary(index)) => $(summary(value))"))
        end,
        (
            f(g(index)) => value
            for f in (identity, rotr90, rot180, rotl90),
            g in (identity, x -> reverse(x; dims = 1))
        )...
    )
    value
end

function Base.iterate(rd::RuleDict, iter = Iterators.Stateful(Iterators.flatten((rd.rules2to3, rd.rules3to4))))
    (popfirst!(iter), iter)
end
Base.eltype(::Type{RuleDict}) = Pair{BitMatrix,BitMatrix}

function parsebitmatrix(str::AbstractString)
    BitMatrix(permutedims(hcat(([c == '#' ? true : false for c in s] for s in split(str, '/'))...)))
end

glider() = BitMatrix([0 1 0; 0 0 1; 1 1 1])

@assert parsebitmatrix("../.#") == BitMatrix([0 0; 0 1])
@assert parsebitmatrix(".#./..#/###") == glider()
@assert parsebitmatrix("#..#/..../#..#/.##.") == [
    1 0 0 1
    0 0 0 0
    1 0 0 1
    0 1 1 0
]

function loadrules(filename::AbstractString)
    rules = RuleDict()
    for line in eachline(filename)
        left, right = split(line, " => ")
        rules[parsebitmatrix(left)] = parsebitmatrix(right)
    end
    rules
end

function hashdot(bits::BitMatrix)
    '\n' * join((join(bit ? '#' : '.' for bit in row) for row in eachrow(bits)), '\n')
end

function advance(rules::RuleDict, bits::BitMatrix)
    offset = 1 + isodd(size(bits, 1))
    inds = [x:x+offset for x in firstindex(bits):(offset+1):size(bits, 1)]
    vcat((hcat((rules[bits[ys, xs]] for xs in inds)...) for ys in inds)...)
end

part1(filename::AbstractString, iters = 5) = part1(loadrules(filename), iters)
function part1(rules::RuleDict, iters = 5)
    bits = glider()
    @debug hashdot(bits)
    for _ in 1:iters
        bits = advance(rules, bits)
        @debug hashdot(bits)
    end
    count(bits)
end

@assert part1("example.txt", 2) == 12
part1("input.txt", 5)
part1("input.txt", 18)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
tests = ["deque",
		 "stack_and_queue",
		 "accumulator",
		 "classifiedcollections",
		 "disjoint_set",
		 "binheap",
		 "mutable_binheap",
		 "defaultdict",
		 "ordereddict",
		 "orderedset",
		 "trie",
		 "list"]

for t in tests
	fp = joinpath(Pkg.dir("DataStructures"), "test", "test_$t.jl")
	println("$fp ...")
	include(fp)
end

ext: jl


text:
struct GroupStats{N, U, S}
    keyvars::NTuple{N, StatVariate}
    covariate::Union{AbstractCovariate, Nothing}
    selector::Union{AbstractBoolVariate, Nothing}
    stats::Dict{NTuple{N, U}, S}
end

function getgroupstats(factors::NTuple{N, Union{AbstractFactor, AbstractIntVariate}}; selector::AbstractBoolVariate = BoolVariate("", BitArray{1}()), slicelength::Integer = SLICELENGTH) where {N}
    U = promote_type(map((s -> eltype(s)), factors)...)
    u = zero(U)
    fromobs = 1
    toobs = length(factors[1])
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    noselector = length(selector) == 0
    if N == 1
        slices = slice(factors[1], fromobs, toobs, slicelength)
    elseif N == 2
        slices = zip(slice(factors[1], fromobs, toobs, slicelength), slice(factors[2], fromobs, toobs, slicelength))
    elseif N == 3
        slices = zip(slice(factors[1], fromobs, toobs, slicelength), slice(factors[2], fromobs, toobs, slicelength), slice(factors[3], fromobs, toobs, slicelength))
    elseif N == 4
        slices = zip(slice(factors[1], fromobs, toobs, slicelength), slice(factors[2], fromobs, toobs, slicelength), slice(factors[3], fromobs, toobs, slicelength), slice(factors[4], fromobs, toobs, slicelength))
    else
        slices = zip(map((s -> slice(s, fromobs, toobs, slicelength)), factors))
    end
    if noselector
        dict = fold(Dict{NTuple{N, U}, Int64}(), slices) do d, slice
            if N == 1
                @inbounds for i in 1:length(slice)
                    v = (slice[i], )
                    d[v] = get(d, v, 0) + 1
                end
            elseif N == 2
                slice1, slice2 = slice
                @inbounds for i in 1:length(slice1)
                    v = oftype(u, slice1[i]), oftype(u, slice2[i])
                    d[v] = get(d, v, 0) + 1
                end
            elseif N == 3
                slice1, slice2, slice3 = slice
                @inbounds for i in 1:length(slice1)
                    v = oftype(u, slice1[i]), oftype(u, slice2[i]), oftype(u, slice3[i])
                    d[v] = get(d, v, 0) + 1
                end
            elseif N == 4
                slice1, slice2, slice3, slice4 = slice
                @inbounds for i in 1:length(slice1)
                    v = oftype(u, slice1[i]), oftype(u, slice2[i]), oftype(u, slice3[i]), oftype(u, slice4[i])
                    d[v] = get(d, v, 0) + 1
                end
            else
                for i in 1:length(slice[1])
                    v = map((x -> oftype(u, x[i])), slice)
                    d[v] = get(d, v, 0) + 1
                end
            end
            d
        end
        GroupStats{N, U, Int64}(factors, nothing, nothing, dict)
    else
        selslices = slice(selector, fromobs, toobs, slicelength)
        zipslices = zip(selslices, slices)
        dict = fold(Dict{NTuple{N, U}, Int64}(), zipslices) do d, zipslice
            selslice, slice = zipslice
            if N == 1
                @inbounds for i in 1:length(slice)
                    if selslice[i]
                        v = (slice[i], )
                        d[v] = get(d, v, 0) + 1
                    end
                end
            elseif N == 2
                slice1, slice2 = slice
                @inbounds for i in 1:length(slice1)
                    if selslice[i]
                        v = oftype(u, slice1[i]), oftype(u, slice2[i])
                        d[v] = get(d, v, 0) + 1
                    end
                end
            elseif N == 3
                slice1, slice2, slice3 = slice
                @inbounds for i in 1:length(slice1)
                    v = oftype(u, slice1[i]), oftype(u, slice2[i]), oftype(u, slice3[i])
                    d[v] = get(d, v, 0) + 1
                end
            elseif N == 4
                slice1, slice2, slice3, slice4 = slice
                @inbounds for i in 1:length(slice1)
                    if selslice[i]
                        v = oftype(u, slice1[i]), oftype(u, slice2[i]), oftype(u, slice3[i]), oftype(u, slice4[i])
                        d[v] = get(d, v, 0) + 1
                    end
                end
            else
                for i in 1:length(slice[1])
                    if selslice[i]
                        v = map((x -> oftype(u, x[i])), slice)
                        d[v] = get(d, v, 0) + 1
                    end
                end
            end
            d
        end
        GroupStats{N, U, Int64}(factors, nothing, selector, dict)
    end
end

function getgroupstats(factors::Union{AbstractFactor, AbstractIntVariate}...; selector::AbstractBoolVariate = BoolVariate("", BitArray{1}()), slicelength::Integer = SLICELENGTH)
    getgroupstats(factors; selector = selector, slicelength = slicelength)
end

function getgroupstats(cov::AbstractCovariate{S}, factors::NTuple{N, AbstractFactor}; slicelength::Integer = SLICELENGTH) where {N} where {S<:AbstractFloat}
    U = promote_type(map((s -> eltype(s)), factors)...)
    u = zero(U)
    fromobs = 1
    toobs = length(factors[1])
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    if N == 1
        slices = slice(factors[1], fromobs, toobs, slicelength)
    elseif N == 2
        slices = zip(slice(factors[1], fromobs, toobs, slicelength), slice(factors[2], fromobs, toobs, slicelength))
    elseif N == 3
        slices = zip(slice(factors[1], fromobs, toobs, slicelength), slice(factors[2], fromobs, toobs, slicelength), slice(factors[3], fromobs, toobs, slicelength))
    elseif N == 4
        slices = zip(slice(factors[1], fromobs, toobs, slicelength), slice(factors[2], fromobs, toobs, slicelength), slice(factors[3], fromobs, toobs, slicelength), slice(factors[4], fromobs, toobs, slicelength))
    else
        slices = zip(map((s -> slice(s, fromobs, toobs, slicelength)), factors))
    end

    covslices = slice(cov, fromobs, toobs, slicelength)
    zipslices = zip(slices, covslices)
    dict = fold(Dict{NTuple{N, U}, CovariateStats}(), zipslices) do d, zipslice
        slice, covslice = zipslice

        if N == 1
            @inbounds for i in 1:length(slice)
                v = covslice[i]
                y = (slice[i], )
                if !(y in keys(d))
                    covstats = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
                    d[y] = covstats
                else
                    covstats = d[y]
                end
                covstats.obscount += 1
                if isnan(v)
                    covstats.nancount += 1
                else
                    if isnan(covstats.sum)
                        covstats.sum = v
                        covstats.sum2 = v * v
                        covstats.min = v
                        covstats.max = v
                    else
                        covstats.sum += v
                        covstats.sum2 += v * v
                        if v < covstats.min
                            covstats.min = v
                        end
                        if v > covstats.max
                            covstats.max = v
                        end
                    end
                end
            end
        elseif N == 2
            slice1, slice2 = slice
            @inbounds for i in 1:length(slice1)
                v = covslice[i]
                y = oftype(u, slice1[i]), oftype(u, slice2[i])
                if !(y in keys(d))
                    covstats = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
                    d[y] = covstats
                else
                    covstats = d[y]
                end
                covstats.obscount += 1
                if isnan(v)
                    covstats.nancount += 1
                else
                    if isnan(covstats.sum)
                        covstats.sum = v
                        covstats.sum2 = v * v
                        covstats.min = v
                        covstats.max = v
                    else
                        covstats.sum += v
                        covstats.sum2 += v * v
                        if v < covstats.min
                            covstats.min = v
                        end
                        if v > covstats.max
                            covstats.max = v
                        end
                    end
                end
            end
        elseif N == 3
            slice1, slice2, slice3 = slice
            @inbounds for i in 1:length(slice1)
                v = covslice[i]
                y = oftype(u, slice1[i]), oftype(u, slice2[i]), oftype(u, slice3[i])
                if !(y in keys(d))
                    covstats = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
                    d[y] = covstats
                else
                    covstats = d[y]
                end
                covstats.obscount += 1
                if isnan(v)
                    covstats.nancount += 1
                else
                    if isnan(covstats.sum)
                        covstats.sum = v
                        covstats.sum2 = v * v
                        covstats.min = v
                        covstats.max = v
                    else
                        covstats.sum += v
                        covstats.sum2 += v * v
                        if v < covstats.min
                            covstats.min = v
                        end
                        if v > covstats.max
                            covstats.max = v
                        end
                    end
                end
            end
        elseif N == 4
            slice1, slice2, slice3, slice4 = slice
            @inbounds for i in 1:length(slice1)
                v = covslice[i]
                y = oftype(u, slice1[i]), oftype(u, slice2[i]), oftype(u, slice3[i]), oftype(u, slice4[i])
                if !(y in keys(d))
                    covstats = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
                    d[y] = covstats
                else
                    covstats = d[y]
                end
                covstats.obscount += 1
                if isnan(v)
                    covstats.nancount += 1
                else
                    if isnan(covstats.sum)
                        covstats.sum = v
                        covstats.sum2 = v * v
                        covstats.min = v
                        covstats.max = v
                    else
                        covstats.sum += v
                        covstats.sum2 += v * v
                        if v < covstats.min
                            covstats.min = v
                        end
                        if v > covstats.max
                            covstats.max = v
                        end
                    end
                end
            end
        else
            for i in 1:length(slice[1])
                v = covslice[i]
                y = map((x -> oftype(u, x[i])), slice)
                if !(y in keys(d))
                    covstats = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
                    d[y] = covstats
                else
                    covstats = d[y]
                end
                covstats.obscount += 1
                if isnan(v)
                    covstats.nancount += 1
                else
                    if isnan(covstats.sum)
                        covstats.sum = v
                        covstats.sum2 = v * v
                        covstats.min = v
                        covstats.max = v
                    else
                        covstats.sum += v
                        covstats.sum2 += v * v
                        if v < covstats.min
                            covstats.min = v
                        end
                        if v > covstats.max
                            covstats.max = v
                        end
                    end
                end
            end
        end
        d
    end
    for (_, stats) in dict
        stats.nanpcnt = 100.0 * stats.nancount / stats.obscount
        stats.mean = stats.sum / (stats.obscount - stats.nancount)
        stats.std = sqrt(((stats.sum2 - stats.sum * stats.sum / (stats.obscount - stats.nancount)) / (stats.obscount - stats.nancount - 1)))
    end
    GroupStats{N, U, CovariateStats}(factors, cov, nothing, dict)
end

function getgroupstats(cov::AbstractCovariate{S}, factors::AbstractFactor...; slicelength::Integer = SLICELENGTH) where {S<:AbstractFloat}
    getgroupstats(cov, factors; slicelength = slicelength)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using Formatting

function scratchmax(arr)
    max_value=arr[1]
    for a in arr
        if max_value < a
            max_value=a
        end
    end
    return max_value

end

function main()
    N=100000000
    arr=[x for x in 0:N-1]
    @time maximum(arr)
    @time max_value=scratchmax(arr)
    printfmt("max value of the array is {}\n", max_value)
end

main()

max_stars_repo_name: terasakisatoshi/pythonCodes

max_issues_repo_name: terasakisatoshi/pythonCodes

max_forks_repo_name: terasakisatoshi/pythonCodes


text:
## Makes an adaptive Metropolis kernel proposed by:
## Haario, H., Saksman, E. and Tamminen, J., 2001. An adaptive Metropolis
## algorithm. Bernoulli, 7(2), pp.223-242.
## Essentially provides an adaptive mechanism to make use of the optimal
## scaling results for random walk Metropolis on "Gaussian-like" targets
## described by:
## Roberts, G.O. and Rosenthal, J.S., 2001. Optimal scaling for various
## Metropolis--Hastings algorithms. Statistical Science, 16(4), pp.351-367.
function makeAMKernel(logTargetDensity::F, Σ::SMatrix{d, d, Float64},
  rng::RNG = Random.GLOBAL_RNG, updateFrequency::Int64 = 1, ϵ::Float64 = 1.0) where
  {F<:Function, d, RNG<:AbstractRNG}

  S::MMatrix{d, d, Float64} = Σ
  A::MMatrix{d, d, Float64} = cholesky(S).L

  scratchv::MVector{d, Float64} = MVector{d, Float64}(undef)
  scratchz::MVector{d, Float64} = MVector{d, Float64}(undef)
  prevx::MVector{d, Float64} = MVector{d, Float64}(undef)
  ldprevx = Ref(-Inf)

  accepts = Ref(0)
  calls = Ref(0)
  covEstimate::MMatrix{d, d, Float64} = Σ
  meanEstimate::MVector{d, Float64} = zeros(MVector{d, Float64})

  @inline function retuneSigma()
    if covEstimate[1,1] == 0.0
      S .= Σ  * ϵ / calls.x
    else
      S .= 5.6644 / d * covEstimate
    end
    try
      A .= cholesky(S).L
      any(isnan, A) && throw(DomainError()) # fix for not throwing exception
    catch e
      S .= Σ  * ϵ / calls.x
      A .= cholesky(S).L
    end
  end

  @inline function P(x::SVector{d, Float64})
    calls.x += 1
    randn!(rng, scratchv)
    mul!(scratchz, A, scratchv)
    z::SVector{d, Float64} = scratchz + x

    if x == prevx
      lpi_x = ldprevx.x
    else
      lpi_x = logTargetDensity(x)
      prevx .= x
    end
    lpi_z = logTargetDensity(z)
    if -randexp(rng) < lpi_z - lpi_x
      prevx .= z
      ldprevx.x = lpi_z
      accepts.x += 1
      rval = z
    else
      rval = x
    end
    t::Int64 = calls.x
    covEstimate .= (t-1)/t * (covEstimate +
      (rval - meanEstimate) * (rval - meanEstimate)' / t)
    meanEstimate .= (t-1)/t.*meanEstimate + rval/t
    mod(t, updateFrequency) == 0 && retuneSigma()
    return rval
  end

  @inline function P(s::Symbol)
    s == :acceptanceRate && return accepts.x / calls.x
    s == :meanEstimate && return meanEstimate
    s == :covEstimate && return covEstimate * calls.x /(calls.x-1)
  end

  return P
end

function makeAMKernel(logTargetDensity::F, Σ::Matrix{Float64},
  rng::RNG = Random.GLOBAL_RNG, updateFrequency::Int64 = 1, ϵ::Float64 = 1.0) where
  {F<:Function, RNG<:AbstractRNG}

  d = size(Σ)[1]
  S::Matrix{Float64} = copy(Σ)
  A::Matrix{Float64} = cholesky(S).L

  scratchv::MVector{d, Float64} = MVector{d, Float64}(undef)
  scratchz::MVector{d, Float64} = MVector{d, Float64}(undef)
  prevx::MVector{d, Float64} = MVector{d, Float64}(undef)
  ldprevx = Ref(-Inf)

  accepts = Ref(0)
  calls = Ref(0)
  covEstimate::Matrix{Float64} = copy(Σ)
  meanEstimate::MVector{d, Float64} = zeros(MVector{d, Float64})

  @inline function retuneSigma()
    if covEstimate[1,1] == 0.0
      S .= Σ  * ϵ / calls.x
    else
      S .= 5.6644 / d * covEstimate
    end
    try
      A .= cholesky(S).L
      any(isnan, A) && throw(DomainError()) # fix for not throwing exception
    catch e
      S .= Σ  * ϵ / calls.x
      A .= cholesky(S).L
    end
  end

  @inline function P(x::SVector{d, Float64}) where d
    calls.x += 1
    randn!(rng, scratchv)
    mul!(scratchz, A, scratchv)
    z::SVector{d, Float64} = scratchz + x

    if x == prevx
      lpi_x = ldprevx.x
    else
      lpi_x = logTargetDensity(x)
      prevx .= x
    end
    lpi_z = logTargetDensity(z)
    if -randexp(rng) < lpi_z - lpi_x
      prevx .= z
      ldprevx.x = lpi_z
      accepts.x += 1
      rval = z
    else
      rval = x
    end
    t::Int64 = calls.x
    covEstimate .= (t-1)/t * (covEstimate +
      (rval - meanEstimate) * (rval - meanEstimate)' / t)
    meanEstimate .= (t-1)/t.*meanEstimate + rval/t
    mod(t, updateFrequency) == 0 && retuneSigma()
    return rval
  end

  @inline function P(s::Symbol)
    s == :acceptanceRate && return accepts.x / calls.x
    s == :meanEstimate && return meanEstimate
    s == :covEstimate && return covEstimate * calls.x /(calls.x-1)
  end

  return P
end

function makeAMKernel(logTargetDensity::F, d::Int64, rng::RNG=Random.GLOBAL_RNG,
  updateFrequency::Int64 = 1, ϵ::Float64 = 1.0) where
  {F<:Function, RNG<:AbstractRNG}
  Id = SMatrix{d, d, Float64}(Matrix(1.0I, d, d))
  return makeAMKernel(logTargetDensity::F, Id, rng, updateFrequency, ϵ)
end


text:
#####
##### generate plot for readme
#####
##### - its directory should be the current one
##### - PGFPlotsX should be available in the environment

using NonhomotheticCES, PGFPlotsX, StaticArrays

U = NonhomotheticCESUtility(0.5, SVector(0.0, 0.0), SVector(1.0, 2.0))

E = range(0.001, 0.5; length = 100)
p̂ = SVector(0.1, 0.0)
Ĉ = [log_consumption_aggregator(U, p̂, log(E)) for E in E]
ĉs = [log_sectoral_consumptions(U, p̂, log(E), Ĉ) for (E, Ĉ) in zip(E, Ĉ)]

pgfsave(joinpath(@__DIR__, "example.png"),
        @pgf Axis({ xlabel = "expenditure", legend_pos = "north west",
                    yticklabel_style={ "/pgf/number format/fixed",
                                       "/pgf/number format/precision" = 2 },
                    scaled_y_ticks = false,
                    },
                  Plot({ no_marks, black }, Table(E, exp.(Ĉ))),
                  LegendEntry(raw"$C$"),
                  Plot({ no_marks, red }, Table(E, exp.(first.(ĉs)))),
                  LegendEntry(raw"$c_1$"),
                  Plot({ no_marks, blue }, Table(E, exp.(last.(ĉs)))),
                  LegendEntry(raw"$c_2$"));
        dpi = 600)


text:
"""
    Ai2{T}

`Ai2` performs over-approximated reachability analysis to compute the over-approximated
output reachable set for a network. `T` can be `Hyperrectangle`, `Zonotope`, or
`HPolytope`, and determines the amount of over-approximation (and hence also performance
tradeoff). The original implementation (from [1]) uses Zonotopes, so we consider this
the "benchmark" case. The `HPolytope` case is more precise, but slower, and the opposite
is true of the `Hyperrectangle` case.

Note that initializing `Ai2()` defaults to `Ai2{Zonotope}`.
The following aliases also exist for convenience:

```julia
const Ai2h = Ai2{HPolytope}
const Ai2z = Ai2{Zonotope}
const Box = Ai2{Hyperrectangle}
```

# Problem requirement
1. Network: any depth, ReLU activation (more activations to be supported in the future)
2. Input: AbstractPolytope
3. Output: AbstractPolytope

# Return
`ReachabilityResult`

# Method
Reachability analysis using split and join.

# Property
Sound but not complete.

# Reference
[1] T. Gehr, M. Mirman, D. Drashsler-Cohen, P. Tsankov, S. Chaudhuri, and M. Vechev,
"Ai2: Safety and Robustness Certification of Neural Networks with Abstract Interpretation,"
in *2018 IEEE Symposium on Security and Privacy (SP)*, 2018.

## Note
Efficient over-approximation of intersections and unions involving zonotopes relies on Theorem 3.1 of

[2] Singh, G., Gehr, T., Mirman, M., Püschel, M., & Vechev, M. (2018). Fast
and effective robustness certification. In Advances in Neural Information
Processing Systems (pp. 10802-10813).
"""
struct Ai2{T<:Union{Hyperrectangle, Zonotope, HPolytope}} <: Solver end

Ai2() = Ai2{Zonotope}()
const Ai2h = Ai2{HPolytope}
const Ai2z = Ai2{Zonotope}
const Box = Ai2{Hyperrectangle}

function solve(solver::Ai2, problem::Problem)
    reach = forward_network(solver, problem.network, problem.input)
    return check_inclusion(reach, problem.output)
end

forward_layer(solver::Ai2, L::Layer, inputs::Vector) = forward_layer.(solver, L, inputs)

function forward_layer(solver::Ai2h, L::Layer{ReLU}, input::AbstractPolytope)
    Ẑ = affine_map(L, input)
    relued_subsets = forward_partition(L.activation, Ẑ) # defined in reachability.jl
    return convex_hull(UnionSetArray(relued_subsets))
end

# method for Zonotope and Hyperrectangle, if the input set isn't a Zonotope
function forward_layer(solver::Union{Ai2z, Box}, L::Layer{ReLU}, input::AbstractPolytope)
    return forward_layer(solver, L, overapproximate(input, Hyperrectangle))
end

function forward_layer(solver::Ai2z, L::Layer{ReLU}, input::AbstractZonotope)
    Ẑ = affine_map(L, input)
    return overapproximate(Rectification(Ẑ), Zonotope)
end

function forward_layer(solver::Box, L::Layer{ReLU}, input::AbstractZonotope)
    Ẑ = approximate_affine_map(L, input)
    return rectify(Ẑ)
end

function forward_layer(solver::Ai2, L::Layer{Id}, input)
    return affine_map(L, input)
end

function convex_hull(U::UnionSetArray{<:Any, <:HPolytope})
    tohrep(VPolytope(LazySets.convex_hull(U)))
end


text:
struct Network
    gateway::String
    ipaddress::String
    netmask::Union{String, Integer}  # ipv6 netmasks are integers
    networktype::String

    function Network(data::Dict{String})
        new(
            data["gateway"],
            data["ip_address"],
            data["netmask"],
            data["type"]
        )
    end
end

function show(io::IO, n::Network)
    print(io, "Network ($(n.ipaddress))")
end

struct Kernel
    id::Integer
    name::String
    version::String

    function Kernel(data::Dict{String})
        new(
            data["id"],
            data["name"],
            data["version"]
        )
    end
end

function show(io::IO, k::Kernel)
    print(io, "Kernel ($(k.name)-$(k.version))")
end

struct Droplet
    id::Integer
    name::String
    memory::Integer
    vcpus::Integer
    disk::Integer
    locked::Bool
    created_at::DateTime
    status::String
    backupids::Array{Integer, 1}
    snapshotids::Array{Integer, 1}
    features::Array{String, 1}
    region::Union{Nothing, Region}
    image::Union{Nothing, Image}
    size::Union{Nothing, Size}
    sizeslug::String
    networks::Dict{String, Array{Network, 1}}
    kernel::Union{Nothing, Kernel}
    nextbackupwindow::Union{Nothing, Dict{String, DateTime}}
    tags::Array{String, 1}
    volumeids::Array{Integer, 1}

    function Droplet(data::Dict{String})
        # we assume all DO datetimes are in UTC
        data["created_at"] = DateTime(data["created_at"][1:end-1])

        if isempty(data["region"])
            data["region"] = nothing
        else
            data["region"] = Region(data["region"])
        end
        if isempty(data["image"])
            data["image"] = nothing
        else
            data["image"] = Image(data["image"])
        end
        if isempty(data["size"])
            data["size"] = nothing
        else
            data["size"] = Size(data["size"])
        end

        if haskey(data["networks"], "v4")
            networks = Array{Network, 1}(UndefInitializer(), length(data["networks"]["v4"]))
            for (j, network) in enumerate(data["networks"]["v4"])
                networks[j] = Network(network)
            end
            data["networks"]["v4"] = networks
        end
        if haskey(data["networks"], "v6")
            networks = Array{Network, 1}(UndefInitializer(), length(data["networks"]["v6"]))
            for (j, network) in enumerate(data["networks"]["v6"])
                networks[j] = Network(network)
            end
            data["networks"]["v6"] = networks
        end
        if data["kernel"] != nothing
            data["kernel"] = Kernel(data["kernel"])
        end
        if !haskey(data, "next_backup_window")
            data["next_backup_window"] = nothing
        end
        if data["next_backup_window"] != nothing
            # we assume all DO datetimes are in UTC
            bst  = DateTime(data["next_backup_window"]["start"][1:end-1])
            bend = DateTime(data["next_backup_window"]["end"][1:end-1])
            data["next_backup_window"] = Dict("start" => bst, "end" => bend)
        end
        if !haskey(data, "tags")
            data["tags"] = String[]
        end

        new(
            data["id"],
            data["name"],
            data["memory"],
            data["vcpus"],
            data["disk"],
            data["locked"],
            data["created_at"],
            data["status"],
            data["backup_ids"],
            data["snapshot_ids"],
            data["features"],
            data["region"],
            data["image"],
            data["size"],
            data["size_slug"],
            data["networks"],
            data["kernel"],
            data["next_backup_window"],
            data["tags"],
            data["volume_ids"]
        )
    end
end

function show(io::IO, d::Droplet)
    print(io, "Droplet ($(d.name))")
end

struct DropletBackup
    id::Integer
    name::String
    backuptype::String
    distribution::String
    slug::Union{Nothing, String}
    public::Bool
    regions::Array{String, 1}
    mindisksize::Integer
    created_at::DateTime

    function DropletBackup(data::Dict{String})
        data["created_at"] = DateTime(data["created_at"][1:end-1])

        new(
            data["id"],
            data["name"],
            data["type"],
            data["distribution"],
            data["slug"],
            data["public"],
            data["regions"],
            data["min_disk_size"],
            data["created_at"]
        )
    end
end

function show(io::IO, d::DropletBackup)
    if d.backuptype == "snapshot"
        print(io, "Snapshot ($(d.name))")
    else
        print(io, "Backup ($(d.name))")
    end
end

# List all droplets
function getalldroplets!(client::AbstractClient)
    uri = joinpath(ENDPOINT, "droplets?per_page=$MAXOBJECTS")
    getalldata!(client, uri, Droplet)
end

# List droplets by tag
function getdropletsbytag!(client::AbstractClient; tag::String)
    uri = joinpath(ENDPOINT, "droplets?tag_name=$(tag)&per_page=$MAXOBJECTS")
    getalldata!(client, uri, Droplet)
end

# Retrieve an existing droplet by ID
function getdroplet!(client::AbstractClient, droplet_id::Integer)
    uri = joinpath(ENDPOINT, "droplets", "$(droplet_id)")
    Droplet(getdata!(client, uri))
end

function getdroplet!(client::AbstractClient, droplet::Droplet)
    getdroplet!(client, droplet.id)
end

# List all available kernels for a droplet
function getalldropletkernels!(client::AbstractClient, dropletid::Integer)
    uri = joinpath(ENDPOINT, "droplets", "$(dropletid)", "kernels")
    getalldata!(client, uri, Kernel)
end

function getalldropletkernels!(client::AbstractClient, droplet::Droplet)
    getalldropletkernels!(client, droplet.id)
end

# Create a new droplet
function createdroplet!(client::AbstractClient; name::String, region::Union{Region, String},
                        size::String, image::String, kwargs...)
    if isa(region, Region)
        region = region.slug
    end
    if isa(size, Size)
        size = size.slug
    end
    if isa(image, Image)
        if image.slug != nothing
            image = image.slug
        else
            image = image.id
        end
    end

    postbody = Dict{String, Any}([String(k[1]) => k[2] for k in kwargs])
    merge!(postbody, Dict{String, Any}("region" => region, "size" => size, "image" => image))

    uri = joinpath(ENDPOINT, "droplets")
    Droplet(postdata!(client, uri, postbody))
end

# Create multiple droplets
function createdroplets!(client::AbstractClient; names::Array{String}, region::Union{Region, String},
                        size::String, image::String, kwargs...)
    if isa(region, Region)
        region = region.slug
    end
    if isa(size, Size)
        size = size.slug
    end
    if isa(image, Image)
        if image.slug != nothing
            image = image.slug
        else
            image = image.id
        end
    end

    postbody = Dict{String, Any}([String(k[1]) => k[2] for k in kwargs])
    merge!(postbody, Dict{String, Any}("region" => region, "size" => size, "image" => image))

    uri = joinpath(ENDPOINT, "droplets")
    Droplet(postdata!(client, uri, postbody))
end

# List snapshots for a droplet
function getalldropletsnapshots!(client::AbstractClient, dropletid::Integer)
    uri = joinpath(ENDPOINT, "droplets", "$dropletid", "snapshots?per_page=$MAXOBJECTS")
    getalldata!(client, uri, DropletBackup)
end

function getalldropletsnapshots!(client::AbstractClient, droplet::Droplet)
    getalldropletsnapshots!(client, droplet.id)
end

# List snapshots for a droplet
function getalldropletbackups!(client::AbstractClient, dropletid::Integer)
    uri = joinpath(ENDPOINT, "droplets", "$dropletid", "backups?per_page=$MAXOBJECTS")
    getalldata!(client, uri, DropletBackup)
end

function getalldropletbackups!(client::AbstractClient, droplet::Droplet)
    getalldropletbackups!(client, droplet.id)
end

# Delete a droplet
function deletedroplet!(client::AbstractClient, dropletid::Integer)
    uri = joinpath(ENDPOINT, "droplets", "$(dropletid)")
    deletedata!(client, uri)
end

function deletedroplet!(client::AbstractClient, droplet::Droplet)
    deletedroplet!(client, droplet.id)
end

# Delete droplets by tag
function deletedroplet!(client::AbstractClient, tag::String)
    uri = joinpath(ENDPOINT, "droplets?tag_name=$(tag)")
    deletedata!(client, uri)
end

# List neighbors for a droplet
function getdropletneighbors!(client::AbstractClient, dropletid::Integer)
    uri = joinpath(ENDPOINT, "droplets", "$dropletid", "neighbors")
    getalldata!(client, uri, Droplet)
end

function getdropletneighbors!(client::AbstractClient, droplet::Droplet)
    getdropletneighbors!(client, droplet.id)
end

# List all droplet neighbors
function getalldropletneighbors!(client)
    uri = joinpath(ENDPOINT, "reports", "droplet_neighbors")
    getalldata!(client, uri, Array{Droplet, 1})
end


text:
start = map(str -> [parse(Int,str[l]) for l in findall(r"[[:digit:]]+", str)], readlines("inputs/6-10/06.txt"))[1].+1

function day6(start, days)
  c = zeros(Int, 9)
  for x in start
    c[x] += 1
  end
  for _ in 1:days
    n = zeros(Int, 9)
    for (age, cnt) in enumerate(c)
      if age == 1
        n[7] += cnt
        n[9] += cnt
      else
        n[age-1] += cnt
      end
    end
    c = n
  end
  sum(c)
end

println("Day 6 Part 1: ",day6(start, 80))
println("Day 6 Part 2: ",day6(start, 256))


text:
#using Revise, JLD
include("GL2017Replication.jl")
using .GL2017Replication

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

# 1) Load model structure
gl = ModelGL()
describe(gl)

# 2) solve model with default parameters  -- requires 1)
compute_steady_state!(gl)
describe(gl)

# 3) calibrate to target values for initial ss  -- requires 1)
calibrate!(gl,load_params=true) #pre-load parameters to speed up calibration
describe(gl)

# 4) calibrate to target values for terminal ss -- requires 3)
gl_tss = calibrate_terminal(gl,load_params=true) #pre-load parameters to speed up calibration
describe(gl,gl_tss)

# 5) compute transition dynamics -- requires 3) and 4)
Tgl      = TransGL()
gl_trans = transition!(gl,gl_tss,Tgl)

# 6) Plotting
(F1,_,_)        = plots_figure_1(gl )
(F3,_,_,_,_)    = plots_figure_3(gl,Tgl )
(F4,_,)         = plots_figure_4(gl,gl_tss )

# -----------------------------------------------------------------------------
# plot results (same as above)

using Plots
plotly(linewidth = 2.)
using Parameters, Statistics

@unpack b_grid,pr = gl                  # objects that are identical across both steady states
@unpack r_t,ϕ_t,Y_t,D_4Y_t = Tgl        # transition objects

# Figure I (from paper)

Y1 = gl.Y_actual
plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.c_pol[2,b_grid .>= -gl.ϕ],label = "θ = 2")
plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.c_pol[8,b_grid .>= -gl.ϕ],
        linestyle = :dash,label = "θ = 8")
title!("consumption")
xaxis!([-gl.ϕ,12.5])
#savefig("images\\c_pol_iss.png")

plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.n_pol[2,gl.b_grid .>= -gl.ϕ],label =  "θ = 2")
plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.n_pol[8,gl.b_grid .>= -gl.ϕ],
        linestyle = :dash,label = "θ = 8")
title!("labor supply")
xaxis!([-gl.ϕ,12.5])
#savefig("images\\l_pol_iss.png")

# Figure IV (from paper)

b_acc1 =  ((pr'*gl.b_pol[:,b_grid     .>= -gl.ϕ])'     .- b_grid[b_grid .>= -gl.ϕ])/(4*Y1)
b_acc2 =  ((pr'*gl_tss.b_pol[:,b_grid .>= -gl.ϕ])' .- b_grid[b_grid .>= -gl.ϕ])/(4*Y1)
plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1), b_acc1)
plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1), b_acc2,
        linestyle = :dash)
title!("bond accumulation policy")
xaxis!([-2,14],-2:2:14)
yaxis!([-0.4,0.6],-0.4:0.2:0.6)
#savefig("b_acc.png")

bond_distribution1 = dropdims( pr'    *gl.JD[:,b_grid .>= -gl.ϕ]     ,dims=1)
bond_distribution2 = dropdims( pr'*gl_tss.JD[:,b_grid .>= -gl.ϕ] ,dims=1)
plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1),  bond_distribution1,label = "initial")
plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1),  bond_distribution2,label = "terminal",
        linestyle = :dash)
title!("bond distribution")
xaxis!([-2,14],-2:2:14)
yaxis!([0,0.004])
xaxis!([-gl.ϕ,12.5])
#savefig("b_dist.png")

# Figure III (from paper)
Tp = 24    # number of periods plotted

p1 = plot(0:Tp, Tgl.ϕ_t[1:Tp+1]./(4*Tgl.Y_t[1:Tp+1]),ylims = (0.5,1), yticks = 0.5:0.1:10, legend = false, title = "borrowing limit")         # borrowing limit
p2 = plot(0:Tp, Tgl.D_4Y_t[1:Tp+1], ylims = (0.08,0.2), yticks = 0:0.02:0.2, legend = false, title = "household debt-to-GDP ratio")             # debt2gdp ratio
p3 = plot(0:Tp, [gl.r;Tgl.r_t[1:Tp]].*400, ylims = (-2,2.5), yticks = -2:0.5:2, legend = false, title = "interest rate")       # annualized interest rate
p4 = plot(0:Tp, [0, 100*(Tgl.Y_t[1:Tp+1]./Y1.-1)], ylims = (-1.2,0), yticks = -1:0.2:0, legend = false, title = "output")   # output deviation from steady state
plot(p1, p2, p3, p4, layout = (2, 2), legend = false)
savefig("images\\trans.png")


text:
module SFGAnalysis

using Distributions
using QuadGK

include("types.jl")
include("general.jl")
include("fresnel.jl")
include("susceptibility.jl")
include("spectra.jl")

export # types.jl
       Setup,
       # general.jl
       sfgintensity,
       sfangle,
       wl2freq,
       freq2wl,
       # fresnel.jl
       angleofrefraction,
       fresnel_x,
       fresnel_y,
       fresnel_z,
       ninterface,
       # susceptibility.jl
       susceptibility,
       effective_susceptibility_ppp,
       effective_susceptibility_ssp,
       effective_susceptibility,
       # spectra.jl
       sfspectrum,
       sfspectrum!,
       sfspectrum_voigt

end


text:
using SpinWaveTheory
using Plots: plot, plot!, savefig
using QuantumLattices: Lattice, Point, PID, Hilbert, Spin, SpinTerm, @heisenberg_str, Algorithm, ReciprocalPath, @rectangle_str, atol
using TightBindingApproximation: EnergyBands

@testset "SquareFM" begin
    lattice = Lattice(:Square,
        [Point(PID(1), [0.0, 0.0])],
        vectors=[[1.0, 0.0], [0.0, 1.0]],
        neighbors=1
        )
    hilbert = Hilbert(pid=>Spin{1//2}(1) for pid in lattice.pids)
    J = SpinTerm(:J, -1.0, 1, heisenberg"xyz")
    ms = MagneticStructure(lattice, Dict(pid=>[0, 0, 1] for pid in lattice.pids))
    lswt = Algorithm(:FM, LSWT(lattice, hilbert, (J,), ms))

    path = ReciprocalPath(lattice.reciprocals, rectangle"Γ-X-M-Γ", length=8)
    data = lswt(:EBS, EnergyBands(path))[2].data[2]
    A(; k) = 2-cos(k[1])-cos(k[2])
    for (i, params) in enumerate(pairs(path))
        @test isapprox(A(; params...), data[i, 1], atol=atol)
        @test isapprox(A(; params...), data[i, 2], atol=atol)
    end

    path = ReciprocalPath(lattice.reciprocals, rectangle"Γ-X-M-Γ", length=100)
    ebs = lswt(:EBS, EnergyBands(path))
    ins = lswt(:INS, InelasticNeutronSpectra(path, range(0.0, 5.0, length=501); η=0.3, log=true))
    plt = plot()
    plot!(plt, ins)
    plot!(plt, ebs, color=:white, linestyle=:dash)
    display(plt)
    savefig("spectra.png")
end


text:
using PtFEM

data = Dict(
  # GenericSolid(ndim, nst, nels, nn, nip, finite_element(nod, nodof), axisymmentrix)
  :struc_el => GenericSolid(3, 6, 6, 8, 1, Tetrahedron(4, 3), false),
  :properties => [100.0 0.3 0.0;],
  :g_coord => [
    0.0 0.0 0.0; 1.0 0.0 0.0; 0.0 0.0 -1.0; 1.0 0.0 -1.0;
    0.0 1.0 0.0; 1.0 1.0 0.0; 0.0 1.0 -1.0; 1.0 1.0 -1.0],
  :g_num => [
    1, 3, 4, 7,
    1, 4, 2, 7,
    1, 2, 5, 7,
    6, 4, 8, 7,
    6, 2, 4, 7,
    6, 5, 2, 7
  ],
  :support => [
    (1, [0 0 1]), (2, [1 0 1]), (3, [0 0 0]), (4, [1 0 0]),
    (5, [0 1 1]), (7, [0 1 0]), (8, [1 1 0])
  ],
  :loaded_nodes => [
    ( 1, [0.0 0.0 -0.1667]),
    ( 2, [0.0 0.0 -0.3333]),
    ( 5, [0.0 0.0 -0.3333]),
    ( 6, [0.0 0.0 -0.1667])
  ]
)

data |> display
println()

@time m, dis_df, fm_df = p54(data)
println()

dis_df |> display
println()

println("\nThe integration point (nip = $(data[:struc_el].nip)) stresses are:\n")

fm_df |> display
println()


text:
SXType = UInt8 # stores the type of SEXPREC

##############################################################################
##
## A table of R's SEXPREC tags from RInternals.h
##
##############################################################################

const NILSXP =      0x00 # nil = NULL
const SYMSXP =      0x01 # symbols
const LISTSXP =     0x02 # lists of dotted pairs
const CLOSXP =      0x03 # closures
const ENVSXP =      0x04 # environments
const PROMSXP =     0x05 # promises: [un]evaluated closure arguments
const LANGSXP =     0x06 # language constructs (special lists)
const SPECIALSXP =  0x07 # special forms
const BUILTINSXP =  0x08 # builtin non-special forms
const CHARSXP =     0x09 # "scalar" string type (internal only)
const LGLSXP =      0x0A # logical vectors
# 11 and 12 were factors and ordered factors in the 1990s
const INTSXP =      0x0D # integer vectors
const REALSXP =     0x0E # real variables
const CPLXSXP =     0x0F # complex variables
const STRSXP =      0x10 # string vectors
const DOTSXP =      0x11 # dot-dot-dot object
const ANYSXP = 	    0x12 # make "any" args work. Used in specifying types for symbol registration to mean anything is okay
const VECSXP =      0x13 # generic vectors
const EXPRSXP =     0x14 # expressions vectors
const BCODESXP =    0x15 # byte code
const EXTPTRSXP =   0x16 # external pointer
const WEAKREFSXP =  0x17 # weak reference
const RAWSXP =      0x18 # raw bytes
const S4SXP =       0x19 # S4, non-vector

# used for detecting PROTECT issues in memory.c
const NEWSXP =      0x1E # fresh node creaed in new page
const FREESXP =     0x1F # node released by GC

const FUNSXP =      0x63 # Closure or Builtin or Special

#=
 = Administrative SXP values
 =#
const REFSXP =            0xFF
const NILVALUE_SXP =      0xFE
const GLOBALENV_SXP =     0xFD
const UNBOUNDVALUE_SXP =  0xFC
const MISSINGARG_SXP =    0xFB
const BASENAMESPACE_SXP = 0xFA
const NAMESPACESXP =      0xF9
const PACKAGESXP =        0xF8
const PERSISTSXP =        0xF7
const CLASSREFSXP =       0xF6
const GENERICREFSXP =     0xF5
const BCREPDEF =          0xF4
const BCREPREF =          0xF3
const EMPTYENV_SXP =      0xF2
const BASEENV_SXP =       0xF1
const ATTRLANGSXP =       0xF0
const ATTRLISTSXP =       0xEF

##############################################################################
##
## R objects in the file are preceded by a UInt32 giving the type and
## some flags.  These functions unpack bits in the flags.  The isobj
## bit might be useful for distinguishing an RInteger from a factor or
## an RList from a data.frame.
##
##############################################################################

typealias RDATag UInt32

isobj(fl::RDATag) = (fl & 0x00000100) != 0
hasattr(fl::RDATag) = (fl & 0x00000200) != 0
hastag(fl::RDATag) = (fl & 0x00000400) != 0
sxtype(fl::RDATag) = fl % UInt8

##############################################################################
##
## Julia types encapsulating various types of R SEXPRECs
##
##############################################################################

"""
    Base class for RData internal representation of all R types.
    `SEXPREC` stands for S (R predecessor) expression record.
"""
abstract RSEXPREC{S}

"""
    R symbol.
    Not quite the same as a Julia symbol.
"""
immutable RSymbol <: RSEXPREC{SYMSXP}
    displayname::RString
end

"""
    Base class for all R types (objects) that can have attributes.
"""
abstract ROBJ{S} <: RSEXPREC{S}

"""
   Base class for all R vector-like objects.
"""
abstract RVEC{T, S} <: ROBJ{S}

"""
    R vector object.
"""
type RVector{T, S} <: RVEC{T, S}
    data::Vector{T}
    attr::Hash                   # collection of R object attributes

    RVector(v::Vector{T} = T[], attr::Hash = Hash()) = new(v, attr)
end

typealias RLogicalVector RVector{Int32, LGLSXP}
typealias RIntegerVector RVector{Int32, INTSXP}
typealias RNumericVector RVector{Float64, REALSXP}
typealias RComplexVector RVector{Complex128, CPLXSXP}

"""
    R vector object with explicit NA values.
"""
immutable RNullableVector{T, S} <: RVEC{T, S}
    data::Vector{T}
    na::BitVector                # mask of NA elements
    attr::Hash                   # collection of R object attributes
end

typealias RStringVector RNullableVector{RString,STRSXP}
typealias RList RVector{RSEXPREC,VECSXP}  # "list" in R == Julia cell array

"""
    Representation of R's paired list-like structures (`LISTSXP`, `LANGSXP`).
    Unlike R that represents it as singly-linked list,
    `RPairList` uses vector representation.
"""
immutable RPairList <: ROBJ{LISTSXP}
    items::Vector{RSEXPREC}
    tags::Vector{RString}
    attr::Hash

    RPairList(attr::Hash = Hash()) = new(RSEXPREC[], RString[], attr)
end

function Base.push!(pl::RPairList, item::RSEXPREC, tag::RString)
    push!(pl.tags, tag)
    push!(pl.items, item)
end

type RClosure <: ROBJ{CLOSXP}
    formals
    body
    env
    attr::Hash

    RClosure(attr::Hash = Hash()) = new(nothing, nothing, nothing, attr)
end

immutable RBuiltin <: RSEXPREC{BUILTINSXP}
    internal_function::RString
end

type RPromise <: ROBJ{PROMSXP}
    value
    expr
    env
    attr::Hash

    RPromise(attr::Hash = Hash()) = new(nothing, nothing, nothing, attr)
end

type REnvironment <: ROBJ{ENVSXP}
    enclosed
    frame
    hashtab
    attr::Hash

    REnvironment() = new(nothing, nothing, nothing, Hash())
end

immutable RRaw <: ROBJ{RAWSXP}
    data::Vector{UInt8}
    attr::Hash
end

immutable RS4Obj <: ROBJ{S4SXP}
    attr::Hash
end

type RExtPtr <: ROBJ{EXTPTRSXP}
    protected
    tag
    attr::Hash

    RExtPtr() = new(nothing, nothing, Hash())
end

type RBytecode <: ROBJ{BCODESXP}
    attr::Hash
    tag
    car
    cdr
    RBytecode(code = nothing, consts = nothing, attr::Hash = Hash()) =
        new(attr, nothing, code, consts)
end

immutable RPackage <: RSEXPREC{PACKAGESXP}
    name::Vector{RString}
end

immutable RNamespace <: RSEXPREC{NAMESPACESXP}
    name::Vector{RString}
end

# R objects without body (empty environments, missing args etc)
immutable RDummy{S} <: RSEXPREC{S}
end

##############################################################################
##
## Utilities for working with basic properties of R objects:
##    attributes, class inheritance, etc
##
##############################################################################

const emptystrvec = RString[]

hasattr(ro::ROBJ, attrnm) = haskey(ro.attr, attrnm)
hasnames(ro::ROBJ) = hasattr(ro, "names")
hasdim(ro::ROBJ) = hasattr(ro, "dim")
hasdimnames(ro::ROBJ) = hasattr(ro, "dimnames")
getattr(ro::ROBJ, attrnm) = getindex(ro.attr, attrnm).data
getattr(ro::ROBJ, attrnm, default) = hasattr(ro, attrnm) ? getindex(ro.attr, attrnm).data : default

Base.names(ro::ROBJ) = getattr(ro, "names")

class(ro::ROBJ) = getattr(ro, "class", emptystrvec)
class(x) = emptystrvec
inherits(x, clnm) = any(class(x) .== clnm)

isdataframe(rl::RList) = inherits(rl, "data.frame")
isfactor(ri::RIntegerVector) = inherits(ri, "factor")

Base.length(rl::RVEC) = length(rl.data)
Base.size(rv::RVEC) = length(rv.data)
Base.size(rl::RList) = isdataframe(rl) ? (length(rl.data[1]), length(rl.data)) : length(rl.data)

row_names(ro::ROBJ) = getattr(ro, "row.names", emptystrvec)


text:
ENV["GKS_WSTYPE"] = "mov"  # or "webm", "mp4"
ENV["GKS_VIDEO_OPTS"] = "600x450@25@2x"

using GR

x = LinRange(0, 2π, 200)
for i = 0:100
    plot(x, sin.(x .+ i / 10.0))
end


text:
module RulesBernoulliMarginalsTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_marginalrules

@testset "marginalrules:Bernoulli" begin

     @testset "out_p: (m_out::PointMass, m_p::Beta)" begin

        @test_marginalrules [ with_float_conversions = true ] Bernoulli(:out_p) [
            (input = (m_out = PointMass(1.0), m_p = Beta(2.0, 1.0)), output = (out = PointMass(1.0), p = Beta(3.0, 1.0))),
            (input = (m_out = PointMass(1.0), m_p = Beta(4.0, 2.0)), output = (out = PointMass(1.0), p = Beta(5.0, 2.0))),
            (input = (m_out = PointMass(0.0), m_p = Beta(1.0, 2.0)), output = (out = PointMass(0.0), p = Beta(1.0, 3.0))),

        ]

    end
end
end


text:
"""
    benchmark(pkg, baseline, candidate="HEAD") -> Tuple{BenchmarkGroup,BenchmarkGroup}

Runs the package benchmarks for the Git revisions `baseline` and `candidate` and returns the
results. Package are expected to have a "bench/benchmarks.jl" file which returns the
benchmark suite.
"""
function benchmark(pkg::AbstractString, baseline::AbstractString, candidate::AbstractString="HEAD")
    repo = LibGit2.GitRepo(Pkg.dir(pkg))
    org_head = LibGit2.head_oid(repo)

    info("Benchmarking baseline ($baseline)")
    checkout_safe!(repo, baseline)
    println(string(LibGit2.head_oid(repo)))  # Show SHA
    trial_baseline = benchmark(pkg)

    info("Benchmarking candidate ($candidate)")
    checkout_safe!(repo, candidate)
    println(string(LibGit2.head_oid(repo)))  # Show SHA
    trial_candidate = benchmark(pkg)

    checkout_safe!(repo, org_head)
    return trial_candidate, trial_baseline
end

function benchmark(pkg::AbstractString)
    results_file = tempname()
    touch(results_file)
    code = """
        using BenchmarkTools
        open("$(escape_string(results_file))", "w") do f
            suite = include(Pkg.dir("$(escape_string(pkg))", "bench", "benchmarks.jl"))
            results = run(suite, verbose=true)
            serialize(f, results)
        end
    """
    results = try
        # TODO: We may want to allow re-building the package between runs
        run(`$(Base.julia_cmd()) --compilecache=no --eval $code`)
        open(results_file, "r") do f
            deserialize(f)
        end
    finally
        isfile(results_file) && rm(results_file)
    end

    return results
end


text:
## Nearfield dynamics
#==========================================================================================#

kutta_joukowsky(ρ, Γ, V, l) = ρ * Γ * V × l
trailing_velocity(r, horseshoe :: Horseshoe, Γ, V) = trailing_legs_velocities(r1(r, horseshoe), r2(r, horseshoe), Γ, V)

"""
    midpoint_velocity(r, Ω, horseshoes, Γs, U)

Evaluate the induced velocity by the trailing legs at a given location ``r``, by summing over the velocities of Horseshoes with vortex strengths ``\\Gamma``s, rotation rates ``\\Omega``, and a freestream flow vector ``U`` in the aircraft reference frame.
"""
midpoint_velocity(r, horseshoes, Γs, U, Ω) = sum(x -> trailing_velocity(r, x[1], x[2], -normalize(U)), zip(horseshoes, Γs)) - U - Ω × r

nearfield_drag(force, U) = -dot(force, normalize(U))
horseshoe_moment(horseshoe :: Horseshoe, force, r_ref) = (bound_leg_center(horseshoe) - r_ref) × force
moments(horseshoes, forces, r_ref) = horseshoe_moment.(horseshoes, forces, Ref(r_ref))

"""
    nearfield_forces(Γ_comp, hs_comp, Γs, horseshoes, U, Ω, ρ)

Compute the nearfield forces via the local Kutta-Jowkowski theorem given an array of horseshoes `hs_comp` to compute the forces on a component, their associated vortex strengths `Γ_comp`, the arrays of horseshoes and vortex strengths `Γs`  of the entire aircraft, the freestream flow vector ``U``, rotation rates ``\\Omega``, and a density ``\\rho``. The velocities are evaluated at the midpoint of the bound leg of each horseshoe, excluding the contribution of the bound leg.
"""
nearfield_forces(Γ_focus, hs_focus, Γs, horseshoes, U, Ω, ρ) = kutta_joukowsky.(Ref(ρ), Γ_focus, midpoint_velocity.(bound_leg_center.(hs_focus), Ref(horseshoes), Ref(Γs), Ref(U), Ref(Ω)), bound_leg_vector.(hs_focus))

function nearfield_dynamics(Γ_focus, hs_focus, Γs, horseshoes, U, Ω, ρ, r_ref)
    geom_forces  = nearfield_forces(Γ_focus, hs_focus, Γs, horseshoes, U, Ω, ρ)
    geom_moments = moments(hs_focus, geom_forces, r_ref)

    geom_forces, geom_moments
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
using Groonga
using Base.Test

Groonga.init()
context = Groonga.ctx_open()

@test Groonga.default_encoding() == "utf8"

@test_throws ErrorException Groonga.db_create(context, "a/b/c/d/e")
@test_throws ErrorException Groonga.db_open(context, "a/b/c/d/e")

dir = mktempdir()
database = Groonga.db_create(context, "$dir/test.db")
println(Groonga.execute_command(context, "status"))

Groonga.execute_command(context, "table_create --name Site --flags TABLE_HASH_KEY --key_type ShortText")
println(Groonga.execute_command(context, "select --table Site"))

Groonga.execute_command(context, "column_create --table Site --name title --type ShortText")
println(Groonga.execute_command(context, "select --table Site"))

Groonga.execute_command(context, "load --table Site [{\"_key\":\"http://example.org/\",\"title\":\"This is test record 1!\"}, {\"_key\":\"http://example.net/\",\"title\":\"test record 2.\"}]")
println(Groonga.execute_command(context, "select --table Site --query _id:1 --output_columns _key,title,_score"))

Groonga.fin()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "OpenrestyBuilder"
version = v"1.15.8"

# Collection of sources required to build OpenrestyBuilder
sources = [
    "https://openresty.org/download/openresty-1.15.8.2.tar.gz" =>
    "436b4e84d547a97a18cf7a2522daf819da8087b188468946b5a89c0dd1ca5d16",

    "https://ftp.pcre.org/pub/pcre/pcre-8.43.tar.gz" =>
    "0b8e7465dc5e98c757cc3650a20a7843ee4c3edf50aaf60bb33fd879690d2c73",

    "https://www.openssl.org/source/openssl-1.0.2t.tar.gz" =>
    "14cb464efe7ac6b54799b34456bd69558a749a4931ecfd9cf9f71d7881cac7bc",

    "https://www.zlib.net/zlib-1.2.11.tar.gz" =>
    "c3e5e9fdd5004dcb542feda5ee4f0ff0744628baf8ed2dd5d66f8ca1197cb1a1",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd openresty-1.15.8.2/
./configure --prefix=$prefix --with-cc=$CC --with-zlib=$WORKSPACE/srcdir/zlib-1.2.11 --with-openssl=$WORKSPACE/srcdir/openssl-1.0.2t --with-pcre=$WORKSPACE/srcdir/pcre-8.43 --with-pcre-jit
make install
rm $prefix/bin/openresty
cd $prefix/bin/
ln -fs ../nginx/sbin/nginx ./openresty
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:musl),
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "openresty", :openresty)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: tanmaykm/OpenrestyBuilder

max_issues_repo_name: tanmaykm/OpenrestyBuilder

max_forks_repo_name: tanmaykm/OpenrestyBuilder


text:
using PointDefects
using Test

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


text:
struct Order{N} end
Order(N) = Order{N}()

function get_process(::Order{1}, T)
    1
end

function get_process(::Order{2}, T)
    @SMatrix [
        1 T
        0 1
    ]
end

function get_process(::Order{3}, T)
    @SMatrix [
        1 T T^2 / 2
        0 1 T
        0 0 1
    ]
end

function get_process_covariance(::Order{1}, T)
    @SMatrix [T] # SMatrix more efficient than scalar for cholesky
end

function get_process_covariance(::Order{2}, T)
    @SMatrix [
        T^3/3 T^2/2
        T^2/2 T
    ]
end

function get_process_covariance(::Order{3}, T)
    @SMatrix [
        T^5/20 T^4/8 T^3/6
        T^4/8 T^3/3 T^2/2
        T^3/6 T^2/2 T
    ]
end

function soft_bound(value, noise, upper_bound, lower_bound, relative_soft_bounding = 0.2)
    value += value[1] > upper_bound * (1 - sign(upper_bound) * relative_soft_bounding) && value[end] > 0 && noise[end] > 0 ?
        SVector(noise[1], -noise[end]) : noise
    value += value[1] < lower_bound * (1 + sign(lower_bound) * relative_soft_bounding) && value[end] < 0 && noise[end] < 0 ?
        SVector(noise[1], -noise[end]) : noise
    value
end


text:
using IterativeSolvers
using FactCheck

include("getDivGrad.jl")

facts("cg") do

context("Small full system") do
    N=10
    A = randn(N,N)
    A = A'*A
    rhs = randn(N)
    tol = 1e-12
    x,ch = cg(A,rhs;tol=tol, maxiter=2*N)

    @fact norm(A*x - rhs) --> less_than(cond(A)*√tol)
    @fact ch.isconverged --> true

    # If you start from the exact solution, you should converge immediately
    x2,ch2 = cg!(A\rhs, A, rhs; tol=tol*10)
    @fact length(ch2[:resnorm]) --> less_than_or_equal(1)

    # Test with cholfact should converge immediately
    F = cholfact(A)
    x2,ch2 = cg(A, rhs, F)
    @fact length(ch2[:resnorm]) --> less_than_or_equal(2)
end

context("Sparse Laplacian") do
    A = getDivGrad(32,32,32)
    L = tril(A)
    D = diag(A)
    U = triu(A)
    JAC(x) = D.\x
    SGS(x) = L\(D.*(U\x))

    rhs = randn(size(A,2))
    rhs/= norm(rhs)
    tol = 1e-5

    context("matrix") do
        xCG, = cg(A,rhs;tol=tol,maxiter=100)
        xJAC, = cg(A,rhs,JAC;tol=tol,maxiter=100)
        xSGS, = cg(A,rhs,SGS;tol=tol,maxiter=100)
        @fact norm(A*xCG - rhs) --> less_than_or_equal(tol)
        @fact norm(A*xSGS - rhs) --> less_than_or_equal(tol)
        @fact norm(A*xJAC - rhs) --> less_than_or_equal(tol)
    end

    Af = MatrixFcn(A)
    context("function") do
        xCG, = cg(Af,rhs;tol=tol,maxiter=100)
        xJAC, = cg(Af,rhs,JAC;tol=tol,maxiter=100)
        xSGS, = cg(Af,rhs,SGS;tol=tol,maxiter=100)
        @fact norm(A*xCG - rhs) --> less_than_or_equal(tol)
        @fact norm(A*xSGS - rhs) --> less_than_or_equal(tol)
        @fact norm(A*xJAC - rhs) --> less_than_or_equal(tol)
    end

    context("function with specified starting guess") do
        tol = 1e-4
        x0 = randn(size(rhs))
        xCG, hCG = cg!(copy(x0),Af,rhs,identity;tol=tol,maxiter=100)
        xJAC, hJAC = cg!(copy(x0),Af,rhs,JAC;tol=tol,maxiter=100)
        xSGS, hSGS = cg!(copy(x0),Af,rhs,SGS;tol=tol,maxiter=100)
        @fact norm(A*xCG - rhs) --> less_than_or_equal(tol)
        @fact norm(A*xSGS - rhs) --> less_than_or_equal(tol)
        @fact norm(A*xJAC - rhs) --> less_than_or_equal(tol)

        iterCG = length(hCG[:resnorm])
        iterJAC = length(hJAC[:resnorm])
        iterSGS = length(hSGS[:resnorm])
        @fact iterJAC --> iterCG
        @fact iterSGS --> less_than_or_equal(iterJAC) "Preconditioner increased the number of iterations"
    end
end

end


text:
module BlockMatching

using OffsetArrays
using Distances
using Requires

export FullSearch, best_match, multi_match

abstract type AbstractBlockMatchingStrategy end

"""
    best_match(S::AbstractBlockMatchingStrategy, ref, [frame=ref]; offset=false)
    best_match(S::AbstractBlockMatchingStrategy, ref, [frame=ref], p; offset=false)

For given pixel `p` in `frame`, find the best match pixel `q` in `ref` with block matching strategy
`S`. If `p` is not provided, it operates on all pixels `p ∈ CartesianIndices(ref)` and returns an array.

!!! note
    For clarification, "pixel" here means the position `CartesianIndex` instead of the value `Colorant`.

# Arguments

* `S::AbstractBlockMatchingStrategy`: required
  The concrete block matching strategy and relevant configs, e.g., patch size, search window size.
  See `subtypes(BlockMatching.AbstractBlockMatchingStrategy)` for possible strategies.
* `ref`::AbstractArray: required
  Reference image where the matched pixel `q` belongs to.
* `frame::AbstractArray`: optional
  The current image/frame where the input pixel `p` belongs to. If `frame` is not provided, it will
  be `ref`.
* `p::CartesianIndex`: optional
  The given pixel `p` that block matching operates on. If not provided, block matching is operating
  on the whole image and thus returns an array of pixels `q`s.

# Parameters

* `offset::Bool`: optional
  if `true`, it returns the motion vector `q-p` instead of the pixel `q`. The default value is
  `false`.

# Output

If `p` is given, it returns `CartesianIndex{N}`. Otherwise the block matching operates on
the whole image and returns `Array{CartesianIndex{N}, N}`, i.e., the motion field.

# Examples

If pixel `p` is provided, block matching only operates on the given pixel, the output is the matched
pixels `q`s for `p`.

```jldoctest best_match
using Images, TestImages
using BlockMatching

ref = imresize(testimage("cameraman"), (64, 64))
img = imrotate(ref, 0.2, CartesianIndices(ref).indices)

S = FullSearch(SqEuclidean(), ndims(img), patch_radius=5, search_radius=11)
p = CartesianIndex(17, 17)
best_match(S, img, ref, p)

# output

CartesianIndex(17, 14)
```

If `p` is not provided, the block matching operates on the whole image. Each item of the output array
is the matched pixels `q`s for `p`, i.e., `best_match(S, ref, frame)[p] = best_match(S, ref, frame, p)`.

```jldoctest best_match
julia> matches = best_match(S, img, ref);

julia> summary(matches)
"64×64 $(Matrix{CartesianIndex{2}})"

julia> matches[p] == best_match(S, img, ref, p)
true
```

!!! tip

    `best_match(S, ref, frame)` is usually more performant than `map(p->best_match(S, ref, frame,
    p), CartesianIndices(frame))` because intermediate results can be cached/pre-allocated so as
    to reduce unnecessary computation and memory allocation.

# See also

* [`multi_match`](@ref) can be used when multiple candidates are wanted.
"""
best_match(S::AbstractBlockMatchingStrategy, ref; kwargs...) = best_match(S, ref, ref; kwargs...)
best_match(S::AbstractBlockMatchingStrategy, ref, p::CartesianIndex; kwargs...) = best_match(S, ref, ref, p; kwargs...)

"""
    multi_match(S::AbstractBlockMatchingStrategy, ref, [frame=ref]; num_patches, offset=true)
    multi_match(S::AbstractBlockMatchingStrategy, ref, [frame=ref], p::CartesianIndex; num_patches, offset=true)

For given pixel `p` in `frame`, find matched pixels `q`s in `ref` with block matching strategy `S`.

# Arguments

* `S::AbstractBlockMatchingStrategy`: required
  The concrete block matching strategy and relevant configs, e.g., similarity mesaure, patch size,
  search window size. See `subtypes(BlockMatching.AbstractBlockMatchingStrategy)` for possible
  strategies.
* `ref`::AbstractArray: required
  Reference image where the matched pixel `q` belongs to.
* `frame::AbstractArray`: optional
  The current image/frame where the input pixel `p` belongs to. If `frame` is not provided, it will
  be `ref`.
* `p::CartesianIndex`: optional
  The given pixel `p` that block matching operates on. If not provided, block matching is operating
  on the whole image and thus returns an array of pixels `q`s.

# Parameters

* `num_patches`: required
  The number of matched pixels. If the candidates number is smaller than `num_patches`, the behavior
  is undefined; it may or may not throw errors.
* `offset::Bool`: optional
  if `true`, it returns the motion vector `q-p` instead of the pixel `q`. The default value is
  `false`.

# Output

If `p` is given, it returns `Vector{CartesianIndex{N}}`. Otherwise the block matching operates on
the whole image and returns `Array{CartesianIndex{N}, N+1}`, where the extra N+1 dimension stores
the matched vector.

# Examples

If pixel `p` is provided, block matching only operates on the given pixel, the output is the matched
pixels `q`s for `p`.

```jldoctest multi_match
using Images, TestImages
using BlockMatching

ref = imresize(testimage("cameraman"), (64, 64))
img = imrotate(ref, 0.2, CartesianIndices(ref).indices)

S = FullSearch(SqEuclidean(), ndims(img), patch_radius=5, search_radius=11)
p = CartesianIndex(17, 17)
multi_match(S, img, ref, p; num_patches=2)

# output

2-element $(Vector{CartesianIndex{2}}):
 CartesianIndex(17, 14)
 CartesianIndex(18, 13)
```

Otherwise, the block matching operates on the whole image. Each item of the output array is the
matched pixels `q`s for `p`, i.e., `multi_match(S, ref, frame)[p] = multi_match(S, ref, frame, p)`.

```jldoctest multi_match
julia> matches = multi_match(S, img, ref; num_patches=2);

julia> summary(matches)
"64×64 $(Matrix{Vector{CartesianIndex{2}}})"

julia> matches[p, :] == multi_match(S, img, ref, p; num_patches=2)
true
```

!!! tip

    `multi_match(S, ref, frame)` is usually more performant than `map(p->multi_match(S, ref, frame,
    p), CartesianIndices(frame))` because intermediate results can be cached so as to reduce
    computation and memory allocation.

# See also

[`best_match`](@ref) is more performant a choice when `num_patches==1`. It also supports more
sophisticated strategies, e.g., ThreeStepSearch.
"""
multi_match(S::AbstractBlockMatchingStrategy, ref; kwargs...) = multi_match(S, ref, ref; kwargs...)
multi_match(S::AbstractBlockMatchingStrategy, ref, p::CartesianIndex; kwargs...) = multi_match(S, ref, ref, p; kwargs...)

include("utils.jl")
include("strategies/full_search.jl")

function __init__()
    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
        include("strategies/full_search_cuda.jl")
    end
end

end


text:
module GeneticScreens

using MatrixLM
import MatrixLM.calc_preds, MatrixLM.calc_preds!,
    MatrixLM.calc_resid, MatrixLM.calc_resid!

using DataFrames
using Loess
using Statistics
import StatsBase.iqr
using CSV

export Response, Predictors, RawData, get_X, get_Z, get_Y,
    add_intercept, remove_intercept, shuffle_rows, shuffle_cols, contr,
    Mlm, mlm, t_stat, perm_pvals, mlm_perms,
    get_powers, get_dose_slopes, read_plate,
    backest_sum!, mlm_backest_sum, mlm_backest_sum_perms,
    coef, predict, fitted, resid,
    S_score, S_score_perms

# Get powers and cross products (for spatial effects)
include("get_powers.jl")

# Dosage slopes
include("get_dose_slopes.jl")

# Read in plate data and convert to RawData object
include("read_plate.jl")

# MLM that back-estimates sum contrasts
include("mlm_backest_sum.jl")
# Permutations
include("mlm_backest_sum_perms.jl")
# Predictions and residuals
include("predict.jl")

# Collins S score and permutations
include("S_score.jl")
include("S_score_perms.jl")

end


text:
function constraint_voltage_product_converter(pm::AbstractWRMModel, wr, wi, w_fr, w_to)
    # # @constraint(pm.model, norm([ 2*wr;  2*wi; w_fr-w_to ]) <= w_fr+w_to )
    # @constraint(pm.model, [ 2*wr,  2*wi, w_fr-w_to, w_fr+w_to] in JuMP.RotatedSecondOrderCone())
    @constraint(pm.model, [w_fr/sqrt(2), w_to/sqrt(2), wr, wi] in JuMP.RotatedSecondOrderCone())
end

"""
Links converter power & current

```
pconv_ac[i]^2 + pconv_dc[i]^2 <= wc[i] * iconv_ac_sq[i]
pconv_ac[i]^2 + pconv_dc[i]^2 <= (Umax)^2 * (iconv_ac[i])^2
```
"""
function constraint_converter_current(pm::AbstractWRMModel, n::Int, cnd::Int, i::Int, Umax, Imax)
    wc = PowerModels.var(pm, n, cnd, :wc_ac, i)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qconv_ac = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)
    iconv_sq = PowerModels.var(pm, n, cnd, :iconv_ac_sq, i)

    # PowerModels.con(pm, n, cnd, :conv_i)[i]  = @constraint(pm.model, norm([2pconv_ac;2qconv_ac; (wc-iconv_sq)]) <= (wc+iconv_sq))
    # PowerModels.con(pm, n, cnd, :conv_i_sqrt)[i] = @constraint(pm.model, norm([2*pconv_ac; 2*qconv_ac]) <= 2*Umax*iconv)
    # @constraint(pm.model, norm([2*iconv; iconv_sq-1]) <= iconv_sq + 1)

    PowerModels.con(pm, n, cnd, :conv_i)[i]  = @constraint(pm.model, [wc/sqrt(2), iconv_sq/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    PowerModels.con(pm, n, cnd, :conv_i_sqrt)[i] = @constraint(pm.model, [Umax * iconv/sqrt(2), Umax * iconv/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    # @constraint(pm.model, [iconv_sq/(2*sqrt(2)), iconv_sq/(2*sqrt(2)), iconv/sqrt(2), iconv/sqrt(2)] in JuMP.RotatedSecondOrderCone())
    @constraint(pm.model, iconv_sq <= iconv*Imax)
end


text:
module Curry

using InteractiveUtils

export
curry,
🍛,
partial,
swap

"""
  swap(f)

Return a function with a reordered argument list.  Reordered so first argument becomes the last.
"""
swap(f::Function) = (a, b...)->f(b..., a)

"""
  partial(f, x...)

Takes a function f and fewer than the normal arguments to f, and
returns a function that takes a variable number of additional args. When
called, the returned function calls f with args + additional args.
"""
partial = (f::Function, y...)->(z...)->f(y..., z...)

∂(f::Function, y...) = partial(f, y...)

function curryn(f, n, args)
  if n == 1
    x->f(args...,x)
  else
    x->curryn(f, n-1, vcat(args,x))
  end
end

curryn(f, n) = curryn(f,n,[])

function curryn(f, n, args; kwargs...)
  if n == 1
    x->f(args..., x; kwargs...)
  else
    x->curryn(f, n-1, vcat(args, x); kwargs...)
  end
end

curryn(f, n; kwargs...) = curryn(f,n,[]; kwargs...)

function curry2(f)
  x->y->f(x, y)
end

function curry3(f)
  x->y->z->f(x, y, z)
end

function curry4(f)
  r->x->y->z->f(r, x, y, z)
end

"""
   curry(f)

Convert a function that accepts multiple arguments into an equivalent chain of single argument functions.
Functions are chained as single argument anonymous functions that return single argument functions.

If multiple arities of a function exist, the one with the least number of arguments greater than 1 will be curried.

# Examples
```jldoctest
julia> curry(map)(sqrt)(1:4)
4-element Array{Float64,1}:
 1.0
 1.4142135623730951
 1.7320508075688772
 2.0

# arities 1 and greater than 2 exist for -, but arity 2 will be choosen for currying
julia> curry((-))(9)
#119 (generic function with 1 method)

julia> curry((-))(9)(1)
8
```

"""
@generated function curry(f)
  m = methodswith(f)
  minargs = mapreduce(x->x.nargs, min, m)-1
  maxargs = mapreduce(x->x.nargs, max, m)-1
  if maxargs > 1
    n = max(2, minargs)
    if n == 2
      :(curry2(f))
    elseif n == 3
      :(curry3(f))
    elseif n == 4
      :(curry4(f))
    else
      :(curryn(f, $n))
    end
  else
    :(f)
  end
end

@generated function curry(f; kwargs...)
  m = methodswith(f)
  minargs = mapreduce(x->x.nargs, min, m)-1
  maxargs = mapreduce(x->x.nargs, max, m)-1
  n = max(2, minargs)
  :(curryn(f, $n; kwargs...))
end

🍛(f) = curry(f)
🍛(f; kwargs...) = curry(f; kwargs...)

end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mc01td(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    DPMAX = 10
    P = Array{Float64,1}(undef, DPMAX+1)
    DWORK = Array{Float64,1}(undef, 2*DPMAX+2)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    DP = parse(BlasInt, vs[1])
    DICO =  vs[2][1]
    if ( DP<=-1 || DP>DPMAX )
        @error "Illegal DP=$DP"
    end

    DPP = DP

    vs = String[]
    _isz = DP+1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    P[1:_isz] .= parsex.(Float64, vs)
    close(f)

# interp call 1

    DP, STABLE, NZ, INFO, IWARN = SLICOT.mc01td!(DICO, DP, P)
    @test INFO == 0
    INFO == 0 || return
    println(io, "DP = $DP")
    println(io, "STABLE = $STABLE")
    println(io, "NZ = $NZ")

    if ( IWARN!=0 )
        @warn "mc01td returned IWARN=$IWARN"
        return
    end # if
    stable = Bool(STABLE)
    if ( stable )
        println(io, "P is stable")
    else
        println(io, "P is unstable")
        println(io, "unstable zeros: $NZ")
    end
end # run_mc01td()


text:
using Random
using Printf

include("Algencan/AlgencanWrapper.jl")
include("sh.jl")

import Voronoi
import Covering

drawing = (haskey(ENV, "DRAW") && ENV["DRAW"] != "" && ENV["DRAW"] != "0")
drawingAll = (haskey(ENV, "DRAWALL") && ENV["DRAWALL"] != "" && ENV["DRAWALL"] != "0")

TRIES = 100
WIDTH = 1.0
HEIGHT = 1.0
INSTANCE = 1

As = [
      Array{Covering.Section, 1}([
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.1, 0.1),
                                                                                (0.9, 0.1),
                                                                                (0.9, 0.9),
                                                                                (0.1, 0.9),
                                                                               ])),
                                 ]),
      Array{Covering.Section, 1}([
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.5, 0.5),
                                                                                (0.4, 0.4),
                                                                                (0.5, 0.1),
                                                                                (0.6, 0.4),
                                                                               ])),
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.5, 0.5),
                                                                                (0.6, 0.4),
                                                                                (0.9, 0.5),
                                                                                (0.6, 0.6),
                                                                               ])),
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.5, 0.5),
                                                                                (0.6, 0.6),
                                                                                (0.5, 0.9),
                                                                                (0.4, 0.6),
                                                                               ])),
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.5, 0.5),
                                                                                (0.4, 0.6),
                                                                                (0.1, 0.5),
                                                                                (0.4, 0.4),
                                                                               ])),
                                 ]),
      Array{Covering.Section, 1}([
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0, 0.125),
                                                                                (0.75, 0.5),
                                                                                (1 ,0.875),
                                                                                (0.25, 0.5),
                                                                               ])),
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0, 1),
                                                                                (0, 0.5),
                                                                                (0.5, 1),
                                                                               ]))
                                 ]),
      Array{Covering.Section, 1}([
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.1, 0.1),
                                                                                (0.2, 0.2),
                                                                                (0.2, 0.8),
                                                                                (0.1, 0.9),
                                                                               ])),
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.1, 0.1),
                                                                                (0.9, 0.1),
                                                                                (0.8, 0.2),
                                                                                (0.2, 0.2),
                                                                               ])),
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.9, 0.1),
                                                                                (0.9, 0.9),
                                                                                (0.8, 0.8),
                                                                                (0.8, 0.2),
                                                                               ])),
                                  Covering.Section(Array{Tuple{Real, Real}, 1}([
                                                                                (0.9, 0.9),
                                                                                (0.1, 0.9),
                                                                                (0.2, 0.8),
                                                                                (0.8, 0.8),
                                                                               ])),
                                 ]),
     ]

if haskey(ENV, "INSTANCE")
    INSTANCE = parse(Int, ENV["INSTANCE"])
end
A = As[INSTANCE]

function calculateArea()
    total_area = 0
    for Aⱼ in A
        t, _, _ = Covering.areaAndGradient(Aⱼ, 0)
        total_area += t
    end
    return total_area
end
total_area = calculateArea()

function draw(WIDTH, HEIGHT, r, points, n)
   if !drawing
       return
   end

   if length(points) == 0
       return
   end

   V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
   Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

   Draw.init(WIDTH, HEIGHT)

   for Aⱼ in A
       W = SH.intersect(V, Aⱼ, r)
       #Draw.coveringPartition(W, "xkcd:pink")
       Draw.coveringPartition(W, nothing)
       Draw.filledPolygon(Aⱼ, "blue")
   end

   Draw.circles(r, points)
   Draw.commit()
   Draw.savefig("bateria/" * string(INSTANCE) * "/" * string(n))
   println("...")
end

function coverWithCircles(n, WIDTH, HEIGHT, r, points)
   function pack(r, points)
       x = zeros(2n + 1)
       x[1] = r
       for (i, p) in enumerate(points)
           x[2i] = p[1]
           x[2i + 1] = p[2]
       end

       return x
   end

   function unpack(x)
       r = x[1]

       points = Array{Tuple{Real, Real}}([])

       for i in 1:n
           append!(points, [(x[2i], x[2i + 1])])
       end

       return r, points
   end

   function draw(r, points) if !drawing
           return
       end

       V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
       Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

       Draw.init(WIDTH, HEIGHT)

       for Aⱼ in A
           W = SH.intersect(V, Aⱼ, r)
           Draw.coveringPartition(W, "xkcd:pink")
           Draw.filledPolygon(Aⱼ, "xkcd:black")
       end

       Draw.circles(r, points)
       Draw.commit()
       print("")
   end

   function f(x)
       r, points = unpack(x)

       if drawing && drawingAll
           draw(r, points)
       end

       #println("R: ", r)
       return r
   end

   function ∇f(x)
       grad = zeros(2n + 1)
       grad[1] = 1
       #println("GRAD: ", grad)
       return grad
   end

   function ∇∇f(x)
       return [], [], []
   end

   function c(ind, x)
       r, points = unpack(x)

       #points, repeats = avoidRepeats(points)

       V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
       Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

       uncovered_area = total_area
       for Aⱼ in A
           W = SH.intersect(V, Aⱼ, r)
           area, _, _ = Covering.areaAndGradient(W)
           uncovered_area -= area
       end

       #println("C: ", uncovered_area)
       return uncovered_area
   end

   function ∇c(ind, x)
       r, points = unpack(x)

       #points, repeats = avoidRepeats(points)

       V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
       Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

       gᵣ = 0
       gₛ = fill((0.0, 0.0), n)
       for Aⱼ in A
           W = SH.intersect(V, Aⱼ, r)
           _, gᵣⱼ, gₛⱼ = Covering.areaAndGradient(W)
           gᵣ += gᵣⱼ
           for i in 1:n
               gₛ[i] = (gₛ[i][1] + gₛⱼ[i][1], gₛ[i][2] + gₛⱼ[i][2])
           end
       end

       #for i in repeats
       #    push!(gₛ, (0, 0))
       #end
       #println("GRAD C.: ", pack(gᵣ, gₛ))
       #println("Aperte Enter pra continuar")
       #readline(stdin)

       return collect(1:2n+1), pack(gᵣ, gₛ)
   end

   l = zeros(2n + 1)
   l[1] = 1e-14

   u = zeros(2n + 1)
   u[1] = 10e20
   for i in 2:2n+1
       u[i] = 1.0
   end

   #draw(r, points)
   #println("Aperte Enter pra continuar")
   #readline(stdin)

   try
       x, fx = AlgencanWrapper.optimize(n = 2n + 1, m = 1,

                                        f = f, g = ∇f,# h = ∇∇f,
                                        equatn = [1],
                                        c = c, jac = ∇c,

                                        x = pack(r, points),
                                        l = l,
                                        u = u,

                                        nvparam = 2,
                                        vparam = [
                                                  "ITERATIONS-OUTPUT-DETAIL 10",
                                                  "PENALTY-PARAMETER-INITIAL-VALUE 1000"
                                                 ],
                                        #checkder = 1,
                                        epsopt = 1.0e-6,
                                        epsfeas = 1.0e-6
                                        )

       r, points = unpack(x)
       #draw(r, points)

       println("\n\n=================\nRESULTADO\n=================\n")
       println("r = ", r)
       println("s = ", 1/r)
       println("pontos: ")
       for p in points
           println("(", p[1], ", ", p[2], ")")
       end
       #if drawing
       #    println("Digite ENTER para sair.")
       #    readline(stdin)
       #end

       return r, points
   catch exc
       println("Exception: ", exc)
       return 1e10, []
   end
end

function randf(start, finish, n)
    v = rand(n)
    return map(x -> start + x*(finish-start), v)
end

if drawing
    import Draw
end

function main()
    n = tryparse(Int64, ARGS[1])
    WIDTH = 1.0
    HEIGHT = 1.0

    rmin = 1e10
    pointsmin = []

    #total_area, _, _ = Covering.areaAndGradient(A, 0)
    ϵ = 10e-4

    Random.seed!(1)

    for k in 1:TRIES
        #println("Tentativa ", k)
        points = convert(Array{Tuple{Real, Real}}, collect(zip(randf(1, WIDTH-1, n), randf(1, HEIGHT-1, n))))

        function insideHalfSpace(point, p, q)
            pq = Voronoi.Geometry.subVector(q, p)
            v = Voronoi.Geometry.rotateVectorCCW(pq)

            function dot(a, b)
                return a[1]*b[1] + a[2]*b[2]
            end

            u = Voronoi.Geometry.subVector(point, p)

            return (dot(u, v) >= 0)
        end

        function pointIsInsideConvexPolygon(point, polygon)
            p = polygon.borderHead

            first = true
            while p != polygon.borderHead || first
                first = false
                q = p.next
                if !insideHalfSpace(point, p.origin, q.origin)
                    return false
                end

                p = p.next
            end

            return true
        end

        function pointIsInsideSomeConvexPolygon(point, polygons)
            for polygon in polygons
                if pointIsInsideConvexPolygon(point, polygon)
                    return true
                end
            end

            return false
        end

        for (i, p) in enumerate(points)
            while !pointIsInsideSomeConvexPolygon(points[i], A)
                #println("Ponto $i ainda está fora. Re-sorteando.")
                points[i] = (randf(1, WIDTH-1, 1)[1], randf(1, HEIGHT-1, 1)[1])
            end
        end

        r = (0.5 + (rand(1)[1]))/n

        try
            r, points = coverWithCircles(n, WIDTH, HEIGHT, r, points)
            #@info "points: $(points)"
        catch exc
            println("Exception: ", exc)
            return 1e10, []
        end

        #draw(WIDTH, HEIGHT, r, points, n)

        @info "n = $(n); k = $(k); rmin = $(rmin); r = $(r)"
        if r < rmin
           V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
           Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

           uncovered_area = total_area
           for Aⱼ in A
               W = SH.intersect(V, Aⱼ, r)
               area, _, _ = Covering.areaAndGradient(W)
               uncovered_area -= area
           end

           if uncovered_area < ϵ
               rmin = r
               pointsmin = points
           else
               @warn "UNFEASIBLE"
           end
        end
    end

    println("\n\n=================\nRESULTADO (n = " * string(n) * ")\n=================\n")
    println("r = ", rmin)
    println("s = ", 1/rmin)
    println("pontos: ")
    for p in pointsmin
        println("(", p[1], ", ", p[2], ")")
    end
    println("")

    if drawing
        draw(WIDTH, HEIGHT, rmin, pointsmin, n)
        #readline(stdin)
    end
end

main()

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
module FairnessInprocessing

# Write your package code here.
include("algorithms/algorithms.jl")

end


text:
module Percolation

export dfspercolate, dfspercolate!, colorpercolate, hkpercolate

using DataStructures

"""
Find out whether or not percolation is possible in a 2D grid using depth-first search

Percolation is assumed to occur from top to bottom.

# Argument

`grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)

# Returns

boolean indicating whether percolation is possible
"""
function dfspercolate(grid::BitMatrix)
    height, width = size(grid)

    undiscovered = BitMatrix(undef, height + 2, width + 2)
    # "pad" the array to avoid bound checking for neigboring
    undiscovered[:, 1] .= false
    undiscovered[1, :] .= false
    undiscovered[end, :] .= false
    undiscovered[:, end] .= false
    # insert grid into undiscovered array
    undiscovered[2:end-1, 2:end-1] = grid

    stack = Stack{Tuple{Integer, Integer}}()
    for root in 2:width+1
        if undiscovered[2, root] # top row
            # search this path
            push!(stack, (2, root))
            if dfs!(height + 1, stack, undiscovered)
                return true
            end
        end
    end
    return false
end

"""
Depth-first search algorithm adopted for the dfspercolate function
"""
@inline function dfs!(destination::Integer, stack::Stack, undiscovered::BitMatrix)
    while !isempty(stack)
        vertex = pop!(stack)
        undiscovered[vertex...] = false

        # check if destination is reached
        if vertex[1] == destination
            return true
        end

        # check neighboring vertices
        for step in ((-1, 0), (0, -1), (0, 1), (1, 0))
            neighbor = vertex .+ step
            if undiscovered[neighbor...]
                push!(stack, neighbor)
            end
        end
    end
    return false
end

"""
Find out whether or not percolation is possible in a 2D grid using depth-first search

Percolation is assumed to occur from top to bottom. this version marks points visited by
the dfs in the visited

# Arguments

- `grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)
- `visited::BitMatrix`: whether each point is visited by the dfs

# Returns

boolean indicating whether percolation is possible
"""
function dfspercolate!(grid::BitMatrix, visited::BitMatrix)
    height, width = size(grid)

    undiscovered = BitMatrix(undef, height + 2, width + 2)
    # "pad" the array to avoid bound checking for neigboring
    undiscovered[:, 1] .= false
    undiscovered[1, :] .= false
    undiscovered[end, :] .= false
    undiscovered[:, end] .= false
    # insert grid into undiscovered array
    undiscovered[2:end-1, 2:end-1] = grid

    stack = Stack{Tuple{Integer, Integer}}()
    for root in 2:width+1
        if undiscovered[2, root] # top row
            # search this path
            push!(stack, (2, root))
            if dfs!(height + 1, stack, undiscovered, visited)
                return true
            end
        end
    end
    return false
end

"""
Depth-first search algorithm adopted for the dfspercolate function

this version modifies the `visited` matrix
"""
@inline function dfs!(destination::Integer, stack::Stack, undiscovered::BitMatrix,
        visited::BitMatrix)
    while !isempty(stack)
        vertex = pop!(stack)
        undiscovered[vertex...] = false
        visited[(vertex .- (1, 1))...] = true

        # check if destination is reached
        if vertex[1] == destination
            return true
        end

        # check neighboring vertices
        for step in ((-1, 0), (0, -1), (0, 1), (1, 0))
            neighbor = vertex .+ step
            if undiscovered[neighbor...]
                push!(stack, neighbor)
            end
        end
    end
    return false
end

"""
labels the clusters in a 2D grid for solving the percolation problem

Uses and inefficient O((length * width)^2) coloring algorithm.
Percolation is assumed to occur from top to bottom.

# Argument

`grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)

# Returns

a boolean indicating whether percolation is possible and a matrix "colored" with integers.
Points with value 1 are connected to the starting row.
"""
function colorpercolate(grid::BitMatrix)
    height, width = size(grid)
    colors = zeros(Integer, height + 1, width)
    colors[1, :] .= 1
    counter = 2

    for row in 2:height+1
        for col in 1:width
            if grid[row - 1, col]
                neighbors = neighborhood(colors, row, col)
                nonzero = filter(n -> colors[n...] != 0, neighbors)
                if length(nonzero) == 0
                    colors[row, col] = counter
                    counter += 1
                elseif length(nonzero) == 1
                    colors[row, col] = colors[nonzero[1]...]
                else
                    minneighbor = minimum([colors[n...] for n in nonzero])
                    colors[row, col] = minneighbor
                    for neighbor in nonzero
                        if colors[neighbor...] != minneighbor
                            colors[colors .== colors[neighbor...]] .= minneighbor
                        end
                    end
                end
            end
        end
    end
    return any(x -> x == 1, colors[end, :]), colors[2:end, :]
end

"""
Retruns a vector of tuples of the neighbors of the point at [row, col]
"""
@inline function neighborhood(grid::Matrix, row::Integer, col::Integer)
    neighbors = [(row + 1, col), (row - 1, col), (row, col + 1), (row, col - 1)]
    return filter(n -> checkbounds(Bool, grid, n...), neighbors)
end

"""
Uses the Hoshen-Kopelman algorithm to label the clusters in a 2D grid for percolation

Percolation is assumed to occur from top to bottom.

# Arguments

- `grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)
- `color::Bool` (default = `false`): whether to color the grid after finishing the algorithm

# Returns

[a boolean indicating whether percolation is possible or, if qinfty is true,
the sum of the sizes of the open connecting clusters (integer)],
an integer matrix where clusters, are marked by unique integers, and the cluster lengths.
"""
function hkpercolate(grid::BitMatrix; color::Bool=false, qinfty::Bool=false)
    height, width = size(grid)
    labels = zeros(Integer, height, width)
    list = Vector{Integer}()
    clusters = Vector{Integer}()

    # first row, special case because of top boundary
    hkfirstrow!(grid, labels, list, clusters)

    for row in 2:height
        # first column, special case because of left boundary
        if grid[row, 1]
            hkfirstcol!(labels, list, clusters, row)
        end

        # "middle" points, with appropriate neighbors
        for col in 2:width
            if grid[row, col]
                hkiter!(labels, list, clusters, row, col)
            end
        end
    end

    if color
        hkcolor!(labels, list)
    end
    if qinfty
        return hkresult(labels, list, clusters), labels, clusters
    else
        return hkbool(labels, list), labels, clusters
    end
end

"""
Helper function for `hkpercolate`; Adds new label.
"""
@inline function hknew!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer, col::Integer)
    labels[row, col] = length(list) + 1
    push!(list, labels[row, col])
    push!(clusters, 0)
end

"""
Helper function for `hkpercolate`; Applies the Hoshen-Kopelman algorithm to the first row.
"""
@inline function hkfirstrow!(grid::BitMatrix, labels::Matrix{<:Integer},
        list::Vector{<:Integer}, clusters::Vector{<:Integer})
    width = size(labels)[2]
    if grid[1, 1]
        hknew!(labels, list, clusters, 1, 1)
        clusters[labels[1, 1]] += 1
    end
    for col in 2:width
        if grid[1, col]
            if labels[1, col - 1] != 0
                labels[1, col] = hkfind!(list, labels[1, col - 1])
            else
                hknew!(labels, list, clusters, 1, col)
            end
            clusters[labels[1, col]] += 1
        end
    end
end

"""
Helper function for `hkpercolate`; Applies the Hoshen-Kopelman algorithm to the first column.
"""
@inline function hkfirstcol!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer)
    if labels[row - 1, 1] != 0
        labels[row, 1] = hkfind!(list, labels[row - 1, 1])
    else
        hknew!(labels, list, clusters, row, 1)
    end
    clusters[labels[row, 1]] += 1
end

"""
Helper function for `hkpercolate`; Applies one iteration of the Hoshen-Kopelman algorithm.
"""
@inline function hkiter!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer, col::Integer)
    above, left = labels[row - 1, col], labels[row, col - 1]
    if above == 0
        if left == 0
            hknew!(labels, list, clusters, row, col)
        else
            labels[row, col] = hkfind!(list, left)
        end
    elseif left == 0
        labels[row, col] = hkfind!(list, above)

    else
        labels[row, col] = hkunion!(list, clusters, above, left)
    end
    clusters[labels[row, col]] += 1
end

"""
Union function of the union-find algorithm used by the Hoshen-Kopelman algorithm
"""
@inline function hkunion!(list::Vector{<:Integer}, clusters::Vector{<:Integer},
        x::Integer, y::Integer)
    root = hkfind!(list, x)
    dest = hkfind!(list, y)
    if dest != root
        list[dest] = root
        clusters[root] += clusters[dest]
        clusters[dest] = 0
    end
    return root
end

"""
Find function of the union-find algorithm used by the Hoshen-Kopelman algorithm
"""
@inline function hkfind!(list::Vector{<:Integer}, x::Integer)
    y = x;
    while list[y] != y
        y = list[y];
    end
    while list[x] != x
        z = list[x];
        list[x] = y;
        x = z;
    end
    return y;
end

"""
Same as hkfind!, but doesn't update labels
"""
@inline function hkfind(list::Vector{<:Integer}, x::Integer)
    while x != list[x]
        x = list[x];
    end
    return x
end

"""
Color resultant grid of the Hoshen-Kopelman algorithm appropriately
"""
@inline function hkcolor!(labels::Matrix{<:Integer}, list::Vector{<:Integer})
    width, height = size(labels)
    for row in 1:height
        for col in 1:width
            if labels[row, col] != 0
                labels[row, col] = hkfind(list, labels[row, col])
            end
        end
    end
end

"""
Find the sum of the sizes of the open connecting clusters (integer)
"""
@inline function hkresult(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer})
    sumsize = 0
    # to avoid over-counting the clusters, save the label of each cluster that is counted
    counted = []

    for root in labels[1, :]
        if root > 0
            for dest in labels[end, :]
                if dest > 0
                    toplabel = hkfind(list, root)
                    bottomlabel = hkfind(list, dest)
                    if toplabel == bottomlabel && !(toplabel in counted)
                        sumsize += clusters[toplabel]
                        append!(counted, toplabel)
                    end
                end
            end
        end
    end
    return sumsize
end

"""
Find out whether percolation is possible or not in hkpercolate output
"""
@inline function hkbool(labels::Matrix{<:Integer}, list::Vector{<:Integer})
    for root in labels[1, :]
        if root > 0
            for dest in labels[end, :]
                if dest > 0 && hkfind(list, dest) == hkfind(list, root)
                    return true
                end
            end
        end
    end
    return false
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys