text:
"""
    ABCDParams

A structure holding parameters for ABCD graph generator. Fields:
* w::Vector{Int}:             a sorted in descending order list of vertex degrees
* s::Vector{Int}:             a sorted in descending order list of cluster sizes
* μ::Union{Float64, Nothing}: mixing parameter
* ξ::Union{Float64, Nothing}: background graph fraction
* isCL::Bool:                 if `true` a Chung-Lu model is used, otherwise configuration model
* islocal::Bool:              if `true` mixing parameter restriction is cluster local, otherwise
                              it is only global

Exactly one of ξ and μ must be passed as `Float64`. Also if `ξ` is passed then
`islocal` must be `false`.

The base ABCD graph is generated when ξ is passed and `isCL` is set to `false`.
"""
struct ABCDParams
    w::Vector{Int}
    s::Vector{Int}
    μ::Union{Float64, Nothing}
    ξ::Union{Float64, Nothing}
    isCL::Bool
    islocal::Bool

    function ABCDParams(w, s, μ, ξ, isCL, islocal)
        length(w) == sum(s) || throw(ArgumentError("inconsistent data"))
        if !isnothing(μ)
            0 ≤ μ ≤ 1 || throw(ArgumentError("inconsistent data on μ"))
        end
        if !isnothing(ξ)
            0 ≤ ξ ≤ 1 || throw(ArgumentError("inconsistent data ξ"))
            if islocal
                throw(ArgumentError("when ξ is provided local model is not allowed"))
            end
        end
        if isnothing(μ) && isnothing(ξ)
            throw(ArgumentError("inconsistent data: either μ or ξ must be provided"))
        end

        if !(isnothing(μ) || isnothing(ξ))
            throw(ArgumentError("inconsistent data: only μ or ξ may be provided"))
        end

        new(sort(w, rev=true),
            sort(s, rev=true),
            μ, ξ, isCL, islocal)
    end
end

function randround(x)
    d = floor(Int, x)
    d + (rand() < x - d)
end

function populate_clusters(params::ABCDParams)
    w, s = params.w, params.s
    if isnothing(params.ξ)
        mul = 1.0 - params.μ
    else
        n = length(w)
        ϕ = 1.0 - sum((sl/n)^2 for sl in s)
        mul = 1.0 - params.ξ*ϕ
    end
    @assert length(w) == sum(s)
    @assert 0 ≤ mul ≤ 1
    @assert issorted(w, rev=true)
    @assert issorted(s, rev=true)

    slots = copy(s)
    clusters = Int[]
    j = 0
    for (i, vw) in enumerate(w)
        while j + 1 ≤ length(s) && mul * vw + 1 ≤ s[j + 1]
            j += 1
        end
        j == 0 && throw(ArgumentError("could not find a large enough cluster for vertex of weight $vw"))
        wts = Weights(view(slots, 1:j))
        wts.sum == 0 && throw(ArgumentError("could not find an empty slot for vertex of weight $vw"))
        loc = sample(1:j, wts)
        push!(clusters, loc)
        slots[loc] -= 1
    end
    clusters
end

function CL_model(clusters, params)
    @assert params.isCL
    w, s, μ = params.w, params.s, params.μ
    cluster_weight = zeros(Int, length(s))
    for i in axes(w, 1)
        cluster_weight[clusters[i]] += w[i]
    end
    total_weight = sum(cluster_weight)
    if params.islocal
        ξl = @. μ / (1.0 - cluster_weight / total_weight)
        maximum(ξl) >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
    else
        if isnothing(params.ξ)
            ξg = μ / (1.0 - sum(x -> x^2, cluster_weight) / total_weight^2)
            ξg >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
        else
            ξg = params.ξ
        end
    end

    wf = float.(w)
    edges = Set{Tuple{Int, Int}}()
    for i in axes(s, 1)
        local_edges = Set{Tuple{Int, Int}}()
        idxᵢ = findall(==(i), clusters)
        wᵢ = wf[idxᵢ]
        ξ = params.islocal ? ξl[i] : ξg
        m = randround((1-ξ) * sum(wᵢ) / 2)
        ww = Weights(wᵢ)
        while length(local_edges) < m
            a = sample(idxᵢ, ww, m - length(local_edges))
            b = sample(idxᵢ, ww, m - length(local_edges))
            for (p, q) in zip(a, b)
                p != q && push!(local_edges, minmax(p, q))
            end
        end
        union!(edges, local_edges)
    end
    wwt = if params.islocal
        Weights([ξl[clusters[i]]*x for (i,x) in enumerate(wf)])
    else
        Weights(ξg * wf)
    end
    while 2*length(edges) < total_weight
        a = sample(axes(w, 1), wwt, randround(total_weight / 2) - length(edges))
        b = sample(axes(w, 1), wwt, randround(total_weight / 2) - length(edges))
        for (p, q) in zip(a, b)
            p != q && push!(edges, minmax(p, q))
        end
    end
    edges
end

