text:
# Collect Time To First X (TTFX)
#
# A few notes about these compile times benchmarks.
#   1. These benchmarks are meant to show where the biggest problems are and to be able to trace back where some regression was introduced.
#   2. The benchmarks use `@eval` to avoid missing compile time, see the `@time` docstring for more info.
#   3. Only add benchmarks for methods which take more than 1 seconds on the first run to reduce noise.
#   4. Note that some benchmarks depend on disk and network speeds too, so focus on the number of allocations since those are more robust.

module Foo end

using UUIDs

# setup required for run_exporession:
const test_notebook_id = uuid1()
let
    channel = Channel{Any}(10)
    Pluto.PlutoRunner.setup_plutologger(
        test_notebook_id,
        channel;
        make_global=false
    )
end
@timeit TOUT "PlutoRunner.run_expression" @eval Pluto.PlutoRunner.run_expression(Foo, Expr(:toplevel, :(1 + 1)), test_notebook_id, uuid1(), nothing);

function wait_for_ready(notebook::Pluto.Notebook)
    while notebook.process_status != Pluto.ProcessStatus.ready
        sleep(0.1)
    end
end

🍭 = Pluto.ServerSession()
🍭.options.server.disable_writing_notebook_files = true
🍭.options.evaluation.workspace_use_distributed = false

path = joinpath(pkgdir(Pluto), "sample", "Basic.jl")

@timeit TOUT "SessionActions.open" nb = @eval Pluto.SessionActions.open(🍭, path; run_async=false)

wait_for_ready(nb)

Pluto.SessionActions.shutdown(🍭, nb; async=false)

# Compile HTTP get.
HTTP.get("http://github.com")

@timeit TOUT "Pluto.run" server_task = @eval let
    port = 13435
    options = Pluto.Configuration.from_flat_kwargs(; port, launch_browser=false, workspace_use_distributed=false, require_secret_for_access=false, require_secret_for_open_links=false)
    🍭 = Pluto.ServerSession(; options)
    server_task = @async Pluto.run(🍭)

    # Give the async task time to start.
    sleep(1)

    HTTP.get("http://localhost:$port/edit").status == 200
    server_task
end


text:
doc"""
    key_sampling(g)
Sample `min(nv(g), 64)` keys in `[1,nv(g)]`
Returns a vector of randomly sampled keys
"""
function key_sampling(
  g::Graph{T}
  ) where T<:Integer

  n_v = nv(g)
  keys = Set{T}()
  sizehint!(keys, 64)
  NBFS = n_v > 64 ? 64 : n_v
  length_keys = 0
  i = 0  #Exit loop after a fixed number of iterations

  while length_keys < NBFS && i < 1000
      keys_required = (NBFS - length_keys)
      rand_keys = sample(1:n_v,keys_required)
      union!(keys, filter(x->degree(g,x)>0, rand_keys))
      length_keys = length(keys)
      i = i + 1;
  end

  return keys
end

doc"""
    kernel_2(g, s)
Perform bfs on graph `g` starting at vertex `s`
Returns vector containing parent info for each vertex
"""
kernel_2(g::Graph{T}, s::T) where T<:Integer = bfs_parents(g,s)

ext: jl


text:
module Generate

bootfile = """
"""

end


text:
@def componentfileds begin
    name :: Symbol
    nodes :: Tuple{Vararg{Int}}
    value :: U
end

abstract type Component{U<:Union{Expr,Number}} end

struct Resistor{U} <: Component{U} @componentfileds end
resistor(name, nodes, value, parameters_string) = Resistor(name, nodes, value)

struct Inductor{U} <: Component{U} @componentfileds end
inductor(name, nodes, value, parameters_string) = Inductor(name, nodes, value)

struct Capacitor{U} <: Component{U} @componentfileds end
capacitor(name, nodes, value, parameters_string) = Capacitor(name, nodes, value)

struct VoltageSource{U} <: Component{U} @componentfileds end
voltageSource(name, nodes, value, parameters_string) = VoltageSource(name, nodes, value)

struct CurrentSource{U} <: Component{U} @componentfileds end
currentSource(name, nodes, value, parameters_string) = CurrentSource(name, nodes, value)

const nodedict_type = Dict{Symbol,Int}

mutable struct ParsedCircuit{N<:Number}
    titleline :: String
    nodedict :: nodedict_type
    max_node :: Int
    max_element :: Int
    netlist :: Vector{Component}
    ParsedCircuit{N}() where N = new("",Dict(Symbol(0)=>0), 0, 0, Vector{Component}())
end

struct MNA{N<:Number, T<:AbstractArray{N,2}}
    G :: T
    H :: T
    g :: Vector{Function}
    D :: T
    S :: T
    s :: Dict{Int,Function}
    nodedict :: nodedict_type
    currentdict :: nodedict_type #cannot use spase vector here, zero(::Function) is undefined
    function MNA{N,T}(G,H,g,D,S,s,nodedict,currentdict) where {N,T}
        @assert begin (G1,G2) = size(G); G1==G2 end "G not square"
        @assert begin (H1,H2) = size(H); H2==length(g) end "H,g mismatch size(H) $(size(H)), $(length(g))"
        @assert begin (D1,D2) = size(D); D1==D2 end "D not square"
        @assert begin (S1,S2) = size(S); S2==1 end "S must be a columb vector"
        if ~(G1==G2==H1==D1==D2==S1) throw(DimensionMismatch()) end
        new{N,T}(G,H,g,D,S,s,nodedict,currentdict)
    end
end

struct MNAbuilder{N<:Number}
    Gi :: Vector{Int}
    Gj :: Vector{Int}
    Gvalue :: Vector{N}
    g :: Vector{Function}
    Hi :: Vector{Int}
    Hj :: Vector{Int}
    Hvalue :: Vector{N}
    Di :: Vector{Int}
    Dj :: Vector{Int}
    Dvalue :: Vector{N}
    Si :: Vector{N}
    Sj :: Vector{N} # all 1 for column vector
    Svalue :: Vector{N}
    s :: Dict{Int,Function}  #  assume sparse
    nodedict :: nodedict_type
    currentdict :: nodedict_type
    MNAbuilder{N}() where N<:Number =
        new([],[],[],[],[],[],[],[],[],[],[],[],[],Dict(),Dict(),Dict())
end


text:
using WorkMitData
using Test

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


text:
using SurfaceWaves
using Documenter

makedocs(;
    modules=[SurfaceWaves],
    authors="Dylan Mikesell, Boise State University",
    repo="https://github.com/dylanmikesell/SurfaceWaves.jl/blob/{commit}{path}#L{line}",
    sitename="SurfaceWaves.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dylanmikesell.github.io/SurfaceWaves.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/dylanmikesell/SurfaceWaves.jl",
)


text:
function construct_correlated_noisefunc(Γ)
  γ = svd(Γ)
  A = γ.U*Diagonal(sqrt.(γ.S))
  b = Vector{eltype(Γ)}(undef,size(Γ,1))
  dist = function (W,dt,rng)
    if typeof(W.dW) <: AbstractArray
      return A*sqrt.(abs(dt))*wiener_randn(rng,size(W.dW))
    else
      return A*sqrt.(abs(dt))*wiener_randn(rng,typeof(W.dW))
    end
  end
  bridge = function (W,W0,Wh,q,h,u,p,t,rng)
    error("Bridging distribution is unknown. Cannot use adapativity")
  end
  dist,bridge
end
CorrelatedWienerProcess(Γ,t0,W0,Z0=nothing;rng = Xorshifts.Xoroshiro128Plus(rand(UInt64))) = NoiseProcess(t0,W0,Z0,construct_correlated_noisefunc(Γ)...,rswm=RSWM(),rng=rng)

function construct_correlated_noisefunc!(Γ)
  γ = svd(Γ)
  A = γ.U*Diagonal(sqrt.(γ.S))
  b = Vector{eltype(Γ)}(undef,size(Γ,1))
  dist = function (rand_vec,W,dt,u,p,t,rng)
    wiener_randn!(rng,b)
    b .*= sqrt.(abs(dt))
    mul!(rand_vec,A,b)
  end
  bridge = function (rand_vec,W,W0,Wh,q,h,u,p,t,rng)
    error("Bridging distribution is unknown. Cannot use adapativity")
  end
  dist,bridge
end
CorrelatedWienerProcess!(Γ,t0,W0,Z0=nothing;rng = Xorshifts.Xoroshiro128Plus(rand(UInt64))) =
                         NoiseProcess(t0,W0,Z0,construct_correlated_noisefunc!(Γ)...,
                         rswm=RSWM(),rng = rng)


text:
# test hist
@test sum(fit(Histogram,[1,2,3]).weights) == 3
@test fit(Histogram,[]).weights == []
@test fit(Histogram,[1]).weights == [1]
@test fit(Histogram,[1,2,3],[0,2,4]) == Histogram([0,2,4],[2,1])
@test fit(Histogram,[1,2,3],[0,2,4]) != Histogram([0,2,4],[1,1])
@test fit(Histogram,[1,2,3],0:2:4) == Histogram(0:2:4,[2,1])
@test all(fit(Histogram,[1:100]/100,0.0:0.01:1.0).weights .==1)
@test fit(Histogram,[1,1,1,1,1]).weights[1] == 5
@test sum(fit(Histogram,(rand(100),rand(100))).weights) == 100
@test fit(Histogram,1:100,nbins=5,closed=:right).weights == [20,20,20,20,20]
@test fit(Histogram,1:100,nbins=5,closed=:left).weights == [19,20,20,20,20,1]
@test fit(Histogram,0:99,nbins=5,closed=:right).weights == [1,20,20,20,20,19]
@test fit(Histogram,0:99,nbins=5,closed=:left).weights == [20,20,20,20,20]

@test fit(Histogram,(1:100,1:100),nbins=5).weights == diagm([20,20,20,20,20])
@test fit(Histogram,(1:100,1:100),nbins=(5,5)).weights == diagm([20,20,20,20,20])

@test fit(Histogram,1:100,weights(ones(100)),nbins=5).weights == [20,20,20,20,20]
@test fit(Histogram,1:100,weights(2*ones(100)),nbins=5).weights == [40,40,40,40,40]

@test eltype(fit(Histogram,1:100,weights(ones(Int,100)),nbins=5).weights) == Int
@test eltype(fit(Histogram,1:100,weights(ones(Float64,100)),nbins=5).weights) == Float64

import StatsBase.midpoints

@test midpoints(1.0:1.0:10.0) == 1.5:1.0:9.5
@test midpoints(1:10) == 1.5:9.5
@test midpoints(Float64[1.0:1.0:10.0]) == Float64[1.5:1.0:9.5]


text:
function read_and_populate_model(node_file, edge_file, range_file,
    num_layers, num_nodes_in_layer; has_final_output_bias = true)

    biases = create_empty_biases(num_layers, num_nodes_in_layer)
    weights = create_empty_weights(num_layers, num_nodes_in_layer)
    ranges = []

    bias_data = readdlm(node_file)
    weight_data = readdlm(edge_file)
    range_data = readdlm(range_file)

    for row in bias_data[2:end]
        values = map(x -> parse(Float64, x), split(row, ","))
        layer_id = Int(values[1])
        node_id = Int(values[2])
        b = values[3]
        biases[layer_id][node_id] = b
    end

    for row in weight_data[2:end]
        values = map(x -> parse(Float64, x), split(row, ","))
        layer_id = Int(values[1])
        f = Int(values[3])
        t = Int(values[4])
        w = values[5]
        weights[layer_id][EdgeId(f, t)] = w
    end

    for row in range_data[2:end]
        values = map(x -> parse(Float64, x), split(row, ","))
        push!(ranges, (values[1], values[2]))
    end

    input_range = ranges[1:6]
    output_range = ranges[7]

    @info "reading completed"

    m, var  = get_milp(num_layers, num_nodes_in_layer, weights, biases,
        input_range, output_range)

    @info "MILP created"
    return m, var

end

function get_overlap_percentage(args::Dict)::Float64
    node_file = args["node_file"]
    edge_file = args["edge_file"]
    range_file = args["range_file"]
    num_layers = args["num_layers"]
    num_nodes_in_layer = args["num_nodes_in_layer"]
    case_file = args["casefile"]
    pmax_value = args["pmax"]

    m, var = read_and_populate_model(node_file,
        edge_file, range_file,
        num_layers,
        num_nodes_in_layer)

    input = var[:input]
    output = var[:output]

    case_data = readdlm(case_file)

    @constraint(m, input[1] == case_data[1])
    @constraint(m, input[2] == case_data[2])
    @constraint(m, input[3] == case_data[3])
    @constraint(m, input[4] == case_data[4])
    @constraint(m, input[6] == case_data[5])

    JuMP.set_lower_bound(input[5], 0.05)
    JuMP.set_upper_bound(input[5], 0.75)
    JuMP.set_upper_bound(output, 4.5)

    @objective(m, Min, (output - pmax_value) * (output - pmax_value))

    solve_model(m)

    return JuMP.value(input[5])