function config_model(clusters, params)
    @assert !params.isCL
    w, s, μ = params.w, params.s, params.μ

    cluster_weight = zeros(Int, length(s))
    for i in axes(w, 1)
        cluster_weight[clusters[i]] += w[i]
    end
    total_weight = sum(cluster_weight)
    if params.islocal
        ξl = @. μ / (1.0 - cluster_weight / total_weight)
        maximum(ξl) >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
        w_internal_raw = [w[i] * (1 - ξl[clusters[i]]) for i in axes(w, 1)]
    else
        if isnothing(params. ξ)
            ξg = μ / (1.0 - sum(x -> x^2, cluster_weight) / total_weight^2)
            ξg >= 1 && throw(ArgumentError("μ is too large to generate a graph"))
        else
            ξg = params.ξ
        end
        w_internal_raw = [w[i] * (1 - ξg) for i in axes(w, 1)]
    end

    clusterlist = [Int[] for i in axes(s, 1)]
    for i in axes(clusters, 1)
        push!(clusterlist[clusters[i]], i)
    end

    edges = Set{Tuple{Int, Int}}()

    unresolved_collisions = 0
    w_internal = zeros(Int, length(w_internal_raw))
    for cluster in clusterlist
        maxw_idx = argmax(view(w_internal_raw, cluster))
        wsum = 0
        for i in axes(cluster, 1)
            if i != maxw_idx
                neww = randround(w_internal_raw[cluster[i]])
                w_internal[cluster[i]] = neww
                wsum += neww
            end
        end
        maxw = floor(Int, w_internal_raw[cluster[maxw_idx]])
        w_internal[cluster[maxw_idx]] = maxw + (isodd(wsum) ? iseven(maxw) : isodd(maxw))

        stubs = Int[]
        for i in cluster
            for j in 1:w_internal[i]
                push!(stubs, i)
            end
        end
        @assert sum(w_internal[cluster]) == length(stubs)
        shuffle!(stubs)
        local_edges = Set{Tuple{Int, Int}}()
        recycle = Tuple{Int,Int}[]
        for i in 1:2:length(stubs)
            e = minmax(stubs[i], stubs[i+1])
            if (e[1] == e[2]) || (e in local_edges)
                push!(recycle, e)
            else
                push!(local_edges, e)
            end
        end
        last_recycle = length(recycle)
        recycle_counter = last_recycle
        while !isempty(recycle)
            recycle_counter -= 1
            if recycle_counter < 0
                if length(recycle) < last_recycle
                    last_recycle = length(recycle)
                    recycle_counter = last_recycle
                else
                    break
                end
            end
            p1 = popfirst!(recycle)
            from_recycle = 2 * length(recycle) / length(stubs)
            success = false
            for _ in 1:2:length(stubs)
                p2 = if rand() < from_recycle
                    used_recycle = true
                    recycle_idx = rand(axes(recycle, 1))
                    recycle[recycle_idx]
                else
                    used_recycle = false
                    rand(local_edges)
                end
                if rand() < 0.5
                    newp1 = minmax(p1[1], p2[1])
                    newp2 = minmax(p1[2], p2[2])
                else
                    newp1 = minmax(p1[1], p2[2])
                    newp2 = minmax(p1[2], p2[1])
                end
                if newp1 == newp2
                    good_choice = false
                elseif (newp1[1] == newp1[2]) || (newp1 in local_edges)
                    good_choice = false
                elseif (newp2[1] == newp2[2]) || (newp2 in local_edges)
                    good_choice = false
                else
                    good_choice = true
                end
                if good_choice
                    if used_recycle
                        recycle[recycle_idx], recycle[end] = recycle[end], recycle[recycle_idx]
                        pop!(recycle)
                    else
                        pop!(local_edges, p2)
                    end
                    success = true
                    push!(local_edges, newp1)
                    push!(local_edges, newp2)
                    break
                end
            end
            success || push!(recycle, p1)
        end
        old_len = length(edges)
        union!(edges, local_edges)
        @assert length(edges) == old_len + length(local_edges)
        @assert 2 * (length(local_edges) + length(recycle)) == length(stubs)
        for (a, b) in recycle
            w_internal[a] -= 1
            w_internal[b] -= 1
        end
        unresolved_collisions += length(recycle)
    end

    if unresolved_collisions > 0
        println("Unresolved_collisions: ", unresolved_collisions,
                "; fraction: ", 2 * unresolved_collisions / total_weight)
    end

    stubs = Int[]
    for i in axes(w, 1)
        for j in w_internal[i]+1:w[i]
            push!(stubs, i)
        end
    end
    @assert sum(w) == length(stubs) + sum(w_internal)
    shuffle!(stubs)
    global_edges = Set{Tuple{Int, Int}}()
    recycle = Tuple{Int,Int}[]
    for i in 1:2:length(stubs)
        e = minmax(stubs[i], stubs[i+1])
        if (e[1] == e[2]) || (e in global_edges) || (e in edges)
            push!(recycle, e)
        else
            push!(global_edges, e)
        end
    end
    last_recycle = length(recycle)
    recycle_counter = last_recycle
    while !isempty(recycle)
        recycle_counter -= 1
        if recycle_counter < 0
            if length(recycle) < last_recycle
                last_recycle = length(recycle)
                recycle_counter = last_recycle
            else
                break
            end
        end
        p1 = pop!(recycle)
        from_recycle = 2 * length(recycle) / length(stubs)
        p2 = if rand() < from_recycle
            i = rand(axes(recycle, 1))
            recycle[i], recycle[end] = recycle[end], recycle[i]
            pop!(recycle)
        else
            x = rand(global_edges)
            pop!(global_edges, x)
        end
        if rand() < 0.5
            newp1 = minmax(p1[1], p2[1])
            newp2 = minmax(p1[2], p2[2])
        else
            newp1 = minmax(p1[1], p2[2])
            newp2 = minmax(p1[2], p2[1])
        end
        for newp in (newp1, newp2)
            if (newp[1] == newp[2]) || (newp in global_edges) || (newp in edges)
                push!(recycle, newp)
            else
                push!(global_edges, newp)
            end
        end
    end
    old_len = length(edges)
    union!(edges, global_edges)
    @assert length(edges) == old_len + length(global_edges)
    if isempty(recycle)
        @assert 2 * length(global_edges) == length(stubs)
    else
        last_recycle = length(recycle)
        recycle_counter = last_recycle
        while !isempty(recycle)
            recycle_counter -= 1
            if recycle_counter < 0
                if length(recycle) < last_recycle
                    last_recycle = length(recycle)
                    recycle_counter = last_recycle
                else
                    break
                end
            end
            p1 = pop!(recycle)
            x = rand(edges)
            p2 = pop!(edges, x)
            if rand() < 0.5
                newp1 = minmax(p1[1], p2[1])
                newp2 = minmax(p1[2], p2[2])
            else
                newp1 = minmax(p1[1], p2[2])
                newp2 = minmax(p1[2], p2[1])
            end
            for newp in (newp1, newp2)
                if (newp[1] == newp[2]) || (newp in edges)
                    push!(recycle, newp)
                else
                    push!(edges, newp)
                end
            end
        end
    end
    if !isempty(recycle)
        unresolved_collisions = length(recycle)
        println("Very hard graph. Failed to generate ", unresolved_collisions,
                "edges; fraction: ", 2 * unresolved_collisions / total_weight)
    end
    return edges
end

"""
    gen_graph(params::ABCDParams)

Generate ABCD graph following parameters specified in `params`.

Return a named tuple containing a set of edges of the graph and a list of cluster
assignments of the vertices.
The ordering of vertices and clusters is in descending order (as in `params`).
"""
function gen_graph(params::ABCDParams)
    clusters = populate_clusters(params)
    edges = params.isCL ? CL_model(clusters, params) : config_model(clusters, params)
    (edges=edges, clusters=clusters)
end


text:
#===============================================================================

    Input.jl

    abstract type for data and analysis parameters
    and methods for both AmpImput and McInput structs

    Author: Tom Price
    Date:   July 2019

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

#===============================================================================
    type definition >>
===============================================================================#

abstract type Input end

#===============================================================================
    function definitions >>
===============================================================================#

ext: jl

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
module TreeShaker

# Deps
using SnoopCompile # for snooping
using Random # for snoop log names
using Pkg # context/manifest info

# Includes
include("meat.jl")

# Exports
export shake

end # module


text:
#=
Direct estimation of the probability of for a diffusion with X(0) = a satisfying
X(T) ∈ (b, ∞) for the double well potential V(x) = (x²-1)².
=#

using Distributed
nw = 4; # number of workers
addprocs(nw);

using Statistics
using HypothesisTests
using Printf
@everywhere include("doublewell_setup.jl")
opts = MDOptions(n_iters=nΔt)

Random.seed!(100)
@everywhere n_trials = 10^4;
n_in_set = @sync @distributed (+) for j in 1:n_trials
    X = copy(x₀);
    sample_trajectory!(X, sampler, options=opts);
    f(X)
end
CI = confint(BinomialTest(n_in_set, n_trials));
@printf("95%% CI = (%g,%g)\n", CI[1], CI[2])


text:
using Diffusion
using Test
using StaticArrays