end


text:
function norm_ast(a::Any)
    if isa(a, Expr)
        for (i, arg) in enumerate(a.args)
            a.args[i] = norm_ast(arg)
        end
        if a.head === :line
            return Expr(:line, a.args[1], :none)
        end
        if a.head === :macrocall
            fa = a.args[1]
            long_enough = length(a.args) >= 3
            if fa === Symbol("@int128_str") && long_enough
                return Base.parse(Int128, a.args[3])
            elseif fa === Symbol("@uint128_str") && long_enough
                return Base.parse(UInt128, a.args[3])
            elseif fa === Symbol("@bigint_str") && long_enough
                return  Base.parse(BigInt, a.args[3])
            elseif fa == Symbol("@big_str") && long_enough
                s = a.args[3]
                n = tryparse(BigInt, s)
                if !(n === nothing)
                    return (n)
                end
                n = tryparse(BigFloat, s)
                if !(n === nothing)
                    return isnan((n)) ? :NaN : (n)
                end
                return s
            end
        elseif length(a.args) >= 2 && Meta.isexpr(a, :call) && a.args[1] == :- && isa(a.args[2], Number)
            return -a.args[2]
        end
        return a
    elseif isa(a, QuoteNode)
        return Expr(:quote, norm_ast(a.value))
    elseif isa(a, AbstractFloat) && isnan(a)
        return :NaN
    end
    return a
end

function meta_parse_file(str)
    pos = 1
    x1 = Expr(:file)
    try
        while pos <= sizeof(str)
            x, pos = Meta.parse(str, pos)
            push!(x1.args, x)
        end
    catch er
        isa(er, InterruptException) && rethrow(er)

        return x1, true
    end
    if length(x1.args) > 0  && x1.args[end] === nothing
        pop!(x1.args)
    end
    x1 = norm_ast(x1)
    CSTParser.remlineinfo!(x1)
    return x1, false
end

function cst_parse_file(str)
    x, ps = CSTParser.parse(CSTParser.ParseState(str), true)
    sp = CSTParser.check_span(x)
    # remove leading/trailing nothings
    if length(x.args) > 0 && CSTParser.is_nothing(x.args[1])
        popfirst!(x.args)
    end

    if !isempty(sp)
        @error "CST spans inconsistent!"
    end

    x0 = norm_ast(to_codeobject(x))
    x0, CSTParser.has_error(ps) || !isempty(sp)
end

@testset "Parsing files in Base" begin
    dir = joinpath(Sys.BINDIR, Base.DATAROOTDIR)
    for (root, _, files) in walkdir(dir)
        for fpath in files
            file = joinpath(root, fpath)
            endswith(file, ".jl") || continue

            str = read(file, String)

            cst = CSTParser.parse(str, true)
            cst_expr, cst_err = cst_parse_file(str)
            meta_expr, meta_err = meta_parse_file(str)
            @test cst_err == meta_err
            @test cst.fullspan == sizeof(str)

            if cst_err || meta_err
                if cst_err && !meta_err
                    @error "CSTParser.parse errored, but Meta.parse didn't." file=file
                elseif !cst_err && meta_err
                    @error "Meta.parse errored, but CSTParser.parse didn't." file=file
                end
            else
                if cst_expr == meta_expr
                    @test true
                else
                    @error "parsing difference" file=file
                    c0, c1 = CSTParser.compare(cst_expr, meta_expr)
                    printstyled(c0, bold = true, color = :light_red)
                    println()
                    printstyled(c1, bold=true, color=:light_green)
                    println()
                    @test false
                end
            end
        end
    end
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module CommQu

using Random
using ExtractMacro
using DataStructures
using ..Interface
using ..Common

using ..DeltaE.ArraySets

export GraphCommQu

import ..Interface: energy, delta_energy, update_cache!, neighbors

function gen_ξ(N::Integer, P::Integer)
    ξ = BitArray(undef, P, N)
    ξv = Vector{BitVector}(undef, P)
    for a = 1:P
        ξv[a] = bitrand(N)
        ξ[a,:] = ξv[a]
    end
    # the ξv representation is for energy
    # the ξ representation is for delta_energy and update_cache!
    y = bitrand(P)
    return ξ, ξv, y
end

mutable struct Stabilities
    K2::Int
    p2::ArraySet
    m2::ArraySet
    Δ1s::Vector{Vector{Int}}
    Δ2s::Vector{Int}
    δΔ2Hsmin::Vector{MutableBinaryMinHeap{Float64}}
    δΔ2Hsmax::Vector{MutableBinaryMaxHeap{Float64}}
    ξsi::BitVector
    last_move::Int
    function Stabilities(K2::Integer, P::Integer)
        p2 = ArraySet(P)
        m2 = ArraySet(P)
        Δ1s = [zeros(P) for k in 1:K2]
        Δ2s = zeros(P)
        δΔ2Hsmin = [MutableBinaryMinHeap{Float64}() for a = 1:P]
        δΔ2Hsmax = [MutableBinaryMaxHeap{Float64}() for a = 1:P]
        ξsi = BitArray(undef, P)
        return new(K2, p2, m2, Δ1s, Δ2s, δΔ2Hsmin, δΔ2Hsmax, ξsi, 0)
    end
end

mutable struct GraphCommQu <: SimpleGraph{Int}
    N::Int
    K1::Int
    K2::Int
    P::Int
    ξ::BitMatrix
    ξv::Vector{BitVector}
    y::BitVector
    stab::Stabilities
    tmps::BitVector
    function GraphCommQu(K2::Int, ξ::BitMatrix, ξv::Vector{BitVector}, y::BitVector)
        P, N = size(ξ)
        @assert N % K2 == 0
        K1 = N ÷ K2
        @assert length(ξv) == P
        @assert all(length(ξ1) == N for ξ1 in ξv)
        @assert length(y) == P
        iseven(K1) || throw(ArgumentError("K1 must be even, given: $K1"))
        iseven(K2) || throw(ArgumentError("K2 must be even, given: $K2"))
        stab = Stabilities(K2, P)
        tmps = BitVector(undef, N)
        return new(N, K1, K2, P, ξ, ξv, y, stab, tmps)
    end
end

"""
    GraphCommQu(K1::Integer, K2::Integer, P::Integer; fc = false) <: SimpleGraph{Int}

A `SimpleGraph` implementing a two-layer binary committee machine with `K2` hidden units,
where each hidden unit has `K1` binary (\$±1\$) synapses, trained on `P` random i.i.d. \$±1\$
patterns. The output of each hidden unit is computed through a quadratic function \$x^2\$,
and then half of the `K2` units have a negative weight associated to their output.
Both `K1` and `K2` must be even.
The default is to generate a tree-like machine, but it can be made fully-connected with the
argument `fc=true`.

The energy of the model is computed as the number of misclassified patterns.
"""
function GraphCommQu(K1::Integer, K2::Integer, P::Integer; fc::Bool = false)
    Kin = fc ? K1 : K1 * K2
    ξ, ξv, y = gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphCommQu(K2, ξ, ξv, y)
end

function Base.empty!(stab::Stabilities)
    @extract stab : p2 m2 Δ1s Δ2s δΔ2Hsmin δΔ2Hsmax ξsi
    empty!(p2)
    empty!(m2)
    foreach(v->fill!(v, 0), Δ1s)
    fill!(Δ2s, 0)
    P = length(δΔ2Hsmin)
    for a = 1:P
        δΔ2Hsmin[a] = MutableBinaryMinHeap{Float64}()
        δΔ2Hsmax[a] = MutableBinaryMaxHeap{Float64}()
    end
    stab.last_move = 0
    return stab
end

function energy(X::GraphCommQu, C::Config)
    @assert X.N == C.N
    @extract C : s
    @extract X : N K1 K2 P ξv y stab tmps
    @extract stab : p2 m2 Δ1s Δ2s δΔ2Hsmin δΔ2Hsmax

    E = 0
    empty!(stab)
    tmps = BitArray(undef, K1)

    for a = 1:P
        o = 2y[a] - 1
        Δ2 = 0
        δΔ2Hmax = MutableBinaryMaxHeap{Float64}()
        δΔ2Hmin = MutableBinaryMinHeap{Float64}()
        for k = 1:K2
            c = o * (2 * (2k ≤ K2) - 1)
            unsafe_copyto!(tmps, 1, s, (k-1)*K1 + 1, K1)
            map!(⊻, tmps, tmps, ξv[a][(k-1)*K1 .+ (1:K1)])
            Δ1 = K1 - 2 * sum(tmps)
            Δ1s[k][a] = Δ1
            Δ2 += c * Δ1^2
            push!(δΔ2Hmax, 4 * (c + abs(c * Δ1)))
            push!(δΔ2Hmin, 4 * (c - abs(c * Δ1)))
        end
        Δ2s[a] = Δ2
        δΔ2Hsmax[a] = δΔ2Hmax
        δΔ2Hsmin[a] = δΔ2Hmin
        if Δ2 > 0
            Δ2 + first(δΔ2Hmin) ≤ 0 && push!(p2, a)
        else
            Δ2 + first(δΔ2Hmax) > 0 && push!(m2, a)
            E += 1
        end
    end

    return E
end

function update_cache!(X::GraphCommQu, C::Config, move::Int)
    @assert X.N == C.N
    @extract C : s
    @extract X : N K1 K2 P ξ y stab
    @extract stab : p2 m2 Δ1s Δ2s δΔ2Hsmin δΔ2Hsmax ξsi last_move
    @assert N == length(s)

    si = s[move]
    k2 = (move - 1) ÷ K1 + 1
    # k1 = (move - 1) % K1 + 1
    last_move ≠ move && unsafe_copyto!(ξsi, 1, ξ, (move-1)*P + 1, P)
    stab.last_move = move
    si && flipbits!(ξsi)
    Δ1k = Δ1s[k2]
    c = 2 * (2k2 ≤ K2) - 1
    # @assert ξsi == ξ[:,move] ⊻ si
    @inbounds for a = 1:P
        co = c * (2y[a] - 1)
        oldΔ2 = Δ2s[a]
        oldΔ1 = Δ1k[a]
        newΔ1 = oldΔ1 + (2 - 4 * ξsi[a])

        newΔ2 = oldΔ2 + co * (newΔ1^2 - oldΔ1^2)

        # oldδΔ2max = first(δΔ2Hsmax[a])
        # oldδΔ2min = first(δΔ2Hsmin[a])

        update!(δΔ2Hsmax[a], k2, 4 * (co + abs(co * newΔ1)))
        update!(δΔ2Hsmin[a], k2, 4 * (co - abs(co * newΔ1)))

        newδΔ2max = first(δΔ2Hsmax[a])
        newδΔ2min = first(δΔ2Hsmin[a])

        oldm2 = a ∈ m2
        newm2 = newΔ2 ≤ 0 && newΔ2 + newδΔ2max > 0
        oldp2 = a ∈ p2
        newp2 = newΔ2 > 0 && newΔ2 + newδΔ2min ≤ 0

        # oldm2 && (@assert oldΔ2 ≤ 0 && oldΔ2 + oldδΔ2max > 0)
        # oldp2 && (@assert oldΔ2 > 0 && oldΔ2 + oldδΔ2min ≤ 0)

        oldm2 && !newm2 && delete!(m2, a)
        !oldm2 && newm2 && push!(m2, a)
        oldp2 && !newp2 && delete!(p2, a)
        !oldp2 && newp2 && push!(p2, a)

        Δ1k[a] = newΔ1
        Δ2s[a] = newΔ2
    end
    si && flipbits!(ξsi)

    return
end

function delta_energy_naive(X::GraphCommQu, C::Config, move::Int)
    @assert C.N == X.N
    @extract C : s
    @extract X : N P ξ stab
    @assert 1 ≤ move ≤ N

    spinflip!(C, move)
    E1 = energy(X, C)
    spinflip!(C, move)
    E0 = energy(X, C)
    return E1 - E0
end

function delta_energy(X::GraphCommQu, C::Config, move::Int)
    @assert C.N == X.N
    @extract C : s
    @extract X : N K1 K2 P ξ y stab
    @extract stab : p2 m2 Δ1s Δ2s ξsi last_move
    @assert 1 ≤ move ≤ N

    si = s[move]
    k2 = (move - 1) ÷ K1 + 1
    # k1 = (move - 1) % K1 + 1
    c = 2 * (2k2 ≤ K2) - 1

    last_move ≠ move && unsafe_copyto!(ξsi, 1, ξ, (move-1)*P + 1, P)
    stab.last_move = move
    si && flipbits!(ξsi)
    Δ1k = Δ1s[k2]
    ΔE = 0
    @inbounds for a in p2
        co = c * (2y[a] - 1)
        if Δ2s[a] + co * 4 * (1 - (1 - 2*ξsi[a]) * Δ1k[a]) ≤ 0
            ΔE += 1
        end
    end
    @inbounds for a in m2
        co = c * (2y[a] - 1)
        if Δ2s[a] + co * 4 * (1 - (1 - 2*ξsi[a]) * Δ1k[a]) > 0
            ΔE -= 1
        end
    end
    si && flipbits!(ξsi)
    # td = delta_energy_naive(X, C, move)
    # ΔE == td || @show ΔE, td
    # @assert ΔE == td
    return ΔE
end

function check_delta(X::GraphCommQu, C::Config, move::Int)
    e0 = energy(X, C)
    delta = delta_energy(X, C, move)
    spinflip!(C, move)
    e1 = energy(X, C)
    spinflip!(C, move)

    (e1-e0) == delta || (@show e1,e0,delta,e1-e0; error())
end

neighbors(X::GraphCommQu, i::Int) = AllButOne(X.N, i)

end


text:
module JavisNB

using Javis

import Javis: AbstractObject, Video

import Interact
import Interact: @map, @layout!, Widget, Widgets, hbox, vbox

export embed

#################
#### Pluto Viewer
#################

"""
    _pluto_viewer(video::Video, frames::Int, actions::Vector)

Creates an interactive viewer in a Pluto Notebook by storing all the frames in-memory.

# Returns

`arr`: List of Javis frames
"""
function _pluto_viewer(video::Video, frames::Int, objects::Vector)
    arr = collect(
        Javis.get_javis_frame(video, objects, frame; layers = video.layers) for
        frame in 1:frames
    )
    return arr
end

"""
    PlutoViewer

Wrapper to assist viewing rendered gifs as cell outputs of Pluto notebooks
when `liveview = false`.
"""
struct PlutoViewer
    filename::String
end

function Base.show(io::IO, ::MIME"image/png", v::PlutoViewer)
    write(io, read(v.filename))
end

###################
#### Jupyter Viewer
###################

"""
    _jupyter_viewer(video::Video, frames::Int, actions::Vector)
Creates an interactive viewer in a Jupyter Notebook.
"""
function _jupyter_viewer(
    video::Video,
    frames::Int,
    objects::Vector{AbstractObject},
    framerate::Int,
)
    t = Interact.textbox(1:frames, value = 1, typ = Int)
    f = Interact.slider(1:frames, label = "Frame", value = t)
    obs = Interact.Observables.throttle(1 / framerate, f)
    output = @map Javis.get_javis_frame(video, objects, &obs; layers = video.layers)
    wdg = Widget(["f" => f, "t" => t], output = output)
    @layout! wdg vbox(hbox(:f, :t), output)
end

##########
#### embed
##########

"""
    embed(pathname::String)

Shows a gif within a notebook by calling the path to the gif.

# Arguments

- `pathname::String`: path of the gif to render in the notebook
"""
function embed(pathname::String)

    if (isdefined(Main, :IJulia) && Main.IJulia.inited) || endswith(@__FILE__, ".ipynb")
        display(MIME("text/html"), """<img src="$(pathname)">""")
    elseif isdefined(Main, :PlutoRunner)
        return PlutoViewer(pathname)
    end

    return pathname
end

function embed(vid::Video, frames::Int, objects::Vector{AbstractObject}, framerate::Int)
    return _jupyter_viewer(vid, frames, objects, framerate)
end

function embed(vid::Video, frames::Int, objects::Vector{AbstractObject})
    return _pluto_viewer(vid, frames, objects)
end

embed(::Nothing) = nothing

"""
    embed(
        video::Video;
        framerate=30,
        pathname="javis_GIBBERISH.gif",
        liveview=false,
        tempdirectory="",
        ffmpeg_loglevel="panic",
        rescale_factor=1.0,
        postprocess_frame=identity,
        postprocess_frames_flow=default_postprocess
    )

The core function of `JavisNB`.
It allows one to show the output of `Javis.render`
from [Javis](https://juliaanimators.github.io/Javis.jl/stable/) directly within a notebook.
Notebooks supported are:

- IJulia
- Pluto

To use `embed` just replace `Javis.render` with `embed`.
Here is an example one can use within a `Pluto` notebook:

```julia
# In one cell
using Javis, JavisNB

# In a new one
begin
    function ground(args...)
        background("white")
        sethue("black")
    end
    vid = Video(500, 500)
    Background(1:50, ground)
    o = Object(JCircle(O, 20, action=:fill))
    act!(o, Action(1:25, anim_translate(Point(20, 0))))
    act!(o, Action(26:50, anim_translate(Point(-20, 0))))

    # In pure Javis here we would have
    # render(args...;kwargs...)
    # Instead, within a notebook one uses
    # embed(args...;kwargs...)
    # For example
    embed(vid, pathname="rendered_with_embed.gif")
end
```
"""
function embed(
    video::Video;
    framerate = 30,
    pathname = "javis_$(randstring(7)).gif",
    liveview = false,
    tempdirectory = "",
    ffmpeg_loglevel = "panic",
    rescale_factor = 1.0,
    postprocess_frames_flow = identity,
    postprocess_frame = Javis.default_postprocess,
)

    out = render(
        video,
        framerate = framerate,
        pathname = pathname,
        liveview = liveview,
        streamconfig = nothing,
        tempdirectory = tempdirectory,
        ffmpeg_loglevel = ffmpeg_loglevel,
        rescale_factor = rescale_factor,
        postprocess_frames_flow = postprocess_frames_flow,
        postprocess_frame = postprocess_frame,
    )
    return if !liveview
        embed(out)
    else
        embed(out...)
    end
end

end


text:
using Documenter, Impute

makedocs(
    modules=[Impute],
    format=:html,
    pages=[
        "Home" => "index.md",
        "Impute" => "api/impute.md",
        "Context" => "api/context.md",
        "Imputors" => "api/imputors.md",
        "Utilities" => "api/utils.md",
    ],
    repo="https://github.com/invenia/Impute.jl/blob/{commit}{path}#L{line}",
    sitename="Impute.jl",
    authors="Invenia Technical Computing Corporation",
    assets=["assets/invenia.css"],
)

deploydocs(
    repo = "github.com/invenia/Impute.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
module Interpolation

export polyint

"""    polyint(xs, ys[, n])

Returns a function that performs polynomial interpolation across the `xs` and
`ys`.

`n`, if given, is the order of the interpolating polynomial; if the length of
the `xs` and `ys` is greater than `n+1`, then multiple interpolants are
returned.
"""
function polyint(xs, ys)
    n = length(xs)
    polyint(xs, ys, n-1)
end

function polyint(xs::Array{T, 1}, ys::Array{T, 1}, n) where T
    triangle = Array{T,1}[ys]
    for i in 2:n+1
        push!(triangle, zeros(T, length(triangle[i-1])-1))
    end

    full =  n == length(xs)-1

    function interpolate(x)
        for i in 2:n+1
            k = length(triangle[i])
            for j in 1:k
                triangle[i][j] = (x-xs[j+i-1])*triangle[i-1][j]/(xs[j]-xs[j+i-1]) + (x-xs[j])*triangle[i-1][j+1]/(xs[j+i-1]-xs[j])
            end
        end

        if full
            triangle[end][1]
        else
            Tuple(triangle[end])
        end
    end
end

end # module


text:
using Distributed # Standard library for distributed memory parallel computing

@everywhere using Ipopt # Library for non-linear optimization of continuous systems

include("code-2-lib/parsers.jl") # Insert the contents of parsers.jl here
include("code-2-lib/lib.jl") # Insert the contents of lib.jl here.

function compute_solution2(con_file::String, inl_file::String, raw_file::String, rop_file::String, time_limit::Int, scoring_method::Int, network_model::String; output_dir::String="", scenario_id::String="none") # Function named "compute_solution2"
    time_data_start = time() # To note the start time
    PowerModels.silence() # Suppresses information and warning messages output by PowerModels
    goc_data = parse_goc_files(con_file, inl_file, raw_file, rop_file, scenario_id=scenario_id) # Parse the input files
    network = build_pm_model(goc_data) # Calls "build_pm_model" function to make model based on input data

    sol = read_solution1(network, output_dir=output_dir) # Reads the base solution
    PowerModels.update_data!(network, sol) # To update network with the values from sol

    check_network_solution(network) # Checks feasibility criteria of network solution, produces an error if a problem is found

    network_tmp = deepcopy(network) # Copies all fields from "network"
    balance = compute_power_balance_deltas!(network_tmp) # Compute power balance deltas for active & reactive power

    if balance.p_delta_abs_max > 0.01 || balance.q_delta_abs_max > 0.01 # Check the delta limit for active & reactive power
        error(LOGGER, "solution1 power balance requirements not satified (all power balance values should be below 0.01). $(balance)")
    end
    load_time = time() - time_data_start # Find difference between the start & end time

    ###### Prepare Solution 2 ######

    time_contingencies_start = time() # To note the contingencies start time

    gen_cont_total = length(network["gen_contingencies"]) # Finds how many generator contingencies
    branch_cont_total = length(network["branch_contingencies"])  # Finds how many branch contingencies
    cont_total = gen_cont_total + branch_cont_total # Finds total contingencies

    cont_order = contingency_order(network) # To retrieve the order of contingency

    workers = Distributed.workers() # Retrieve the workers

    process_data = [] # An empty dic.

    cont_per_proc = cont_total/length(workers) # Number of contingencies per procs

    for p in 1:length(workers)
        cont_start = trunc(Int, ceil(1+(p-1)*cont_per_proc))
        cont_end = min(cont_total, trunc(Int,ceil(p*cont_per_proc)))
        pd = (
            pid = p,
            processes = length(workers),
            con_file = con_file,
            inl_file = inl_file,
            raw_file = raw_file,
            rop_file = rop_file,
            scenario_id = scenario_id,
            output_dir = output_dir,
            cont_range = cont_start:cont_end,
        )
        push!(process_data, pd) # Push pd to process_data dic.
    end

    for (i,pd) in enumerate(process_data)
        info(LOGGER, "worker task $(pd.pid): $(length(pd.cont_range)) / $(pd.cont_range)")
    end

    solution2_files = pmap(solution2_solver, process_data, retry_delays = zeros(3)) # Parallel mapping
    sort!(solution2_files) # In-place sorting: over-write the existing file after sorting
    #println("pmap result: $(solution2_files)")

    time_contingencies = time() - time_contingencies_start # Calculates time difference
    info(LOGGER, "contingency eval time: $(time_contingencies)") # Total evaluation time for contingencies
    info(LOGGER, "time per contingency: $(time_contingencies/cont_total)") # How much evaluation time per contingency

    info(LOGGER, "combine $(length(solution2_files)) solution2 files")
    combine_files(solution2_files, "solution2.txt"; output_dir=output_dir) # Combines the solution files and adds solution2.txt to output directory
    remove_files(solution2_files) # Removes the redundant files

    println("")

    data = [
        "----",
        "scenario id",
        "bus",
        "branch",
        "gen_cont",
        "branch_cont",
        "runtime (sec.)",
    ]
    println(join(data, ", "))

    data = [
        "DATA_SSS",
        goc_data.scenario,
        length(network["bus"]),
        length(network["branch"]),
        length(network["gen_contingencies"]),
        length(network["branch_contingencies"]),
        time_contingencies,
    ]
    println(join(data, ", "))
end

@everywhere function solution2_solver(process_data) # Function "solution2_solver" takes "process_data" as input argument
    #println(process_data)
    time_data_start = time() # To note the start time
    PowerModels.silence() # Suppresses information and warning messages output by PowerModels
    goc_data = parse_goc_files(
        process_data.con_file, process_data.inl_file, process_data.raw_file,
        process_data.rop_file, scenario_id=process_data.scenario_id) # Parsing the input files
    network = build_pm_model(goc_data) # Builds the model based on the input data (goc_data)

    sol = read_solution1(network, output_dir=process_data.output_dir) # Reads the base solution
    PowerModels.update_data!(network, sol) # To update network with the values from sol
    correct_voltage_angles!(network) # Correct the voltage angles. Assumes there is one reference bus and one connected component and adjusts voltage
# angles to be centered around zero at the reference bus.

    time_data = time() - time_data_start # Calculate the time difference between current time and start time

    for (i,bus) in network["bus"]
        if haskey(bus, "evhi")
            bus["vmax"] = bus["evhi"]
        end
        if haskey(bus, "evlo")
            bus["vmin"] = bus["evlo"]
        end
    end

    for (i,branch) in network["branch"]
        if haskey(branch, "rate_c")
            branch["rate_a"] = branch["rate_c"]
        end
    end

    contingencies = contingency_order(network)[process_data.cont_range] # Build a static ordering of all contigencies

    for (i,branch) in network["branch"]
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)
        branch["g"] = g
        branch["b"] = b
        branch["tr"] = tr
        branch["ti"] = ti
    end

    bus_gens = gens_by_bus(network) # Function to check the bus status

    network["delta"] = 0
    for (i,bus) in network["bus"]
        bus["vm_base"] = bus["vm"]
        bus["vm_start"] = bus["vm"]
        bus["va_start"] = bus["va"]
        bus["vm_fixed"] = length(bus_gens[i]) != 0

    end

    for (i,gen) in network["gen"]
        gen["pg_base"] = gen["pg"]
        gen["pg_start"] = gen["pg"]
        gen["qg_start"] = gen["qg"]
        gen["pg_fixed"] = false # To accomodate the contingencies, generator active power is not fixed.
        gen["qg_fixed"] = false # To accomodate the contingencies, generator reactive power is not fixed.
    end

    #nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, mu_init=1e-6, hessian_approximation="limited-memory", print_level=0)
    nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0) # Uses JuMP built-in optimizer
    #nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6)
    #nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, hessian_approximation="limited-memory")

    pad_size = trunc(Int, ceil(log(10,process_data.processes))) # Finding the padding size
    padded_pid = lpad(string(process_data.pid), pad_size, "0") # (lpad - left pad) lpad(string, n, "p") Make a string at least n columns wide when printed, by padding on the left with copies of p.

    solution_filename = "solution2-$(padded_pid).txt" # Assigning the name to the solutin file

    if length(process_data.output_dir) > 0 # Checks the length of output_dir in process_data dic.
        solution_path = joinpath(process_data.output_dir, solution_filename) # If there is content in output_dir then joins path of solutin_filename & output_dir
    else
        solution_path = solution_filename
    end
    if isfile(solution_path) # Checks the existence of solution_path
        warn(LOGGER, "removing existing solution2 file $(solution_path)")
        rm(solution_path)
    end
    open(solution_path, "w") do sol_file # creates an empty file in the case of workers without contingencies

    end

    #network_tmp = deepcopy(network)
    for cont in contingencies
        if cont.type == "gen" # Checks if there is contingency of generator
            info(LOGGER, "working on: $(cont.label)") # Gives the current status of what is going on
            time_start = time() # To note the start time
            network_tmp = deepcopy(network) # Copies all fields from "network"
            debug(LOGGER, "contingency copy time: $(time() - time_start)") # How much time for contingency copy
            network_tmp["cont_label"] = cont.label

            cont_gen = network_tmp["gen"]["$(cont.idx)"]
            cont_gen["contingency"] = true
            cont_gen["gen_status"] = 0 # gen_status during contingency
            pg_lost = cont_gen["pg"]

            time_start = time() # To note the start time
            result = run_fixpoint_pf_v2_2!(network_tmp, pg_lost, ACRPowerModel, nlp_solver, iteration_limit=5) # Core algorithm that computes starting guess for delta and runs initial AC-PF
            debug(LOGGER, "second-stage contingency solve time: $(time() - time_start)") # Finds the time spent for 2nd contingency solution

            result["solution"]["label"] = cont.label # Mentions the contingency label
            result["solution"]["feasible"] = (result["termination_status"] == LOCALLY_SOLVED) # If status is true then solutin is feasible
            result["solution"]["cont_type"] = "gen" #  Mention the contingency type
            result["solution"]["cont_comp_id"] = cont.idx # Mentions the contingency index

            result["solution"]["gen"]["$(cont.idx)"]["pg"] = 0.0 # During contingency, generator active power is zero
            result["solution"]["gen"]["$(cont.idx)"]["qg"] = 0.0 # During contingency, generator reactive power is zero

            correct_contingency_solution!(network, result["solution"]) #  Used for post-processing. Various fallbacks if any step fails
            open(solution_path, "a") do sol_file # Opens the sol_file
                sol2 = write_solution2_contingency(sol_file, network, result["solution"]) # Writes the contingency solutin 2
            end

            network_tmp["gen"]["$(cont.idx)"]["gen_status"] = 1
        elseif cont.type == "branch" # Checks if there is contingency of branch
            info(LOGGER, "working on: $(cont.label)")  # Gives the current status of what is going on
            time_start = time() # To note the start time
            network_tmp = deepcopy(network)  # Copies all fields from "network"
            debug(LOGGER, "contingency copy time: $(time() - time_start)") # How much time for contingency copy
            network_tmp["cont_label"] = cont.label
            network_tmp["branch"]["$(cont.idx)"]["br_status"] = 0 # branch status during contingency

            time_start = time()  # To note the start time
            result = run_fixpoint_pf_v2_2!(network_tmp, 0.0, ACRPowerModel, nlp_solver, iteration_limit=5) # Core algorithm that computes starting guess for delta and runs initial AC-PF
            debug(LOGGER, "second-stage contingency solve time: $(time() - time_start)") # Finds the time spent for 2nd contingency solution

            result["solution"]["label"] = cont.label # Mentions the contingency label
            result["solution"]["feasible"] = (result["termination_status"] == LOCALLY_SOLVED) # If status is true then solutin is feasible
            result["solution"]["cont_type"] = "branch" # Mention the contingency type
            result["solution"]["cont_comp_id"] = cont.idx # Mentions the contingency index

            correct_contingency_solution!(network, result["solution"]) # Used for post-processing. Various fallbacks if any step fails
            open(solution_path, "a") do sol_file # Opens the sol_file
                sol2 = write_solution2_contingency(sol_file, network, result["solution"]) # Writes the contingency solutin 2
            end

            network_tmp["branch"]["$(cont.idx)"]["br_status"] = 1
        else
            @assert("contingency type $(cont.type) not known") #  If contingency other than generator or branch then it is unknown
        end
    end

    return solution_path