@testset "Diffusion.jl" begin
    # tests for absorbing and reflective boundaries
    @test absorb_boundary(SA[3,4],6) == SA[3,4]
    @test absorb_boundary(SA[3,4],5) === nothing
    @test absorb_boundary(SA[3,4],4) === nothing
    @test reflect_boundary(SA[3,4],6) == SA[3,4]
    @test reflect_boundary(SA[3,4],5) ≈ SA[3,4]
    @test reflect_boundary(SA[4,4],sqrt(18)) ≈ SA[2,2]

    for r in 0.5:0.1:1.2, θ in 0:π/8:2π
        pt = SA[r*cos(θ), r*sin(θ)]
        rbound = 0.8
        rreflect = 2*rbound - r
        @test reflect_boundary(pt, rbound) ≈ (r <= rbound ? pt : SA[rreflect*cos(θ), rreflect*sin(θ)])
    end

    #tests for simulating the movement of molecules
    @test length(simmolreflect(1, 4, SA[0.,0.], 3)) == 5
    @test typeof(simmolreflect(1, 4, SA[1.,1.], 4)) == Vector{SVector{2, Float64}}
    @test length(simmolabsorb(1, 4, SA[0.,0.], 3)) <= 5
    @test typeof(simmolabsorb(1, 4, SA[1.,1.], 4)) == Vector{SVector{2, Float64}}
    @test length(nmolreflect(1, 4, SA[1., 1.], 3, 100)) == 100
    @test typeof(nmolreflect(1, 4, SA[1., 1.], 3, 100)) == Vector{Vector{typeof(SA[1.,1.])}}
    @test length(nmolabsorb(1, 4, SA[1., 1.], 3, 100)) == 100
    @test typeof(nmolabsorb(1, 4, SA[1., 1.], 3, 100)) == Vector{Vector{typeof(SA[1.,1.])}}
end


text:
module MeasureTheory

using Random

using ConcreteStructs
using MLStyle
using NestedTuples
using TransformVariables
const TV = TransformVariables

import Base
import Distributions
const Dists = Distributions

export ≪
export sampletype

export AbstractMeasure
using InfiniteArrays
using ConcreteStructs
using DynamicIterators
using KeywordCalls
using ConstructionBase
using Accessors

const ∞ = InfiniteArrays.∞

export ∞

export as

abstract type AbstractMeasure end

sampletype(μ::AbstractMeasure) = typeof(testvalue(μ))

# sampletype(μ::AbstractMeasure) = sampletype(basemeasure(μ))

"""
    logdensity(μ::AbstractMeasure [, ν::AbstractMeasure], x::X)

Compute the logdensity of the measure μ at the point x. This is the standard way
to define `logdensity` for a new measure. the base measure is implicit here, and
is understood to be `basemeasure(μ)`.
"""
function logdensity end

include("const.jl")
include("exp.jl")
include("domains.jl")
include("utils.jl")
include("traits.jl")
include("absolutecontinuity.jl")
include("parameterized.jl")
include("macros.jl")
include("resettablerng.jl")

include("primitive.jl")
include("primitives/counting.jl")
include("primitives/lebesgue.jl")
include("primitives/dirac.jl")
include("primitives/trivial.jl")

include("combinators/weighted.jl")
include("combinators/superpose.jl")
include("combinators/product.jl")
include("combinators/for.jl")
include("combinators/power.jl")
include("combinators/transforms.jl")
include("combinators/spikemixture.jl")
include("combinators/chain.jl")
include("kernel.jl")
include("combinators/likelihood.jl")
include("combinators/pointwise.jl")

include("distributions.jl")
include("rand.jl")

include("parameterized/normal.jl")
include("parameterized/studentt.jl")
include("parameterized/cauchy.jl")
include("parameterized/laplace.jl")
include("parameterized/uniform.jl")
include("parameterized/beta.jl")
include("parameterized/dirichlet.jl")
include("parameterized/gumbel.jl")
include("parameterized/exponential.jl")
include("parameterized/mvnormal.jl")
# include("parameterized/inverse-gamma.jl")
include("parameterized/bernoulli.jl")
include("parameterized/poisson.jl")
include("parameterized/binomial.jl")
include("parameterized/multinomial.jl")
include("parameterized/lkj-cholesky.jl")
include("parameterized/negativebinomial.jl")

include("transforms/corrcholesky.jl")

include("density.jl")
# include("pushforward.jl")

include("distproxy.jl")
end # module


text:
using OpenTelemetryProto: OpentelemetryClients
using Test

@testset "OpenTelemetryProto.jl" begin
    trace_req =
        OpentelemetryClients.opentelemetry.proto.collector.trace.v1.ExportTraceServiceRequest(
            resource_spans = [
                OpentelemetryClients.opentelemetry.proto.trace.v1.ResourceSpans(
                    resource = OpentelemetryClients.opentelemetry.proto.resource.v1.Resource(
                        dropped_attributes_count = 0,
                    ),
                    instrumentation_library_spans = [
                        OpentelemetryClients.opentelemetry.proto.trace.v1.InstrumentationLibrarySpans(
                            instrumentation_library = OpentelemetryClients.opentelemetry.proto.common.v1.InstrumentationLibrary(;
                                name = "julia",
                                version = "1.0.0",
                            ),
                            spans = [
                                OpentelemetryClients.opentelemetry.proto.trace.v1.Span(
                                    trace_id = reinterpret(UInt8, rand(UInt128, 1)),
                                    span_id = reinterpret(UInt8, rand(UInt64, 1)),
                                    trace_state = "",
                                    parent_span_id = [],
                                    name = "JJJJJJJJJJJJJJJJJJ",
                                    kind = 0,
                                    start_time_unix_nano = 0xe3c6aac9ebe5b6f8,
                                    end_time_unix_nano = 0xe3c6aac9ebe5b6f9,
                                ),
                            ],
                            schema_url = "http://localhost:7777",
                        ),
                    ],
                    schema_url = "http://localhost:8888",
                ),
            ],
        )

    trace_client = OpentelemetryClients.TraceServiceBlockingClient("http://localhost:4317")

    # OpentelemetryClients.Export(client, req)

    metrics_req =
        OpentelemetryClients.opentelemetry.proto.collector.metrics.v1.ExportMetricsServiceRequest(
            resource_metrics = [
                OpentelemetryClients.opentelemetry.proto.metrics.v1.ResourceMetrics(
                    resource = OpentelemetryClients.opentelemetry.proto.resource.v1.Resource(
                        dropped_attributes_count = 0,
                    ),
                    instrumentation_library_metrics = [
                        OpentelemetryClients.opentelemetry.proto.metrics.v1.InstrumentationLibraryMetrics(
                            instrumentation_library = OpentelemetryClients.opentelemetry.proto.common.v1.InstrumentationLibrary(;
                                name = "julia",
                                version = "1.0.0",
                            ),
                            metrics = [
                                OpentelemetryClients.opentelemetry.proto.metrics.v1.Metric(
                                    name = "test_metric",
                                    description = "just for test",
                                    unit = "",
                                    gauge = OpentelemetryClients.opentelemetry.proto.metrics.v1.Gauge(
                                        data_points = [
                                            OpentelemetryClients.opentelemetry.proto.metrics.v1.NumberDataPoint(
                                                start_time_unix_nano = 0xe3c6aac9ebe5b6f8,
                                                time_unix_nano = 0xe3c6aac9ebe5b6f9,
                                                as_int = 123,
                                            ),
                                        ],
                                    ),
                                ),
                            ],
                            schema_url = "http://localhost:7777",
                        ),
                    ],
                    schema_url = "http://localhost:8888",
                ),
            ],
        )

    log_req =
        OpentelemetryClients.opentelemetry.proto.collector.logs.v1.ExportLogsServiceRequest(
            resource_logs = [
                OpentelemetryClients.opentelemetry.proto.logs.v1.ResourceLogs(
                    resource = OpentelemetryClients.opentelemetry.proto.resource.v1.Resource(
                        dropped_attributes_count = 0,
                    ),
                    instrumentation_library_logs = [
                        OpentelemetryClients.opentelemetry.proto.logs.v1.InstrumentationLibraryLogs(
                            instrumentation_library = OpentelemetryClients.opentelemetry.proto.common.v1.InstrumentationLibrary(;
                                name = "julia",
                                version = "1.0.0",
                            ),
                            log_records = [
                                OpentelemetryClients.opentelemetry.proto.logs.v1.LogRecord(
                                    time_unix_nano = 0xe3c6aac9ebe5b6f9,
                                    observed_time_unix_nano = 0xe3c6aac9ebe5b6f8,
                                    severity_number = OpentelemetryClients.opentelemetry.proto.logs.v1.SeverityNumber.SEVERITY_NUMBER_ERROR,
                                    severity_text = "ERROR",
                                    body = OpentelemetryClients.opentelemetry.proto.common.v1.AnyValue(
                                        string_value = "Oh no!",
                                    ),
                                    attributes = [],
                                    dropped_attributes_count = 0,
                                    trace_id = reinterpret(UInt8, rand(UInt128, 1)),
                                    span_id = reinterpret(UInt8, rand(UInt64, 1)),
                                ),
                            ],
                            schema_url = "http://localhost:7777",
                        ),
                    ],
                    schema_url = "http://localhost:8888",
                ),
            ],
        )