end


text:
#= common.jl
Philip Barrett, pobarrett@gmail.com
26may2016, Chicago

Defines polygon type and methods for construction
=#

"""
    type polygon
Defines the polygon type with three entries: the vertices, the normals of the
sides, and the projection of each face onto the normal.
"""
type Polygon
  pts::Matrix{Float64}     # The points representation
  dirs::Matrix{Float64}    # The directions representation
  dists::Vector{Float64}   # The distances associated with the vectors
end

"""
    acwOrder( ptsDirs::Matrix, dists::Vector=[NaN] )
Orients either a set of points, or a a directions/distance pair anit-clockwise.
"""
function acwOrder( ptsDirs::Matrix{Float64}, dists::Vector{Float64}=[NaN] )
    N = size(ptsDirs)[1]
        # Number of points/directions
    if( isnan( dists[1] ) )
        centered = ptsDirs - ones(N) * mean( ptsDirs, 1 )
            # Relative to centre if doing points
            # Enhances numerical stability of orientation
    else
        centered = ptsDirs
    end
    cah = [ centered[i,1] / norm( centered[i,:] ) for i in 1:N ]
        # The cosine
    ord = [ ( ( centered[i,2] >= 0 ) ? cah[i] : - 2 - cah[i] ) for i in 1:N ]
        # Order based on a decreasing transform of the cosine
    ord_idx = sortperm(ord, rev=true)
        # The sorting vector
    ptsDirsOut::Matrix{Float64} = ptsDirs[ ord_idx, : ]
        # Reorder the opints

    if( isnan(dists[1]) )
        return ptsDirsOut
    else
        distsOut::Vector{Float64} = dists[ord_idx]
        return ptsDirsOut, distsOut
    end
end

"""
    deeDoop( pts::Matrix )
De-duplicates a collection of points.  Must already be ordered anti-clockwise.
"""
function deeDoop( pts::Matrix{Float64}, tol::Float64=1e-10 )
  N = size(pts)[1]
      # Number of points
  neighbors = [ pts[ 2:N,: ]; pts[1,:] ]
      # The neighboring points
  diff = [ norm( pts[i,:] - neighbors[i,:] ) for i in 1:N ]
      # The vector of norms
  if all( diff .<= tol )
    return pts[ 1:3, :]
        # The single-point case
  end
  return pts[ diff .> tol, : ]
end

function deeDoop( dirs::Matrix{Float64}, dists::Vector{Float64},
                      tol=1e-8 )
  N = size(dirs)[1]
      # Number of points
  neighbors = [ dirs[ 2:N,: ]; dirs[1,:] ]
      # The neighboring distance vectors
  diff = [ norm( dirs[i,:] - neighbors[i,:] ) for i in 1:N ]
      # The vector of norms
  return dirs[ diff .> tol, : ], dists[ diff .> tol ]
end

function deeDoop( poly::Polygon, tol=1e-10 )
  dirs, dists = deeDoop( poly.dirs, poly.dists )
  return Polygon( dirs=dirs, dists=dists )
end

"""
    isVert( pts::Matrix{Float64} )
Returns a vector of booleans that indicates whether a point is
a vertex
"""
function isVert( pts::Matrix{Float64})
  N = size(pts)[1]
      # Number of points
  neighbors = [ pts[ 2:N,: ]; pts[1,:] ]
      # Neighboring points
  v = pts - neighbors
  v = vcat( v[N,:], v )
      # The vectors along each side
  cosv = [ sum( v[i,:] .* - v[i+1,:] ) /
            ( norm(v[i,:])*norm(v[i+1,:]) ) for i in 1:N ]
      # The cosine of the angle at each point
  keep = [ !isnan(cosv[i]) && cosv[i] > -1 for i in 1:N ]
      # The vertices to retain
  return keep
end

"""
    polygon( ; pts=[ NaN NaN ], dirs=[ NaN NaN ], dists=[NaN]  )
Constructor for poylgon.  Removes duplicates and orders anti-clockwise.
"""
function Polygon( ; pts::Matrix{Float64}=[ NaN NaN ],
                    dirs::Matrix{Float64}=[ NaN NaN ],
                    dists::Vector{Float64}=[NaN]  )

    # Step 1: Order and convert to points if necessary
    if( !isnan( pts[1] ) && isnan( dirs[1] ) )
        pts = acwOrder( pts )
          # If points only
    elseif( !isnan( dirs[1] ) && isnan( pts[1] ) )
        dirs, dists = acwOrder( dirs, dists )
        pts = dirsToPts( dirs, dists )
          # If dists/dirs only
    end
    pts = [ pts[ end, : ]; pts[ 1:(end-1), : ] ]
        # re-initiate the points

    # Step 2: Enforce the convex hull
    pts = chull(pts)

    # Step 3: Remove non-vertex points
    pts = deeDoop( pts, 1e-14 )
        # Strict de-duplication
    pts = pts[ isVert(pts), : ]
        # Remove the non-vertex (and near-duplicate) points
    dirs, dists = ptsToDirs(pts)
        # Create dirs and dists

    # if size(pts)[1] > 3
    #   dirs, dists = deeDoop( dirs, dists )
    #       # The ultimate test of dee-dooping is on the directions.
    #       # This is because we can have many points satisfying the same
    #       # directional constraint (i.e. on a straight line)
    #       # Question: Can I replace this with some test of whether
    #       # points are on more than two lines? (Or very close to?)
    #   pts = dirsToPts( dirs, dists )
    # end
    return Polygon( pts, dirs, dists )
end

"""
    dirsToPts( dirs::Matrix, dists::Vector )
Given a normal/distance description of a set, returns a set of points mZ which
lie at the vertices.  This assumes that the vector of normals is ordered
counter-clockwise already.
"""
function dirsToPts( dirs::Matrix{Float64}, dists::Vector{Float64} )
  nPts = length(dists)
      # Number of points to be computed
  pts = zeros(Float64, 2, nPts)
      # Initialize the (transpose) output vector
  counter = 1
      # Counter: needed in cases where multiple boundary lines intersect at only
      # one point

  # Create the extended matrices (put the first case at the end)
  dirsExt = [ dirs ; dirs[1,:] ]
  distsExt = [ dists; dists[1] ]

  # Main loop: Over all points
  for i in 1:nPts
    A = dirsExt[i:(i+1),:]
        # The directional vectors to be intersected
    b = distsExt[i:(i+1)]
        # And their distances
    sol = A \ b
        # The candidate solution

    # Now compute the error: Detects multiple intersections
    err = norm( A * sol - b )
        # Inversion error
    relerr = ( norm(b) == 0 ) ? err : err / norm(b)
        # Relative error (where possible)

    # If there is no error, then assign to the output
    if( relerr < 1e-14 )
      pts[:,counter] = sol
      counter += 1
    end
  end
  out = transpose(pts)
  return( [ out[ nPts, : ] ; out[ 1:(nPts-1), : ] ] )