end


text:
function gsd_utility(A::Vector{Int}, u::Matrix{T}, pd::Array{T,3}, n, param; tau=0.0, h=0.001) where T <: Real
    obj = 0.0
    B, J, K = size(pd)

    for b in 1:B
        cost = n[b,end]
        reject = false
        for j in 1:J
            cost = n[b,j]
            for k in A
                if pd[b,j,k] > u[j,k]
                    reject = true
                    break
                end
            end
            if reject
                break
            end
        end
        obj -= cost
    end

    obj / B
end

function sponsor_utility(A::Vector{Int}, u::Matrix{T}, pd::Array{T,3}, n, param::Matrix{T}; tau=0.0, h=0.001) where T <: Real
    obj = 0.0
    B, J, K = size(pd)

    for b in 1:B
        pp = [1 - param[b,end], param[b,end]]
        cost = h*n[b,end]
        reject = false
        for j in 1:J
            cost = h*n[b,j]
            for k in A
                if pd[b,j,k] > u[j,k]
                    obj += pp[k] ## reward for reject
                    reject = true
                end
            end
            if reject
                break
            end
        end
        obj -= cost
    end

    obj / B
end

function public_utility(A::Vector{Int}, u::Matrix{T}, pd::Array{T,3}, n, param::Matrix{T}; tau=0.0, h=0.001) where T <: Real
    obj = 0.0
    B, J, K = size(pd)

    for b in 1:B
        pp = [1 - param[b,end], param[b,end]]
        cost = h*n[b,end]
        reject = false
        for j in 1:J
            cost = h*n[b,j]
            for k in A
                if pd[b,j,k] > u[j,k]
                    if param[b,k] < tau
                        obj += pp[k] ## reward for TP
                    end
                    reject = true
                end
            end

            if reject
                break
            end
        end

        ## reward for TN selection
        for k in 1:K
            if (k ∉ A) & (param[b,k] >= tau)
                obj += pp[k]
            end
        end

        obj -= cost
    end

    obj / B
end

function alpha(pd::Array{T,3}, upper::Matrix{T}; A=1:size(pd,3)) where T <: Real
    B, J, K = size(pd)
    reject = 0

    @inbounds for b in 1:B
        stop = false
        @inbounds for j in 1:J
            @inbounds for k in A
                if pd[b,j,k] > upper[j,k]
                    reject += 1
                    stop = true
                    break
                end
            end
            if stop
                break
            end
        end
    end

    reject / B
end

function objective(A::Vector{Int}, x::Matrix, type1, type2, fsim, fsim0, fsim1; h=0.00025, utility=gsd_utility)
    ## actual type I and type II error for intersection hypothesis
    real_t1 = alpha(fsim0.post_pd, x; A=A)
    real_t2 = 1-alpha(fsim1.post_pd, x; A=A)

    if (real_t1 <= type1) & (real_t2 <= type2)
        utility(A, x, fsim.post_pd, fsim.n, fsim.param; h=h), real_t1, real_t2
    else
        -Inf, real_t1, real_t2
    end
end


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

function part1(filename::AbstractString)
    registers = Dict{AbstractString,Int}()
    for line in eachline(filename)
        target, instr, val, _, lc, cond, rc = split(line)
        if eval(Meta.parse(cond))(get(registers, lc, 0), parse(Int, rc))
            registers[target] = (instr == "inc" ? (+) : (-))(get(registers, target, 0), parse(Int, val))
        end
    end
    maximum(values(registers))
end

@assert part1("example.txt") == 1
part1("input.txt")

function part2(filename::AbstractString)
    registers = Dict{AbstractString,Int}()
    biggest = 0
    for line in eachline(filename)
        target, instr, val, _, lc, cond, rc = split(line)
        if eval(Meta.parse(cond))(get(registers, lc, 0), parse(Int, rc))
            registers[target] = (instr == "inc" ? (+) : (-))(get(registers, target, 0), parse(Int, val))
            biggest = max(registers[target], biggest)
        end
    end
    biggest
end

@assert part2("example.txt") == 10
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using Test, Arrow, Tables, Dates, PooledArrays

include(joinpath(dirname(pathof(Arrow)), "../test/testtables.jl"))
include(joinpath(dirname(pathof(Arrow)), "../test/integrationtest.jl"))

@testset "Arrow" begin

@testset "table roundtrips" begin

for case in testtables
    testtable(case...)
end

end # @testset "table roundtrips"

@testset "arrow json integration tests" begin

for file in readdir(joinpath(dirname(pathof(Arrow)), "../test/arrowjson"))
    jsonfile = joinpath(joinpath(dirname(pathof(Arrow)), "../test/arrowjson"), file)
    println("integration test for $jsonfile")
    df = ArrowJSON.parsefile(jsonfile);
    io = IOBuffer()
    Arrow.write(io, df)
    seekstart(io)
    tbl = Arrow.Table(io; convert=false);
    @test isequal(df, tbl)
end

end # @testset "arrow json integration tests"

@testset "misc" begin

# multiple record batches
t = Tables.partitioner(((col1=Union{Int64, Missing}[1,2,3,4,5,6,7,8,9,missing],), (col1=Union{Int64, Missing}[1,2,3,4,5,6,7,8,9,missing],)))
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == 1
@test isequal(tt.col1, vcat([1,2,3,4,5,6,7,8,9,missing], [1,2,3,4,5,6,7,8,9,missing]))
@test eltype(tt.col1) === Union{Int64, Missing}

# Arrow.Stream
seekstart(io)
str = Arrow.Stream(io)
state = iterate(str)
@test state !== nothing
tt, st = state
@test length(tt) == 1
@test isequal(tt.col1, [1,2,3,4,5,6,7,8,9,missing])

state = iterate(str, st)
@test state !== nothing
tt, st = state
@test length(tt) == 1
@test isequal(tt.col1, [1,2,3,4,5,6,7,8,9,missing])

@test iterate(str, st) === nothing

t = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
meta = Dict("key1" => "value1", "key2" => "value2")
Arrow.setmetadata!(t, meta)
meta2 = Dict("colkey1" => "colvalue1", "colkey2" => "colvalue2")
Arrow.setmetadata!(t.col1, meta2)
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == length(t)
@test tt.col1 == t.col1
@test eltype(tt.col1) === Int64
@test Arrow.getmetadata(tt) == meta
@test Arrow.getmetadata(tt.col1) == meta2

end # @testset "misc"

end

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
module Fortran90Namelists

include("FortranToJulia.jl")
include("JuliaToFortran.jl")  # The order has to be like this!

include("Tokenize.jl")

end # module


text:
struct Metropolis
    steps::Int
    perturb

    Metropolis(steps::Int, perturb) = new(steps, perturb)
    Metropolis(steps::Int) = new(steps, genericperturb)
end

# TODO: this is bad
function genericperturb(node::Node)
    if node.distribution isa Normal
        node.value + std(node.distribution) * randn()
    elseif node.distribution isa Beta
        node.value + .1 * randn()
    elseif node.distribution isa BiMix
        node.value + randn()
    else
        rand(node.distribution)
    end
end

function run(config::Metropolis, model::Model)
    runmodel = compile(model)

    observations = Dict{Symbol, Array}()
    loglikelihoods = Array{Float64}(config.steps)
    foreach(var -> observations[var] = Array{Any}(config.steps), model.observed)

    trace = Base.invokelatest(runmodel, nothing, observations, 1)
    lhood = loglikelihood(trace)
    loglikelihoods[1] = lhood

    for i = 2:config.steps
        proposal = propose(trace, config.perturb)
        proposal = Base.invokelatest(runmodel, proposal, observations, i)
        newlhood = loglikelihood(proposal)

        # NOTE: this is assuming that the proposal kernel for each var is symmetrical, bad
        if log(rand()) <= newlhood - lhood
            trace = proposal
            lhood = newlhood
        else
            # overwrite the observations since proposal was rejected
            foreach(var -> observations[var][i] = observations[var][i-1], model.observed)
        end
        loglikelihoods[i] = lhood
    end

    observations, loglikelihoods
end


text:
module XYZFiles

using AtomBase

function read(io::IOStream, properties::Array{T, 1}) where T<:AbstractProperties
    allAtoms = Array{Atoms, 1}(undef, 0)
    while (line = readline(io)) != ""
        atoms = Atoms()
        number = parse(Int64, split(line)[1])

        line = readline(io)
        words = split(line)
        for prop in properties
            reProp = Regex(prop.name*"=.*")
            matchProp(w) = !(match(reProp, w) === nothing)
            if (ind = findfirst(matchProp.(words))) !== nothing
                val = parse(eltype(prop), split(words[ind], "=")[2])
                prop(val)
                addProperty!(atoms, deepcopy(prop))
            end
        end

        for i in 1:number
            line = readline(io)
            words = split(line)
            atom = Atom(String(words[1]), [parse(Float64,words[2]), parse(Float64,words[3]), parse(Float64,words[4])])
            addAtom!(atoms, atom)
        end

        push!(allAtoms, atoms)
    end

    return allAtoms
end

function read(file::String, properties::Array{T, 1}) where T<:AbstractProperties
    io = open(file, "r")
    structures = read(io, properties)
    close(io)
    return structures
end

function read(io::IOStream)
    properties = Array{Properties, 1}(undef, 0)
    read(io, properties)
end

function read(file::String)
    properties = Array{Properties, 1}(undef, 0)
    read(file, properties)
end

function read(io::IOStream, propertiesString::Array{String, 1})
    properties = Array{Properties, 1}(undef, 0)
    for propString in propertiesString
        prop = Properties{Float64}(propString, 0.0)
        push!(properties, prop)
    end

    read(io, properties)
end

function read(file::String, propertiesString::Array{String, 1})
    properties = Array{Properties, 1}(undef, 0)
    for propString in propertiesString
        prop = Properties{Float64}(propString, 0.0)
        push!(properties, prop)
    end

    io = open(file, "r")
    structures = read(io, properties)
    close(io)
    return structures
end

function write(io::IOStream, allAtoms::Array{Atoms, 1})
    for atoms in allAtoms
        println(io, atoms.number)

        for prop in atoms.properties
            print(io, prop)
        end
        print(io, "\n")

        for atom in atoms.atom
            println(io, atom.symbol, "\t", atom.position[1], "\t", atom.position[2], "\t", atom.position[3])
        end
    end
    return nothing
end

function write(file::String, allAtoms::Array{Atoms, 1})
    io = open(file, "w")
    write(io, allAtoms)
    close(io)
    return nothing
end

function write(io::IOStream, atoms::Atoms)
    write(io, [atoms, ])
end

function write(file::String, atoms::Atoms)
    write(file, [atoms,])
end

end

max_stars_repo_name: yaozhenghangma/XYZ

max_issues_repo_name: yaozhenghangma/XYZ

max_forks_repo_name: yaozhenghangma/XYZ


text:
using Test

# Check to see if an environment variable is defined. If so, save plots there.
outputdir = get(ENV, "JL_PLOT_TEMPDIR", nothing)
if isnothing(outputdir)
    # If there is no environment variable grab a temporary directory and save all plots there by default
    outputdir = mktempdir()
else
    # if it doesn't exist, make the directory
    if !ispath(outputdir)
        mkdir(outputdir)
    end
end

println("Saving plots in test/$outputdir")

include("interpolation.jl") # visibility interpolation

include("ensemble_sampler.jl") # sampler test


text:
using Ekztazy

### +----------------------
### | This is basically just a remake of
### | https://github.com/Kyando2/AniBot/blob/master/cogs/Money.py.
### | The commands are just being 1:1 copied with a few additions
### | that I wanted to do anyway.
### | An non exhaustive list of big differences:
### | ーMoney doesn't exist, it's just calculated based on the posessed items.
### | ーーPaying/Buying doesn't exist. You can trade with items.
### | ーーsome items are more suited to an economy.
### | ーComponents/Slash Commands instead of *disgusting* plain old commands.
### | ーRemove the absolutely unnecessary OOP (like ItemHandler and PlayerHandler).
### | ーFix a lot of really just bad code in there (I was like 12 when I wrote it spare me).
### | ーAdd some new features.
### | ーーCharacter, you now get to choose a character which gives you some bonuses.
### | ーーWorld map, I had started that over in the other project, but now its a full fledged map system.
### | ーーEmpires! You get to fight for a country, or even maybe create one!
### | ーーPVP!!!
### | ーーPVE!!!
### | ーThis is made in Julia, not Python. Python is disgusting.
### +--
### | The text here definitely needs some help.
### +--
### | Maybe perform registration check for all commands.
### +--
### | Maybe handlers should be moved to a sub folder,
### | so each series of command->component has its own folder.
### +--
### | Maybe make a function in separate modules for each handler series.
### +--
### | The getClient() function is just to prevent having to do global
### | on each handler, instead you can just have a default argument
### | which shouldn't mess with Ekztazy, unless you try to do this
### | with the new mode of handlers (which would then mess with Ekztazy in 0.7,
### | see https://github.com/Humans-of-Julia/Ekztazy.jl/issues/54).
### +--
### | getId() should probably be implemented by Ekztazy
### | (but like better, like actually check if its in a guild or in DMs),
### | IDK who made Ekztazy but like why do they not have that??
### +--
### | This uses Ekztazy >0.5, and probably doesn't work with 0.4 but IDK tbh.
### +----------------------