end

"""
    ptsToDirs( pts::Matrix )
Given a set of points, returns the normal-distance representation of the edges
Assumes that the points are ordered anti-clockwise already.
"""
function ptsToDirs( pts::Matrix{Float64} )
  nPts = size(pts)[1]
      # Number of points
  flip = [ 0.0 1.0; -1.0 0.0 ]
      # Matrix to convert grdient to normal
  neighbors = [ pts[ 2:nPts, : ] ; pts[1,:] ]
      # The neighboring points
  unscaled = ( pts - neighbors ) * flip
      # The unscaled vaules
  dirs = ( unscaled ./ sqrt( sum( unscaled .* unscaled, 2 ) ) )::Matrix{Float64}
      # The directions
  dists = (vec(sum( pts .* dirs, 2)))::Vector{Float64}
      # The distances in each direction
  return dirs, dists
end

# """
#     ptsCheck( pts::Matrix, dirs::Matrix, dists::Vector )
# Checks that a set of points satisfy all the vector restrictions.  If not,
# returns anly those that do
# """
# function ptsCheck( pts::Matrix, dirs::Matrix, dists::Vector )
#  TODO: This function! :)
# end


text:
function functional_test()
    day_test_functional()
    windowing_test_functional()
    week_test_functional()
    unassigned_test_functional()
    bifurcate_test_functional()

    Logging.info("Scheduler Functional Tests passing. Good job.")
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
using DataStructures
function SimplePeel(A)
    """
    GreedyPeel with parameter p
    """
    n = size(A,1)
    d = vec(round.(Int64,sum(A,dims = 1)))
    numer = sum(d)

    zz = length(findall(x->x==0,d))
    # save graph as adjacency list
    cp = A.colptr
    row = A.rowval
    AdjList = Vector{Vector{Int64}}()
    for i = 1:n
        push!(AdjList,row[cp[i]:cp[i+1]-1])
    end

    D = MutableBinaryMinHeap(d)

    # greedily peel
    ranking = zeros(Int64,n)
    Objs = zeros(n+1)
    Objs[1] = (numer/n)
    next = 2
    denom = n-1
    degrees = zeros(Int64,n)

    # if there are zero values by chance, get rid of them
    for t = 1:zz
        deg, i = top_with_handle(D)
        pop!(D)
        Objs[next] = 0
        denom -= 1
        @assert(deg == 0)
        degrees[next] = 0
        ranking[next-1] = i
        next += 1
    end

    while ~isempty(D)
        deg, i = top_with_handle(D)
        ranking[next-1] = i
        degrees[next-1] = deg
        pop!(D)
        Ni = AdjList[i]
        di = d[i]
        if di > 0
            # Decrease the degree of all neighbors
            d[i] = 0
            for nb = Ni
                removej!(AdjList[nb],i)             # This is no longer a neighboring node, we're removed it
                update!(D,nb,D[nb]-1)  # subtract degree, update the heap
                d[nb] -= 1
            end
        end
        numer -= 2*di
        Objs[next] = numer/denom
        denom -= 1
        next +=1
    end
    Objs = Objs[1:end-1]
    density, w = findmax(Objs)
    S = ranking[w:end]
    degeneracy, w2 = findmax(degrees)
    maxcore = ranking[w2:end]

    return S, density, degeneracy, maxcore, Objs, ranking, degrees
end

function GenPeel(A,p)
    """
    GreedyPeel with parameter p
    """
    if p == 1
        S, density, degeneracy, maxcore, Objs, ranking, degrees = SimplePeel(A)
        return S, density, Objs, ranking
    end

    n = size(A,1)
    d = sum(A,dims = 1)
    numer = sum(d.^p)

    # converting to dictionary helps later
    d2 = Dict{Int64,Int64}()
    for i = 1:n
        d2[i] = d[i]
    end
    d = d2

    # save graph as adjacency list
    cp = A.colptr
    row = A.rowval
    tic = time()
    AdjList = Vector{Vector{Int64}}()
    for i = 1:n
        push!(AdjList,row[cp[i]:cp[i+1]-1])
    end

    # Create initial minheap
    D = zeros(n)
    for i = 1:n
        D[i] = d[i]^p
        for k = AdjList[i]
            D[i] += d[k]^p - (d[k]-1)^p
        end
    end
    H = MutableBinaryMinHeap(D)

    # greedily peel
    ranking = zeros(Int64,n)
    Objs = zeros(n+1)
    Objs[1] = (numer/n)^(1/p)
    next = 2
    denom = n-1
    while ~isempty(H)
        # global numer, denom,next
        v, i = top_with_handle(H)
        # push!(order,i)
        ranking[next-1] = i
        pop!(H)
        Ni = AdjList[i]
        AdjList[i] = Vector{Int64}()
        di = d[i]
        if di > 0
            dichange = di^p - (di-1)^p
            d[i] = 0
            numer -= di^p
            # Remove this from its neighbors list
            for nb = Ni
                removej!(AdjList[nb],i)             # This is no longer a neighboring node, we're removed it
                db = d[nb]
                dbchange = -db^p + (db-1)^p
                Old = H[nb]                         # previous value of D[nb]
                New = Old + dbchange - dichange     # new value now that node i is gone
                update!(H,nb,New)                   # update the heap
                # D[nb] = New
                numer += dbchange
                d[nb] -= 1                          # decrease degree
            end

            # At this point, we have made all changes due to node i disappearing.
            # Now we need to account for changes in D/H
            # that are due to changes in the degree for i's neighbors
            for b = Ni
                db = d[b]
                if db > 0
                    dbchange = 2*db^p - (db-1)^p -(db+1)^p
                    Nb = AdjList[b]
                    for k = Nb
                        update!(H,k,H[k]+dbchange)
                    end
                end
            end

        end

        # if mod(next,100) == 0
        #     println("$next $numer")
        # end
        # if numer < 0
        #     @show numer
        #     break
        # end
        numer = max(numer,0)                # there can be rounding errors at the very end
        Objs[next] = (numer/denom)^(1/p)
        denom -= 1
        next +=1
    end
    Objs = Objs[1:end-1]
    density, w = findmax(Objs)
    S = ranking[w:end]

    return S, density, Objs, ranking

end

# remove element j from vector v, assuming that j appears exactly once
function removej!(v::Vector{Int64},j)
    for i = 1:length(v)
        if v[i] == j
            # v = v[1:i-1,i+1:end]
            deleteat!(v,i)
            return
        end
    end
end

function pdensityobj(A,S,p)

    if length(S) == 0
        return 0
    else
        As = A[S,S]
        d = sum(As,dims = 1)

        return (sum(d.^p)/length(S))^(1/p)
    end

end

max_stars_repo_name: nveldt/GenMeanDSG

max_issues_repo_name: nveldt/GenMeanDSG

max_forks_repo_name: nveldt/GenMeanDSG


text:
#=
Given a string of round, curly, and square open and closing brackets, return whether the brackets are balanced (well-formed).

For example, given the string "([])[]({})", you should return true.

Given the string "([)]" or "((()", you should return false.
=#

#=
For a string of brackets to be balanced, it needs to have as many opening brackets than closing brackets. Also, inside all these pairs of brackets, it needs to be balanced.
=#
function isbalanced(s::String)
    is_balanced = true
    if s == ""
        is_balanced = true
    else
        char_index = 1
        while is_balanced == true && char_index <= length(s)
            if s[char_index] in ['[', '(', '{']
                closingbracket = closingbrackets(s[char_index])
                closingbracket_index = findnext(isequal(closingbracket), s, char_index)
                if isnothing(closingbracket_index)
                    is_balanced = false
                else
                    is_balanced = isbalanced(s[char_index+1:closingbracket_index-1])
                end
            end
            char_index = char_index + 1
        end
    end
    return is_balanced
end

function closingbrackets(opening_bracket::Char)
    if opening_bracket == '['
        return ']'
    elseif opening_bracket == '{'
        return '}'
    elseif opening_bracket == '('
        return ')'
    else
        return error("Not a valid bracket")
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
import HTTP.WebSockets: WebSocket
function handle_rpcs(ws::WebSocket)
    while !eof(ws)
        data = readavailable(ws)
        rpc = JSON2.read(IOBuffer(data))
        dispatch(rpc, ws)
    end
end

const _msg_id = Ref{Int}(0)
msg_id() = _msg_id[] += 1

function dispatch(rpc, ws)
    meth = @match rpc.method begin
        "Runtime.enable"                    => Runtime.enable
        "Runtime.runIfWaitingForDebugger"   => Runtime.runIfWaitingForDebugger
        "Runtime.getHeapUsage"              => Runtime.getHeapUsage
        "Debugger.setPauseOnExceptions"     => Debugger.setPauseOnExceptions
        "Debugger.setAsyncCallStackDepth"   => Debugger.setAsyncCallStackDepth
        "Debugger.setBlackboxPatterns"      => Debugger.setBlackboxPatterns
        "Debugger.setBreakpointsActive"     => Debugger.setBreakpointsActive
        "Debugger.setBreakpointByUrl"       => Debugger.setBreakpointByUrl
        "Debugger.enable" => Debugger.enable
        _                 => nothing
    end
    resp = if meth === nothing
        Base.invokelatest(dispatch_unhandled, rpc)
    else
        try
        if isdefined(rpc, :params)
            Base.invokelatest(meth, rpc.params)
        else
            Base.invokelatest(meth)
        end
        catch e
            Base.showerror(stderr,e,Base.catch_backtrace())
            ()
        end
    end

    # each response has to hve the same id as the source
    resp = (;id = rpc.id, result = resp)
    @show resp
    io = IOBuffer()
    JSON2.write(io, resp)
    write(ws, take!(io))
end

function dispatch_unhandled(rpc)
    @info "Unhandled rpc" rpc
    return ()
end

module Debugger
    enable(params) = (;debuggerId = "0")
    setPauseOnExceptions(params) = ()
    setAsyncCallStackDepth(params) = ()
    setBlackboxPatterns(params) = ()
    setBreakpointsActive(params) = ()
    setBreakpointByUrl(params) = ()
    # pause
    # setSkipAllPauses
end

module Runtime
    enable() = ()
    runIfWaitingForDebugger() = ()
    getIsolateId() = ()
    getHeapUsage(params)= (;usedSize = 1024, totalSize = 4096)
end


text:
@testset "Fast cis" begin
    @test @inferred(GNSSSignals.sin_fast(π / 2)) ≈ sin(π / 2)
    @test @inferred(GNSSSignals.sin_fast(-π / 2)) ≈ sin(-π / 2)
    @test @inferred(GNSSSignals.sin_fast(0.0)) ≈ sin(0.0)
    @test @inferred(GNSSSignals.sin_fast(1π)) ≈ sin(π) atol = eps(Float64)
    @test @inferred(GNSSSignals.sin_fast(-1π)) ≈ sin(-π) atol = eps(Float64)

    @test @inferred(GNSSSignals.cos_fast(π / 2)) ≈ cos(π / 2) atol = eps(Float64)
    @test @inferred(GNSSSignals.cos_fast(-π / 2)) ≈ cos(-π / 2) atol = eps(Float64)
    @test @inferred(GNSSSignals.cos_fast(0.0)) ≈ cos(0.0)
    @test @inferred(GNSSSignals.cos_fast(1π)) ≈ cos(π)
    @test @inferred(GNSSSignals.cos_fast(-1π)) ≈ cos(-π)

    x = -π:0.01:π
    @test sqrt(sum(abs2.(GNSSSignals.sin_vfast.(x) .- sin.(x))) / length(x)) < 0.0359
    @test sqrt(sum(abs2.(GNSSSignals.sin_fast.(x) .- sin.(x))) / length(x)) < 0.000597

    @test sqrt(sum(abs2.(GNSSSignals.cos_vfast.(x) .- cos.(x))) / length(x)) < 0.03589
    @test sqrt(sum(abs2.(GNSSSignals.cos_fast.(x) .- cos.(x))) / length(x)) < 0.000597

    @test sqrt(sum(abs2.(GNSSSignals.cis_vfast.(x) .- cis.(x))) / length(x)) < 0.0507
    @test sqrt(sum(abs2.(GNSSSignals.cis_fast.(x) .- cis.(x))) / length(x)) < 0.00085

    @test get_carrier_fast_unsafe.(x) == GNSSSignals.cis_fast.(x)
    @test get_carrier_vfast_unsafe.(x) == GNSSSignals.cis_vfast.(x)
end