function getClient()
    global c
    c
end

getId(ctx::Context) = string(ctx.interaction.member.user.id)

function init()
    global c

    g = ENV["TESTGUILD"]

    command!(register_command_handler, c, g, "register", "Register as a new player")

    on_ready!(c) do ctx
        @info "Successfully logged in!"
    end

    start(c)
end

function register_command_handler(ctx, c=getClient())
    global REG_DEF_COMP
    id = getId(ctx)
    t = getUserData(id)
    if t[1] == false
        reply(c, ctx, components=[REG_DEF_COMP], content="Hey, welcome to the world of Elvur!", flags=1<<6)
    else
        reply(c, ctx, content="You are already registered!")
    end
end

function register_component_first_default(ctx, c=getClient())
    global REG_KEN_COMP2
    reply(c, ctx, components=[REG_KEN_COMP2], content="Welcome to the world of Elvur... What character you want?", flags=1<<6)
end

function register_component_second_ken(ctx, c=getClient())
    id = getId(ctx)
    # ken is id 0, make const later
    createUser(id, 0)
    reply(c, ctx, content="So you've chosen **Ken Kaneki**. You can now start playing.", flags=1<<6)
end

c = Client()
REG_DEF_COMP = component!(register_component_first_default, c, "default00"; type=2, style=1, label="Start!")
REG_KEN_COMP2 = component!(register_component_second_ken, c, "ken00"; type=2, style=1, label="Ken Kaneki")


text:
"""
    cuparse(raw::CDataRaw; dim="2Dsym")
    cuparse(path::String)

**Keyword Arguments**
- `mode::String`: "dset", "table" _(, "auto")_

!!! note "COMSOL"
    Latest tested version: 5.3.0.223
"""
function cuparse(raw::CDataRaw; dim="2Dsym", mode="auto")
    mode = mode == "auto" ? _autodetectmode(raw.field) : mode

    # switch between the data mode
    return mode == "dset" ? cuparse_dset(raw, dim) :
           mode == "table" ? cuparse_table(raw, dim) :
           error("Keyword mode with unsupported argument \"$mode\". ")
end
cuparse(path::String; dim="2Dsym") = cuparse(curead(path); dim=dim)

"""
    cuparse_dset(raw::CDataRaw, dim::String)
"""
function cuparse_dset(raw::CDataRaw, dim::String)
    n = String[]; u = String[]; t = Float64[]; v = Dict{String, Any}(); # init vectors

    N, _ = _parsedim(dim) # dimension
    for i in N+1:length(raw.field)
        name, unit, time = _readfieldall(raw.field[i]) # extract data from field
        v[name] = haskey(v, name) ? hcat(v[name], raw.v[raw.field[i]]) : raw.v[raw.field[i]] # values (Dict)
        push!(t, time) # time

        if name ∈ n; continue; end # name MUST be unique
        push!(n, name); push!(u, unit); # name and unit
    end

    t = sort(unique(t)) # time
    p = _parseposition(raw.v, dim) # position

    return CDataDSet(raw.path, raw.meta, n, u, t, p, v)
end

"""
    cuparse_table(raw::CDataRaw, dim::String)
"""
function cuparse_table(raw::CDataRaw, dim::String)
    n = String[]; u = String[]; v = Dict{String, Any}(); # init vectors

    old_keys = collect(keys(raw.v)); # old keys from raw dictionary

    for i in 1:length(raw.field)
        new_key, unit = _separatenameunit(old_keys[i]); # split field
        push!(n, new_key); push!(u, unit); # push name and unit
        v[new_key] = raw.v[old_keys[i]]; # add to Dict
    end

    return CDataTable(raw.path, raw.meta, n, u, v)
end

## INTERNAL FUNCTION ======================================================== ##
function _autodetectmode(field::Vector{String})
    sum([occursin("@", field[i]) ? 1 : 0 for i in length(field)]) > 0 ? "dset" : "table"
end

function _readfieldall(s::String)
    n, u, _, t = split(s, " ") # expected format: `<name> (<unit>) @ t=<time>`
    return n, u[2:end-1], parse(Float64, t[3:end])
end

function _separatenameunit(s::String)
    for i in length(s):-1:1
        if s[i] == '('
            return s[1:i-2], s[i+1:end-1]
        end
    end
end

function _parsedim(dim::String)
    if dim == "2Dsym"
        return 2, ("R", "Z")
    end
end

function _parseposition(v::Dict{String, Any}, dim::String)
    if dim == "2Dsym"
        return CPosition2Dsym(v["R"], v["Z"])
    end
end


text:
module Antiresonant
using Reexport
import Logging: @warn
import Printf: @sprintf
import Luna: Capillary
import Luna.PhysData: c, wlfreq
@reexport using Luna.Modes
import Luna.Modes: AbstractMode, dimlimits, neff, field, Aeff, N, α

struct ZeisbergerMode{mT<:Capillary.MarcatiliMode, LT} <: AbstractMode
    m::mT
    wallthickness::Float64
    loss::LT # Val{true}(), Val{false}() or a number (scaling factor)
end

"""
    ZeisbergerMode(args...; wallthickness, kwargs...)

Create a capillary-like mode with the effective index given by eq. (15) in [1].

`wallthickness` (mandatory kwarg) sets the thickness of the anti-resonant struts and
`loss` (optional, defaults to `true`) can be either a `Bool` (to switch on/off loss
completely) or a `Real` (to up/down-scale the loss given by the model).
 Other kwargs are passed on to the constructor of a [`Capillary.MarcatiliMode`](@ref).

[1] Zeisberger, M., Schmidt, M.A. Analytic model for the complex effective index of the
leaky modes of tube-type anti-resonant hollow core fibers. Sci Rep 7, 11761 (2017).
https://doi.org/10.1038/s41598-017-12234-5
"""
function ZeisbergerMode(args...; wallthickness, loss=true, kwargs...)
    return ZeisbergerMode(Capillary.MarcatiliMode(args...; kwargs...),
                          wallthickness, wraptype(loss))
end

function ZeisbergerMode(m::Capillary.MarcatiliMode; wallthickness, loss=true)
    ZeisbergerMode(m, wallthickness, wraptype(loss))
end

wraptype(loss::Bool) = Val(loss)
wraptype(loss::Real) = loss
wraptype(loss) = throw(
    ArgumentError("loss has to be a Bool or Real, not $(typeof(loss))"))

# Effective index is given by eq (15) in [1]
neff(m::ZeisbergerMode, ω; z=0) = _neff(m.m, ω, m.wallthickness, m.loss; z=z)

# All other mode properties are identical to a MarcatiliMode
for fun in (:Aeff, :field, :N, :dimlimits)
    @eval ($fun)(m::ZeisbergerMode, args...; kwargs...) = ($fun)(m.m, args...; kwargs...)
end

function _neff(m::Capillary.MarcatiliMode, ω, wallthickness, loss; z=0)
    nco = m.coren(ω, z=z)
    ncl = m.cladn(ω, z=z)
    ϵ = ncl^2 / nco^2
    k0 = ω / c
    ka = k0*nco
    ϕ = k0*wallthickness*sqrt(ncl^2 - nco^2)
    σ = 1/(ka*Capillary.radius(m, z))
    if (m.kind == :TE) || (m.kind == :TM)
        A = m.unm^2/2;
        if m.kind == :TE
            B = A/(sqrt(ϵ - 1))/tan(ϕ);
            C = m.unm^4/8 + 2*m.unm^2/(ϵ - 1)/tan(ϕ)^2
            D = m.unm^3*(1+1/tan(ϕ)^2)/(ϵ - 1)
        elseif m.kind == :TM
            B = A*ϵ/(sqrt(ϵ - 1))/tan(ϕ)
            C = m.unm^4/8 + 2*m.unm^2*ϵ^2/(ϵ - 1)/tan(ϕ)^2
            D = m.unm^3*ϵ^2*(1+1/tan(ϕ)^2)/(ϵ - 1)
        end
        return __neff(A, B, C, D, σ, nco, loss)
    else
        if m.kind == :EH
            s = 1
        elseif m.kind == :HE
            s = -1
        end
        A = m.unm^2/2
        B = m.unm^2*(ϵ+1)/(sqrt(ϵ-1)*tan(ϕ))
        C = (m.unm^4/8
            + m.unm^2*m.m*s/2
            + (m.unm^2/4*(2+m.m*s)*(ϵ+1)^2/(ϵ-1)
            - m.unm^4*(ϵ-1)/(8*m.m))*1/tan(ϕ)^2)
        D = m.unm^3/2 * (ϵ^2+1)/(ϵ-1) * (1 + 1/tan(ϕ)^2)
        __neff(A, B, C, D, σ, nco, loss)
    end
end

#= If ncl and nco are real, the type of neff depends on whether loss is included.
    Despatching on the type makes sure the type can be inferred by the compiler =#
__neff(A, B, C, D, σ, nco, loss::Val{true}) = nco*(1 - A*σ^2 - B*σ^3 - C*σ^4 + 1im*D*σ^4)
__neff(A, B, C, D, σ, nco, loss::Val{false}) = real(nco*(1 - A*σ^2 - B*σ^3 - C*σ^4))
__neff(A, B, C, D, σ, nco, loss::Number) = nco*(1 - A*σ^2 - B*σ^3 - C*σ^4 + 1im*loss*D*σ^4)

struct VincettiMode{mT<:Capillary.MarcatiliMode, LT} <: AbstractMode
    m::mT
    t::Float64 # tube wall thickness
    r_ext::Float64 # tube radius
    Ntubes::Int # number of tubes
    cladn::Float64 # cladding refractive index (constant, defaults to 1.45)
    Nterms::Int # number of terms (cladding modes) to include in the sum
    loss::LT # Val{true}(), Val{false}() or a number (scaling factor)
end

"""
    VincettiMode(a, args...; wallthickness, tube_radius, Ntubes, cladn, Nterms,
                             loss=true, kwargs...)

Create a mode with Marcatili-like mode fields but loss, dispersion and effective area given
by the semi-empirical model developed by Vincetti et al. in refs [1-3]. Arguments are identical
to `Capillary.MarcatiliMode` but with the following additions/changes as keyword arguments:

# Mandatory keyword arguments
- `wallthickness` : wall thickness of the resonators (cladding capillaries)
- `tube_radius` : radius of the resonators
- `Ntubes` : number of resonators

# Optional keyword arguments
- `cladn` : refractive index of the resonators. This must be a constant for the model.
            Defaults to 1.45
- `Nterms` : number of resonator dielectric modes to include in the model. Defaults to 8.
- `loss` : can be `true` or `false` to switch loss on/off, or a `Real` to scale the loss.

To specify the gap between resonators, calculate the core radius with [`getRco(r_ext, N, δ)`](@ref)
or calculate the external radius of the resonators with [`getr_ext(Rco, N, δ)`](@ref).

# References

[1] L. Vincetti
Empirical formulas for calculating loss in hollow core tube lattice fibers,
Opt. Express, OE, vol. 24, no. 10, pp. 10313-10325, May 2016, doi: 10.1364/OE.24.010313.

[2] L. Vincetti and L. Rosa
A simple analytical model for confinement loss estimation in hollow-core Tube Lattice Fibers
Opt. Express, OE, vol. 27, no. 4, pp. 5230-5237, Feb. 2019, doi: 10.1364/OE.27.005230.

[3]L. Rosa, F. Melli, and L. Vincetti
Analytical Formulas for Dispersion and Effective Area in Hollow-Core Tube Lattice Fibers
Fibers, vol. 9, no. 10, Art. no. 10, Oct. 2021, doi: 10.3390/fib9100058.
"""
function VincettiMode(Rco, args...; wallthickness, tube_radius, Ntubes,
                                    cladn=1.45, Nterms=8, loss=true, kwargs...)
    if getδ(Rco, tube_radius, Ntubes) < 0
        @warn("the given fibre parameters correspond to a negative gap between resonators")
    end
    return VincettiMode(Capillary.MarcatiliMode(Rco, args...; kwargs...),
                          wallthickness, tube_radius, Ntubes, cladn, Nterms, wraptype(loss))
end

# create complex effective index
neff(m::VincettiMode, ω; z=0) = neff_real(m, ω) + 1im*c/ω*α(m, ω; z)

α(m::VincettiMode{mT, Val{true}}, ω; z=0) where mT = log(10)/10 * CL(m, ω; z)
α(m::VincettiMode{mT, Val{false}}, ω; z=0) where mT = zero(ω)
α(m::VincettiMode{mT, <:Number}, ω; z=0) where mT = m.loss * log(10)/10 * CL(m, ω; z)

# All other mode properties are identical to a MarcatiliMode
for fun in (:field, :N, :dimlimits)
    @eval ($fun)(m::VincettiMode, args...; kwargs...) = ($fun)(m.m, args...; kwargs...)
end

function CL(λ, Rco, t, r_ext, N; cladn=1.45, Nterms=8)
    # eq. (6) of [2]
    # confinement loss in dB/m
    F = normfreq(λ, t, cladn)
    pvs = p_ν_sum(F, t, r_ext, cladn, Nterms)
    clm = CLmin(Rco, λ, t, r_ext, cladn)
    return clm*pvs
end

CL(m::VincettiMode, ω; z=0) = CL(wlfreq(ω), m.m.a, m.t, m.r_ext, m.Ntubes;
                                 cladn=m.cladn, Nterms=m.Nterms)

function FcHE(μ::Integer, ν::Integer, t, r_ext, n=1.45)
    # eq. (4) in [2]
    if ν == 1
        return (abs(0.21 + 0.175μ - 0.1/(μ-0.35)^2) * (t/r_ext)^(0.55+5e-3*sqrt(n^4-1))
                + 0.04*sqrt(μ)*t/r_ext)
    else
        return 0.3/n^0.3 * (2/ν)^1.2*abs(μ-0.8)*(t/r_ext) + ν - 1
    end
end