@testset "Generate carrier" begin

    carrier = StructArray(zeros(Complex{Int16}, 2500))

    fpcarrier!(
        carrier,
        1500Hz,
        2.5e6Hz,
        0.25, # π / 2
        start_sample = 111,
        num_samples = 2390,
        bits = Val(7)
    )

    @test sqrt(mean(abs2.(carrier.re[111:2500] ./ 1 << 7 .-
                            cos.(2π * (0:2389) * 1500Hz / 2.5e6Hz .+ π / 2)))) < 6.5e-3
    @test sqrt(mean(abs2.(carrier.im[111:2500] ./ 1 << 7 .-
                            sin.(2π * (0:2389) * 1500Hz / 2.5e6Hz .+ π / 2)))) < 6.5e-3
end


text:
function get_name(obj::TDBScale)
    return jcall(obj, "getName", JString, ())
end

function offset_from_tai(obj::TDBScale, arg0::AbsoluteDate)
    return jcall(obj, "offsetFromTAI", jdouble, (AbsoluteDate,), arg0)
end

function offset_from_tai(obj::TDBScale, arg0::FieldAbsoluteDate)
    return jcall(obj, "offsetFromTAI", RealFieldElement, (FieldAbsoluteDate,), arg0)
end

function to_string(obj::TDBScale)
    return jcall(obj, "toString", JString, ())
end


text:
#=
rect.jl
=#

using MIRT: rect

using Test: @test, @inferred

@test (@inferred rect(Int(3))) === zero(Float32)
@test (@inferred rect(Int(0))) === one(Float32)
@test (@inferred rect(0.25)) === one(Float64)
@test (@inferred rect(Float16(0.25))) === one(Float32)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -0.9705716506753843 <= q <= 1.6199082189427347)

                     add_NL_constraint(m, :(swish(0.27812717705086953 + -0.4967460082986306*swish(-0.6740257057257448 + -0.0007917455797805673*swish(-0.08584163188401472 + 0.6683752749724832*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*swish(-0.7769463282105753 + 0.25325339074048747*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + -0.883304587160544*swish(0.3709718635886192 + -0.47514710764751555*swish(-0.08584163188401472 + 0.6683752749724832*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*swish(-0.7769463282105753 + 0.25325339074048747*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) + swish(0.48450040123340576 + -0.32897580116872405*swish(-0.6740257057257448 + -0.0007917455797805673*swish(-0.08584163188401472 + 0.6683752749724832*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*swish(-0.7769463282105753 + 0.25325339074048747*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + 0.43090905819531766*swish(0.3709718635886192 + -0.47514710764751555*swish(-0.08584163188401472 + 0.6683752749724832*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*swish(-0.7769463282105753 + 0.25325339074048747*swish(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*swish(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Meganet, BenchmarkTools

const history = Pkg.dir("Meganet")*"//benchmarks//micro//bm_batchnorm.jld"

const TYPE = Float64

npixel = 500
nex = 1000
nchannel = 3

const L = getNormLayer(TYPE,[npixel,nchannel,nex],3)
const theta = initTheta(L)
const Y     = randn(TYPE,nFeatIn(L),nex)

function benchmarkJYTmv()
    funcName = "JYTmv" #TODO: pass funcName to history instead of calling it "hist"
    _,_,tmp = apply(L,theta,Y,true)
    Zout = randn(TYPE,nFeatOut(L),nex)

    #Warmup
    Z1 =  JYTmv(L,copy(Zout),(TYPE)[],theta,Y,tmp)
    trial = @benchmark JYTmv(L,copy(Zout),(TYPE)[],theta,Y,tmp)
    Meganet.updatehistory!(history, trial, "hist")
    hist = JLD.load(history, "hist")
    judge(hist)
end

function benchmarkJYmv()
    funcName = "JYmv"
    _,_,tmp = apply(L,theta,Y,true)
    dY    = randn(TYPE,nFeatIn(L),nex)
    #Warmup
    dZ, dZ = JYmv(L,dY,theta,Y,tmp)
    # @code_warntype JYmv(L, dY, theta, Y, tmp)
    trial = @benchmark JYmv(L, dY, theta, Y, tmp)
    Meganet.updatehistory!(history, trial, "hist")
    hist = JLD.load(history, "hist")
    judge(hist)
end

function benchmarkApply(L, theta, Y, history)
    funcName = "apply"
    Q = copy(Y)
    Yout,Yout,tmp = apply(L,theta,Q,true)

    @code_warntype apply(L,theta,Y,true)

    trial = @benchmark apply(L,theta,Y,true)

    Meganet.updatehistory!(history, trial, "hist")
    hist = JLD.load(history, "hist")
    judge(hist)

    return trial
end


text:
export DimensionDiagrams, EstimateDimension, EstimateDimensionMLE, EstimateDimensionANOVA, EstimateDimensionNPCA, EstimateDimensionPCA, EstimateDimensionCorrSum, EstimateDimensionBoxCounting

######################
# Dimension Diagrams #
######################
"""
DimensionDiagrams(
    data::Array{T,2},
    projective::Bool;
    methods  = [:CorrSum, :BoxCounting, :PHCurve, :NPCA, :MLE, :ANOVA, :PHCurve],
    eps_ticks = 25,
    fontsize = 16,
    lw = 4
    ) where {T <: Number}

Produces Dimension Diagrams.
* data is a matrix whose colums are the data points in Ω.
* projective = false makes diagrams in euclidean space.
* projective = true makes diagrams in projective space.
There are some optional arguments.
* methods lists the dimension estimators to be plotted.
* eps_ticks = k puts into k evenly spaces ϵ on [0,1] at which the dimension is computed.
* fontsize sets the fontsize of the axes.
* lw sets the linewidth.
"""
function DimensionDiagrams(
    data::Array{T,2},
    projective::Bool;
    diagrams  = [:CorrSum, :BoxCounting, :NPCA, :MLE, :ANOVA, :PHCurve],
    eps_ticks = 25,
    fontsize = 16,
    lw = 4
    ) where {T <: Number}

    n = size(data,1)

    cols = Colors.colormap("RdBu", mid = 0.5)
    colors = Dict("CorrSum" => cols[5], "BoxCounting"=> cols[20], "PHCurve"=> cols[35], "NPCA" => cols[55], "MLE" => cols[70], "ANOVA" => cols[85], "PHCurve" => cols[100])

    ϵ = Array{Float64}(range(0.1, length = eps_ticks, stop = 0.9))

    p = Plots.plot()
    dimension_estimates = map(diagrams) do m
        d = EstimateDimension(m, data, ϵ, projective)
        Plots.plot!(p, ϵ, d, label = string(m), lw = lw, linecolor = colors["$m"])

        d
    end

    if !projective
        y_upper = n + 1 + 0.1
    else
        y_upper = n + 0.1
    end
    y_upper = max(y_upper, maximum(vcat(dimension_estimates...)))

    Plots.plot(p,
              xlims = [0,1],
              ylims = [0,y_upper],
              titlefont = fontsize,
              tickfont = fontsize,
              labelfontsize = fontsize,
              xlabel = "epsilon",
              ylabel = "d(epsilon)",
              xscale = :none,
              yscale = :none
      )
end

#################
# Wrapper #
#################
function EstimateDimension(method::Symbol, data::Array{T,2}, ϵ::Vector{S}, projective::Bool) where {S,T<:Number}
    if method == :CorrSum
        return EstimateDimensionCorrSum(data, ϵ, projective)
    elseif method == :BoxCounting
        return EstimateDimensionBoxCounting(data, ϵ, projective)
    elseif method == :NPCA
        return EstimateDimensionNPCA(data, ϵ, projective)
    elseif method == :MLE
        return EstimateDimensionMLE(data, ϵ, projective)
    elseif method == :ANOVA
        return EstimateDimensionANOVA(data, ϵ, projective)
    elseif method == :PHCurve
        return EstimateDimensionPHCurve(data, ϵ, projective)
    else
        println("Method not supported")
        return zeros(length(ϵ))
    end
end

#################
# PCA #
#################
function EstimateDimensionPCA(data::Array{T,2}, projective::Bool) where {T <: Number}
    if !projective
        n = size(data, 1)
        data = data .- Statistics.mean.([data[i,:] for i in 1:n])
    else
        data = SafeDehomogenization(data)
        n = size(data, 1)
        data = data .- Statistics.mean.([data[i,:] for i in 1:n])
    end
    if size(data, 2) > 1
        L = log10.(LinearAlgebra.svdvals(data))
        if length(L) == 1
            if L[1] < 1e-15
                return 0
            else
                return 1
            end
        else
            i = findmin(diff(L))
            # return rank(data)
            return i[2]
        end
    else
        return 0
    end
end

#################
# Nonlinear PCA #
#################
function EstimateDimensionNPCA(data::Array{T,2}, ϵ_array::Vector{S}, projective::Bool) where {S,T<:Number}

    if !projective
        dist = ScaledEuclidean(data)
    else
        dist = ScaledFubiniStudy(data)
    end

    H = Clustering.hclust(dist, linkage = :single)
    m = size(data,2)
    return map(ϵ_array) do ϵ
        myclust = Clustering.cutree(H, h = ϵ)
        groups = unique(myclust)
        r = map(groups) do group
            index = findall(myclust .== group)
            return EstimateDimensionPCA(data[:,index], projective) * length(index)
        end
        return sum(r) / m
    end
end
function EstimateDimensionNPCA(data::Array{T,2}, projective::Bool; eps_ticks = 25) where {T<:Number}
    ϵ = Array{Float64}(range(0.1, length =eps_ticks, stop = 0.9))
    EstimateDimensionNPCA(data, ϵ, projective)
end

#########################
# Correlation Dimension #
#########################
function EstimateDimensionCorrSum(data::Array{T,2}, ϵ_array::Vector{S}, projective::Bool) where {S,T<:Number}
    n = size(data, 1)
    m = size(data, 2)
    if !projective
        dist = ScaledEuclidean(data)
    else
        dist = ScaledFubiniStudy(data)
    end

    output =  map(ϵ_array) do ϵ
        count = (sum(dist .<  ϵ) - m) / 2
        if count > 0
            return abs( log(count / (m * (m-1))))
        else
            return 0.0
        end
    end

    if !projective
        output = abs.(diff(output) ./ diff(log.(ϵ_array)))
        return [0.0; output]
    else
        output = abs.(diff(output) ./ diff(log.(sin.(ϵ_array))))
        return [0.0; output]
    end
end
function EstimateDimensionCorrSum(data::Array{T,2}, projective::Bool; eps_ticks = 25) where {T<:Number}
    ϵ = Array{Float64}(range(0.1, length = eps_ticks, stop = 0.9))
    EstimateDimensionCorrSum(data, ϵ, projective)
end

##########################
# Box Counting Dimension #
##########################
function EstimateDimensionBoxCounting(data::Array{T,2}, ϵ_array::Vector{S}, projective::Bool) where {S,T<:Number}
    n = size(data, 1)
    m = size(data, 2)

    if !projective
        u_min = minimum.([data[i,:] for i in 1:n])
        u_max = maximum.([data[i,:] for i in 1:n])

        d = abs.(u_max - u_min)
        λ = maximum(d) * √(n)

        return map(ϵ_array) do ϵ
            R = floor(Int64, λ/ϵ + 1)
            if R == 1
                return 0.0
            else
                qs = [floor.((abs.(data[:,i] - u_min) ./ d) .* R) for i in 1:m]
                return log(length(unique(qs)))/abs(log(R))
            end
        end
    else
        data = SafeDehomogenization(data)
        n = n-1
        u_min = minimum.([data[i,:] for i in 1:n])
        u_max = maximum.([data[i,:] for i in 1:n])

        i = 1
        #compute spherical angle for the box, not projective
        d = map(1:n) do i
            a = [u_min[i]; 1.0]
            b = [u_max[i]; 1.0]
            c = abs(LinearAlgebra.dot(a,b) / (LinearAlgebra.norm(a) * LinearAlgebra.norm(b)))
            if c < 1.0
                return acos(c)
            else
                return 0.0
            end
        end
        λ = maximum(d)
        return map(ϵ_array) do ϵ
            R = floor(Int64, λ/ϵ + 1)

            if R == 1
                return 0.0
            else
                D = map(1:m) do j
                    return map(1:n) do i
                        a = [u_min[i]; 1.0]
                        b = [data[i,j]; 1.0]
                        c = abs(LinearAlgebra.dot(a,b) / (LinearAlgebra.norm(a) * LinearAlgebra.norm(b)))
                        if c < 1.0
                            return acos(c)
                        else
                            return 0.0
                        end
                    end
                end
                qs = [floor.(( D[i] ./ d) .* R) for i in 1:m]
                return log(length(unique(qs)))/abs(log(R))
            end
        end
    end
end
function EstimateDimensionBoxCounting(data::Array{T,2}, projective::Bool; eps_ticks = 25) where {T<:Number}
    ϵ = Array{Float64}(range(0.1, length = 25, stop = 0.9))
    EstimateDimensionBoxCounting(data, ϵ, projective)
end

###########################
# Bickel-Levina-Algorithm #
###########################
function BL_MLE(dist::Array{T,1}, ϵ::Number, projective::Bool) where {T<:Number}
    k = length(dist)
    if !projective
        return float(abs(inv(sum(log.(ϵ .* inv.(dist))) / k)))

    else
        return float(abs(inv(sum(log.(sin(ϵ) .* inv.(sin.(dist)))) / k)))

    end
end
function EstimateDimensionMLE(data::Array{T,2}, ϵ_array::Vector{S}, projective::Bool) where {S,T<:Number}
    m = size(data, 2)
    if !projective
        dist = ScaledEuclidean(data)
    else
        dist = ScaledFubiniStudy(data)
    end
    estimators = map(ϵ_array) do ϵ
            estimator = zeros(Float64, m, 2)
            for i = 1:m
                dist_i = dist[i,[1:i-1;i+1:m]]

                if !projective
                    index = findall(dist_i .<  ϵ)
                else
                    index = findall(dist_i .<  ϵ)
                end

                k = length(index)
                if k > 1
                    estimator[i,:] = [BL_MLE(dist_i[index], ϵ, projective) * k, k]
                else
                    estimator[i,:] = [0.0, 0]
                end
            end
            return estimator
        end

    return map(estimators) do entry
        return sum(entry[:,1]) / sum(entry[:,2])
    end
end
function EstimateDimensionMLE(data::Array{T,2}, projective::Bool; eps_ticks = 25) where {T<:Number}
    ϵ = Array{Float64}(range(0.1, length = 25, stop = 0.9))
    EstimateDimensionMLE(data, ϵ, projective)
end

#################################
# Velasco-Quiroz-Diaz-Algorithm #
#################################
# Main function to compute the estimate of the dimension
# It is the algorithm from Sec. 2.1
function DQV_Estimator(dist::Array{T,2}) where {T<:Number}
    # Compute the U-statistic from Equation (3)
    m = size(dist,1)
    v = vec(dist) .- pi/2
    S = LinearAlgebra.dot(v, v)
    S /= (2 * binomial(m,2))

    # Now we have to compute the variances β from Lemma 2.5
    # Ininitialize the βs.
    # Treat d even and d odd separately
    β_even = pi^2 / 12
    β_odd = pi^2 / 4
    n_even = 0
    n_odd = 0

    # Compute the βs recursively. Since they are decreasing,
    # I only need to check for when U>β happens
    while S < β_even
        n_even = n_even + 1
        β_even = β_even - 2 / (2*n_even)^2
    end

    while S < β_odd
        n_odd = n_odd + 1
        β_odd = β_odd - 2 / (2*n_odd-1)^2
    end

    # Now I have four canditates for the closest β.
    index = findmin([abs(S-β_odd) abs(S-β_odd-2 / (2*n_odd-1)^2) abs(S-β_even) abs(S-β_even-2 / (2*n_even)^2) ])

    # return the corresponding dimension
    if index[2][2] == 1
        return 2 * n_odd + 1
    elseif index[2][2] == 2
        return maximum([2 * n_odd - 1,0])
    elseif index[2][2] == 3
        return 2 * n_even + 2
    elseif index[2][2] == 4
        return 2 * n_even
    end
end
function EstimateDimensionANOVA(data::Array{T,2}, ϵ_array::Vector{U}, projective) where {T,U <: Number}
    n = size(data, 1)
    m = size(data, 2)
    dist = zeros(T, m, m)
    D = zeros(Float64, m, m)

    if !projective
        dist = ScaledEuclidean(data)
    else
        data = hcat([LinearAlgebra.normalize(data[:,i]) for i in 1:m]...)
        dist = ScaledFubiniStudy(data)
    end

    estimators = map(1:m) do i
        P = zeros(T, n, n)
        if projective
            P = LinearAlgebra.diagm(0 => ones(n)) - data[:,i] * LinearAlgebra.transpose(data[:,i])
            data_x = (P*data)[:,[1:i-1;i+1:m]]
        else
            data_x = data[:,[1:i-1;i+1:m]] .- data[:,i]
        end
        cos_dist = FubiniStudyDistances!(D, data_x)
        return map(ϵ_array) do ϵ
            index = findall(dist[i,[1:i-1;i+1:m]] .< ϵ)
            k = length(index)
                if k > 1
                    return [DQV_Estimator(cos_dist[index,index]) * k, k]
                else
                    return [0, 0]
                end
        end
    end

    estimators = map(1:length(ϵ_array)) do i
        return [map(entry -> entry[i][1], estimators), map(entry -> entry[i][2], estimators)]
    end

    return map(1:length(ϵ_array)) do i
        sum(estimators[i][1]) / sum(estimators[i][2])
    end
end
function EstimateDimensionANOVA(data::Array{T,2}, projective::Bool; eps_ticks = 25) where {T<:Number}
    ϵ = Array{Float64}(range(0.1, length = 25, stop = 0.9))
    EstimateDimensionANOVA(data, ϵ, projective)
end

#################################
# PH curve dimension #
#################################
function EstimateDimensionPHCurve(dist::Array{T,2}) where {T <: Number}

    if sum(dist .> 0.0) > 0
        C = Eirene.eirene(dist, maxdim = 0)
        B = Eirene.barcode(C, dim = 0)
        m = size(B,1)

        if m > 1
            integral = sum([abs(B[i,2] - B[i,1]) for i in 1:(m-1)]) / (m-1)
            return abs(log(m) / log(integral))
        else
            return 0.0
        end
    else
        return 0.0
    end
end
function EstimateDimensionPHCurve(data::Array{T,2}, ϵ_array::Vector{S}, projective) where {T,S <: Number}
    if !projective
        dist = ScaledEuclidean(data)
    else
        dist = ScaledFubiniStudy(data)
    end

    H = Clustering.hclust(dist, :single)
    m = size(data,2)
    return map(ϵ_array) do ϵ

        myclust = Clustering.cutree(H, h = ϵ)
        groups = unique(myclust)
        r = map(groups) do group
            index = findall(myclust .== group)
            return EstimateDimensionPHCurve(dist[index,index]) * length(index)
        end
        return sum(r) / m
    end
end

function EstimateDimensionPHCurve(data::Array{T,2}, projective::Bool; eps_ticks = 25) where {T<:Number}
    ϵ = Array{Float64}(range(0.1, length =eps_ticks, stop = 0.9))
    EstimateDimensionPHCurve(data, ϵ, projective)
end


text:
import TracedRandom: Here

## Rewrite the following random primitives as calls to Gen distributions
const randprims = Set([
    rand, randn, randexp, randperm, randcycle, shuffle, sample,
    rand!, randn!, randexp!, randperm!, randcycle!, shuffle!, sample!
])
const primnames = Set(Symbol(fn) for fn in randprims)

## Manual addressing ##

# Forward manually addressed calls for each random primitive
for fn in randprims
    @eval @inline trace(options::Options, state, ::ManualAddress, ::typeof($fn), addr::Address, args...) =
        trace(options, state, addr::Address, $fn, args...)
    @eval @inline trace(options::Options, state, ::ManualAddress, ::typeof($fn), ::AbstractRNG, addr::Address, args...) =
        trace(options, state, addr::Address, $fn, rng, args...)
end

# Forward traced function calls that are manually addressed via rand
@inline trace(options::Options, state, ::ManualAddress, ::typeof(rand), addr::Address, fn::Function, args...) =
    trace(options, state, addr::Address, fn, args...)

# Splice function calls with `here` addresses into the current context
@inline trace(options::Options, state, ::ManualAddress, ::typeof(rand), ::Here, fn::Function, args...) =
    splice(options, state, fn, args...)

## Strip RNGs supplied to primitives

for fn in randprims
    @eval @inline trace(options::Options, state, addr::Address, ::typeof($fn), ::AbstractRNG, args...) =
        trace(options, state, addr::Address, $fn, args...)
end

## rand ##

# Default to uniformly sampling Float64 values from [0, 1]
@inline trace(options::Options, state, addr::Address, ::typeof(rand)) =
    trace(options, state, addr::Address, rand, Float64)
@inline trace(options::Options, state, addr::Address, ::typeof(rand), d::Integer, dims::Integer...) =
    trace(options, state, addr::Address, rand, Float64, d, dims...)

# rand for numeric types
@inline trace(::Options, state, addr::Address, ::typeof(rand), T::Type{<:Real}) =
    Gen.traceat(state, TypedScalarDistribution(T), (), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), T::Type{<:Real}, ::Tuple{}) =
    Gen.traceat(state, TypedScalarDistribution(T), (), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), T::Type{<:Real}, dims::Dims) =
    Gen.traceat(state, TypedArrayDistribution(T, dims), (), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), T::Type{<:Real}, d::Integer, dims::Integer...) =
    Gen.traceat(state, TypedArrayDistribution(T, d, dims...), (), addr)

# rand for indexable collections
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T) where {T <: Indexable} =
    Gen.traceat(state, LabeledUniform{eltype(T)}(), (c,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T, ::Tuple{}) where {T <: Indexable} =
    Gen.traceat(state, LabeledUniform{eltype(T)}(), (c,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T, dims::Dims) where {T <: Indexable} =
    Gen.traceat(state, ArrayedDistribution(LabeledUniform{eltype(T)}(), dims), (c,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T, d::Integer, dims::Integer...) where {T <: Indexable} =
    Gen.traceat(state, ArrayedDistribution(LabeledUniform{eltype(T)}(), d, dims...), (c,), addr)

# rand for set-like collections
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T) where {T <: Setlike} =
    Gen.traceat(state, SetUniform{eltype(T)}(), (c,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T, ::Tuple{}) where {T <: Setlike} =
    Gen.traceat(state, SetUniform{eltype(T)}(), (c,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T, dims::Dims) where {T <: Setlike} =
    Gen.traceat(state, ArrayedDistribution(SetUniform{eltype(T)}(), dims), (c,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), c::T, d::Integer, dims::Integer...) where {T <: Setlike} =
    Gen.traceat(state, ArrayedDistribution(SetUniform{eltype(T)}(), d, dims...), (c,), addr)

# rand for Distributions.jl distributions
@inline trace(::Options, state, addr::Address, ::typeof(rand), dist::D) where {D <: Distributions.Distribution} =
    Gen.traceat(state, WrappedDistribution(dist), params(dist), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), dist::D, ::Tuple{}) where {D <: Distributions.Distribution} =
    Gen.traceat(state, WrappedDistribution(dist), params(dist), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), dist::D, dims::Dims) where {D <: Distributions.Distribution} =
    Gen.traceat(state, ArrayedDistribution(WrappedDistribution(dist), dims), params(dist), addr)
@inline trace(::Options, state, addr::Address, ::typeof(rand), dist::D, d::Integer, dims::Integer...) where {D <: Distributions.Distribution} =
    Gen.traceat(state, ArrayedDistribution(WrappedDistribution(dist), d, dims...), params(dist), addr)

# rand! (simulate in-place allocation)
@inline trace(options::Options, state, addr::Address, ::typeof(rand!), A::AbstractArray{T}) where {T} =
    copyto!(A, trace(options, state, addr::Address, rand, T, size(A)))
@inline trace(options::Options, state, addr::Address, ::typeof(rand!), A::AbstractArray, S) =
    copyto!(A, trace(options, state, addr::Address, rand, S, size(A)))

## randn ##

# Default to sampling Float64 values
@inline trace(options::Options, state, addr::Address, ::typeof(randn)) =
    trace(options, state, addr::Address, randn, Float64)
@inline trace(options::Options, state, addr::Address, ::typeof(randn), dims::Dims) =
    trace(options, state, addr::Address, randn, Float64, dims)
@inline trace(options::Options, state, addr::Address, ::typeof(randn), d::Integer, dims::Integer...) =
    trace(options, state, addr::Address, randn, Float64, d, dims...)

# Forward sampling of all other Float types
@inline trace(::Options, state, addr::Address, ::typeof(randn), T::Type{<:AbstractFloat}, args...) =
    trace(options, state, addr::Address, randn, Float64, args...)

# Sampling of Float64 values
@inline trace(::Options, state, addr::Address, ::typeof(randn), T::Type{Float64}) =
    Gen.traceat(state, Gen.normal, (0, 1), addr)
@inline trace(::Options, state, addr::Address, ::typeof(randn), T::Type{Float64}, ::Tuple{}) =
    Gen.traceat(state, Gen.normal, (0, 1), addr)
@inline trace(::Options, state, addr::Address, ::typeof(randn), T::Type{Float64}, dims::Dims) =
    Gen.traceat(state, ArrayedDistribution(Gen.normal, dims), (0, 1), addr)
@inline trace(::Options, state, addr::Address, ::typeof(randn), T::Type{Float64}, d::Integer, dims::Integer...) =
    Gen.traceat(state, ArrayedDistribution(Gen.normal, d, dims...), (0, 1), addr)

# randn! (simulate in-place allocation)
@inline trace(options::Options, state, addr::Address, ::typeof(randn!), A::AbstractArray{T}) where {T} =
    copyto!(A, trace(options, state, addr::Address, randn, T, size(A)))

## randexp ##

# Default to sampling Float64 values
@inline trace(options::Options, state, addr::Address, ::typeof(randexp)) =
    trace(options, state, addr::Address, randexp, Float64)
@inline trace(options::Options, state, addr::Address, ::typeof(randexp), dims::Dims) =
    trace(options, state, addr::Address, randexp, Float64, dims)
@inline trace(options::Options, state, addr::Address, ::typeof(randexp), d::Integer, dims::Integer...) =
    trace(options, state, addr::Address, randexp, Float64, d, dims...)

# Forward sampling of all other Float types
@inline trace(::Options, state, addr::Address, ::typeof(randexp), T::Type{<:AbstractFloat}, args...) =
    trace(options, state, addr::Address, randexp, Float64, args...)

# Sampling of Float64 values
@inline trace(::Options, state, addr::Address, ::typeof(randexp), T::Type{Float64}) =
    Gen.traceat(state, Gen.exponential, (1,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(randexp), T::Type{Float64}, ::Tuple{}) =
    Gen.traceat(state, Gen.exponential, (1,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(randexp), T::Type{Float64}, dims::Dims) =
    Gen.traceat(state, ArrayedDistribution(Gen.exponential, dims), (1,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(randexp), T::Type{Float64}, d::Integer, dims::Integer...) =
    Gen.traceat(state, ArrayedDistribution(Gen.exponential, d, dims...), (1,), addr)

# randexp! (simulate in-place allocation)
@inline trace(options::Options, state, addr::Address, ::typeof(randexp!), A::AbstractArray{T}) where {T} =
    copyto!(A, trace(options, state, addr::Address, randexp, T, size(A)))

## randperm, randcycle. shuffle ##

# Forward to corresponding Gen distributions
@inline trace(options::Options, state, addr::Address, ::typeof(randperm), n::Integer) =
    Gen.traceat(state, RandomPermutation(n), (), addr)
@inline trace(options::Options, state, addr::Address, ::typeof(randcycle), n::Integer) =
    Gen.traceat(state, RandomCycle(n), (), addr)
@inline trace(options::Options, state, addr::Address, ::typeof(shuffle), v::AbstractArray) =
    Gen.traceat(state, RandomShuffle(v), (), addr)

# Simulate in-place versions
@inline trace(options::Options, state, addr::Address, ::typeof(randperm!), A::AbstractArray{<:Integer}) =
    copyto!(A, trace(options, state, addr::Address, randperm, length(A)))
@inline trace(options::Options, state, addr::Address, ::typeof(randcycle!), A::AbstractArray{<:Integer}) =
    copyto!(A, trace(options, state, addr::Address, randcycle, length(A)))
@inline trace(options::Options, state, addr::Address, ::typeof(shuffle!), v::AbstractArray) =
    copyto!(v, trace(options, state, addr::Address, shuffle, copy(v)))

## sample ##

# Unweighted sample
@inline trace(::Options, state, addr::Address, ::typeof(sample), a::AbstractArray{T}) where {T} =
    Gen.traceat(state, LabeledUniform{T}(), (a,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(sample), a::AbstractArray{T}, n::Integer) where {T} =
    Gen.traceat(state, ArrayedDistribution(LabeledUniform{T}(), n), (a,), addr)
@inline trace(::Options, state, addr::Address, ::typeof(sample), a::AbstractArray{T}, dims::Dims) where {T} =
    Gen.traceat(state, ArrayedDistribution(LabeledUniform{T}(), dims), (a,), addr)

# Weighted sample
@inline trace(::Options, state, addr::Address, ::typeof(sample), w::AbstractWeights) =
    Gen.traceat(state, Categorical(), (Vector(w) ./ sum(w)), addr)
@inline trace(::Options, state, addr::Address, ::typeof(sample), a::AbstractArray{T}, w::AbstractWeights) where {T} =
    Gen.traceat(state, LabeledCategorical{T}(), (a, w), addr)
@inline trace(::Options, state, addr::Address, ::typeof(sample), a::AbstractArray{T}, w::AbstractWeights, n::Integer) where {T} =
    Gen.traceat(state, ArrayedDistribution(LabeledCategorical{T}(), n), (a, w), addr)
@inline trace(::Options, state, addr::Address, ::typeof(sample), a::AbstractArray{T}, w::AbstractWeights, dims::Dims) where {T} =
    Gen.traceat(state, ArrayedDistribution(LabeledCategorical{T}(), dims), (a, w), addr)

# sample! (simulate in-place allocation)
@inline trace(options::Options, state, addr::Address, ::typeof(sample!), a::AbstractArray, x::AbstractArray) =
    copyto!(x, trace(options, state, addr::Address, sample, a, size(x)))
@inline trace(options::Options, state, addr::Address, ::typeof(sample!), a::AbstractArray, w::AbstractWeights, x::AbstractArray) =
    copyto!(x, trace(options, state, addr::Address, sample, a, w, size(x)))


text:
using Test
using Random
using LinearAlgebra
using Combinatorics
using TensorKit
using TensorKit: ProductSector, fusiontensor
using TensorOperations
TensorOperations.disable_cache() # avoids memory overflow during CI?
using TupleTools
using TupleTools: StaticLength
using Base.Iterators: take, product

Random.seed!(1234)

smallset(::Type{G}) where {G<:Sector} = take(values(G), 6)
function smallset(::Type{ProductSector{Tuple{G1,G2}}}) where {G1,G2}
    iter = product(smallset(G1),smallset(G2))
    s = collect(i × j for (i,j) in iter if dim(i)*dim(j) <= 6)
    return length(s) > 6 ? rand(s, 6) : s
end
function smallset(::Type{ProductSector{Tuple{G1,G2,G3}}}) where {G1,G2,G3}
    iter = product(smallset(G1),smallset(G2),smallset(G3))
    s = collect(i × j × k for (i,j,k) in iter if dim(i)*dim(j)*dim(k) <= 6)
    return length(s) > 6 ? rand(s, 6) : s
end
function randsector(::Type{G}) where {G<:Sector}
    s = collect(smallset(G))
    a = rand(s)
    while a == one(a) # don't use trivial label
        a = rand(s)
    end
    return a
end
function hasfusiontensor(G::Type{<:Sector})
    try
        fusiontensor(one(G), one(G), one(G))
        return true
    catch
        return false
    end
end

include("sectors.jl")
include("fusiontrees.jl")
include("spaces.jl")
include("tensors.jl")


text:
function sample!(fin_el::Triangle, s::Matrix{Float64} , wt::Vector{Float64})
  #
  # This subroutine returns the local coordinates and weighting coefficients
  # of the integrating points.
  #

  root3 = 1.0 / sqrt(3.0)
  r15 = 0.2 * sqrt(15.0)
  nip = size(s,1)

  w = [5.0/9.0, 8.0/9.0, 5.0/9.0]
  v = [5.0/9.0*w; 8.0/9.0*w; 5.0/9.0*w]

  if nip == 1
    s[1,1] = 0.333333333333333
    s[1,2] = 0.333333333333333
    wt[1] = 0.500000000000000
  elseif nip == 3
    s[1,1] = 0.500000000000000
    s[1,2] = 0.500000000000000
    s[2,1] = 0.500000000000000
    s[2,2] = 0.000000000000000
    s[3,1] = 0.000000000000000
    s[3,2] = 0.500000000000000
    wt[1:3] .= 0.333333333333333
    wt[:] .= 0.5*wt
  elseif nip == 4
    s[1,1] = 0.6
    s[1,2] = 0.2
    s[2,1] = 0.2
    s[2,2] = 0.6
    s[3,1] = 0.2
    s[3,2] = 0.2
    s[4,1] = 0.333333333333333
    s[4,2] = 0.333333333333333
    wt[1:3] .= 0.520833333333333
    wt[4] = -0.5625
    wt[:] .= 0.5*wt
  elseif nip == 6
    s[1,1] = 0.816847572980459
    s[1,2] = 0.091576213509771
    s[2,1] = 0.091576213509771
    s[2,2] = 0.816847572980459
    s[3,1] = 0.091576213509771
    s[3,2] = 0.091576213509771
    s[4,1] = 0.108103018168070
    s[4,2] = 0.445948490915965
    s[5,1] = 0.445948490915965
    s[5,2] = 0.108103018168070
    s[6,1] = 0.445948490915965
    s[6,2] = 0.445948490915965
    wt[1:3] .= 0.109951743655322
    wt[4:6] .= 0.223381589678011
    wt[:] .= 0.5*wt
  elseif nip == 7
    s[1,1] = 0.333333333333333
    s[1,2] = 0.333333333333333
    s[2,1] = 0.797426985353087
    s[2,2] = 0.101286507323456
    s[3,1] = 0.101286507323456
    s[3,2] = 0.797426985353087
    s[4,1] = 0.101286507323456
    s[4,2] = 0.101286507323456
    s[5,1] = 0.470142064105115
    s[5,2] = 0.059715871789770
    s[6,1] = 0.059715871789770
    s[6,2] = 0.470142064105115
    s[7,1] = 0.470142064105115
    s[7,2] = 0.470142064105115
    wt[1] = 0.225000000000000
    wt[2:4] .= 0.125939180544827
    wt[5:7] .= 0.132394152788506
    wt[:] .= 0.5*wt
  elseif nip == 12
    s[1,1] = 0.873821971016996
    s[1,2] = 0.063089014491502
    s[2,1] = 0.063089014491502
    s[2,2] = 0.873821971016996
    s[3,1] = 0.063089014491502
    s[3,2] = 0.063089014491502
    s[4,1] = 0.501426509658179
    s[4,2] = 0.249286745170910
    s[5,1] = 0.249286745170910
    s[5,2] = 0.501426509658179
    s[6,1] = 0.249286745170910
    s[6,2] = 0.249286745170910
    s[7,1] = 0.053145049844817
    s[7,2] = 0.310352451033784
    s[8,1] = 0.310352451033784
    s[8,2] = 0.053145049844817
    s[9,1] = 0.053145049844817
    s[9,2] = 0.636502499121398
    s[10,1] = 0.310352451033784
    s[10,2] = 0.636502499121398
    s[11,1] = 0.636502499121398
    s[11,2] = 0.053145049844817
    s[12,1] = 0.636502499121398
    s[12,2] = 0.310352451033784
    wt[1:3] .= 0.050844906370207
    wt[4:6] .= 0.116786275726379
    wt[7:12] .= 0.082851075618374
    wt[:] .= 0.5*wt
  elseif nip == 16
    s[1,1] = 0.333333333333333
    s[1,2] = 0.333333333333333
    s[2,1] = 0.658861384496478
    s[2,2] = 0.170569307751761
    s[3,1] = 0.170569307751761
    s[3,2] = 0.658861384496478
    s[4,1] = 0.170569307751761
    s[4,2] = 0.170569307751761
    s[5,1] = 0.898905543365938
    s[5,2] = 0.050547228317031
    s[6,1] = 0.050547228317031
    s[6,2] = 0.898905543365938
    s[7,1] = 0.050547228317031
    s[7,2] = 0.050547228317031
    s[8,1] = 0.081414823414554
    s[8,2] = 0.459292588292723
    s[9,1] = 0.459292588292723
    s[9,2] = 0.081414823414554
    s[10,1] = 0.459292588292723
    s[10,2] = 0.459292588292723
    s[11,1] = 0.008394777409958
    s[11,2] = 0.263112829634638
    s[12,1] = 0.008394777409958
    s[12,2] = 0.728492392955404
    s[13,1] = 0.263112829634638
    s[13,2] = 0.008394777409958
    s[14,1] = 0.263112829634638
    s[14,2] = 0.728492392955404
    s[15,1] = 0.728492392955404
    s[15,2] = 0.008394777409958
    s[16,1] = 0.728492392955404
    s[16,2] = 0.263112829634638
    wt[1] = 0.144315607677787
    wt[2:4] .= 0.103217370534718
    wt[5:7] .= 0.032458497623198
    wt[8:10] .= 0.095091634267284
    wt[11:16] .= 0.027230314174435
    wt[:] .= 0.5*wt
  else
    println("Wrong number of integrating points for a triangle.")
  end
end