function FcEH(μ::Integer, ν::Integer, t, r_ext, n=1.45)
    # eq. (5) in [2]
    if ν == 1
        return (0.73 + 0.57*(μ^0.8 + 1.5)/4 - 0.04/(μ-0.35))*(t/r_ext)^(0.5 - (n-1)/(10*(μ+0.5)^0.1))
    else
        tmp = (11.5/(ν^1.2*(7.75-ν))*(0.34+μ/4*(n/1.2)^1.15)/(μ+0.2/n)^0.15
               * (t/r_ext)^(0.75 + 0.06/n^1.15 + 0.1*sqrt(1.44/n)*(ν-2)))
        return tmp + ν -1
    end
end

function CLmin(Rco, λ, t, r_ext, n)
    # eq. (3) of [1]
    3e-4 * λ^4.5/Rco^4 * (1-t/r_ext)^-12 * sqrt(n^2-1)/(t*sqrt(r_ext)) * exp(2λ/(r_ext*(n^2-1)))
end

γloss = 3e-3
L(F) = γloss^2/(γloss^2 + F^2) # eq. (2) of [2]

A(μ) = 2e3 * exp(-0.05*abs(μ-1)^2.6) # eq. (3) of [2]

function p_ν(F, ν, t, r_ext, n=1.45, Nterms=8)
    # eq. (7) of [2]
    out = A(1) * (L(F-FcHE(1, ν, t, r_ext, n)) + L(F - FcEH(1, ν, t, r_ext, n)))
    for μ in 2:Nterms
        out += A(μ) * (L(F-FcHE(μ, ν, t, r_ext, n)) + L(F - FcEH(μ, ν, t, r_ext, n)))
    end
    out
end

function p_ν_sum(F, t, r_ext, n=1.45, Nterms=8)
    out = p_ν(F, 1, t, r_ext, n, Nterms)
    for ν in 2:Nterms
        out += p_ν(F, ν, t, r_ext, n, Nterms)
    end
    out
end

normfreq(λ, t, n, nco=1) = 2t/λ*sqrt(n^2-nco) # eq. (2) of [1]

function Rco_eff(λ, Rco, t, r_ext, N, n=1.45)
    # eq. (10) in [3]
    δ = getδ(Rco, r_ext, N)
    F = normfreq(λ, t, n)
    t1 = 1.027 + 1e-3*(F + 2/F^4)
    # corrected a probable typo here: |
    #                                 V
    t2 = sqrt(Rco^2 + N/π*3/64* r_ext^2 *(1+(3+20λ/Rco)*δ/r_ext))
    # in the paper there is no square, but this would be dimensionally incorrect
    # and eq. (9) has the square
    return t1*t2
end

Rco_eff(m::VincettiMode, ω) = Rco_eff(wlfreq(ω), m.m.a, m.t, m.r_ext, m.Ntubes, m.cladn)

"""
    getδ(Rco, r_ext, N)

Calculate the gap between resonators in a single-ring antiresonant PCF with inscribed core
radius `Rco` for `N` resonators with external radius `r_ext`.
"""
getδ(Rco, r_ext, N) = 2*(sin(π/N)*(Rco + r_ext) - r_ext) # from eq. (1) in [1]

"""
    getRco(r_ext, N, δ)

Calculate the inscribed core radius of a single-ring antiresonant PCF with `N` resonators
with external radius `r_ext` and a gap between resonators of `δ`.
"""
function getRco(r_ext, N, δ)
    # eq. (1) of [1]
    k = 1 + δ/2r_ext
    Rco = r_ext * (k/sin(π/N) - 1)
    return Rco
end

"""
    getr_ext(Rco, N, δ)

Calculate the external radius of the resonators for a single-ring antiresonant PCF with
core radius `Rco`, `N` resonators and a gap between resonators of `δ`.
"""
getr_ext(Rco, N, δ) = (δ/2 - Rco*sin(π/N))/(sin(π/N) - 1)

γdisp = 3e-2
Li(F, F_0) = (F_0^2 - F^2)/((F^2 - F_0^2)^2 + (γdisp*F)^2) # eq. (2) in [3]

function νsum(F, t, r_ext, n=1.45, Nterms=8)
    # eq. (6) of [3]
    # here only μ=1 is considered
    # ν=1 term:
    out = Li(F, FcHE(1, 1, t, r_ext, n)) + Li(F, FcEH(1, 1, t, r_ext, n))
    for ν in 2:Nterms
        out += Li(F, FcHE(1, ν, t, r_ext, n)) + Li(F, FcEH(1, ν, t, r_ext, n))
    end
    return out * A(1)
end

function Δneff(λ, Rco, t, r_ext, n=1.45, Nterms=8)
    # eq. (6) of [3]
    F = normfreq(λ, t, n)
    ρ = 1 - t/r_ext
    return 4.5e-7/ρ^4 * (λ/Rco)^2 * νsum(F, t, r_ext, n, Nterms)
end

Δneff(m::VincettiMode, ω) = Δneff(wlfreq(ω), m.m.a, m.t, m.r_ext, m.cladn, m.Nterms)

function neff_real(m::VincettiMode, ω; z=0)
    # eq. (21) of [3]
    ng = m.m.coren(ω; z) # gas index
    return (ng
            - 1/2 * (m.m.unm*c/(ω*ng*Rco_eff(m, ω)))^2
            + Δneff(m, ω))
end

Aeff(m::VincettiMode, ω; z=0) = 0.48/8π * (m.m.unm*wlfreq(ω))^2/(m.m.coren(ω; z)-neff_real(m, ω))

end # module


text:
module Utils
    include("./src/Chronometer.jl")
    import .Chronometer

    include("./src/FileReader.jl")
    import .FileReader

    include("./src/ConfigReader.jl")
    import .ConfigReader

    include("./src/Getters.jl")
    import .Getters

    include("./src/Declares.jl")
    import .Declares
end


text:
using Documenter
#using CachedArrays

makedocs(
    modules     = [Documenter],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename    = "CachedArrays",
    doctest     = true,
    pages       = Any[
        "Getting Started"   => "index.md",
        "MemoryAllocation" => Any[
            "blocks.md",
        ],
    ],
)

deploydocs(
    repo        = "github.com/darchr/CachedArrays.jl.git",
    target      = "build",
)


text:
# Pauli Groups
export Pauli

"""
Pauli Matrix
"""
module Pauli
export sigmax, sigmay, sigmaz, sigmai
using Compat.SparseArrays

const sigmax = sparse([0 1;1 0])
const sigmay = sparse([0 -im;im 0])
const sigmaz = sparse([1 0;0 -1])
const sigmai = sparse([1 0;0 1])

export σ, σ₀, σ₁, σ₂, σ₃

const σ₀ = sigmai
const σ₁ = sigmax
const σ₂ = sigmay
const σ₃ = sigmaz

const σ = [sigmai, sigmax, sigmay, sigmaz]

# Pauli Vec

export PauliVector, S, *, ⋅

"""
    PauliVector

PauliVector ``(\sigma_x, \sigma_y, \sigma_z)``
"""
struct PauliVector
end

const S = PauliVector()

import Base: *
*(A::PauliVector, B::PauliVector) = kron(sigmax, sigmax) + kron(sigmay, sigmay) + kron(sigmaz, sigmaz)
⋅(A::PauliVector, B::PauliVector) = A * B

import Base: getindex

function getindex(pv::PauliVector, index::Integer)
    if index == 1
        return sigmax
    elseif index == 2
        return sigmay
    elseif index == 3
        return sigmaz
    else
        throw(BoundsError(pv, index))
    end
end

end