text:
Network = NetLib.SimpleNet

netparams = NetLib.SimpleNetHP(
  width=100,
  depth_common=4,
  use_batch_norm=false)

self_play = SelfPlayParams(
  sim=SimParams(
    num_games=1000,
    num_workers=4,
    use_gpu=false,
    reset_every=16,
    flip_probability=0.,
    alternate_colors=false),
  mcts=MctsParams(
    num_iters_per_turn=50,
    cpuct=1.0,
    temperature=ConstSchedule(0.),
    dirichlet_noise_ϵ=0.,
    dirichlet_noise_α=1.))

arena = ArenaParams(
  sim=SimParams(
    num_games=100,
    num_workers=10,
    use_gpu=false,
    reset_every=1,
    flip_probability=0.,
    alternate_colors=true),
  mcts = self_play.mcts,
  update_threshold=0.00)

learning = LearningParams(
  use_gpu=false,
  use_position_averaging=false,
  samples_weighing_policy=CONSTANT_WEIGHT,
  rewards_renormalization=10,
  l2_regularization=1e-4,
  optimiser=Adam(lr=5e-3),
  batch_size=64,
  loss_computation_batch_size=2048,
  nonvalidity_penalty=1.,
  min_checkpoints_per_epoch=1,
  max_batches_per_checkpoint=5_000,
  num_checkpoints=1)

params = Params(
  arena=arena,
  self_play=self_play,
  learning=learning,
  num_iters=5,
  memory_analysis=nothing,
  ternary_rewards=false,
  use_symmetries=false,
  mem_buffer_size=PLSchedule(80_000))

benchmark_sim = SimParams(
  arena.sim;
  num_games=500,
  num_workers=10)

benchmark = [
  Benchmark.Single(
    Benchmark.Full(self_play.mcts),
    benchmark_sim),
  Benchmark.Single(
    Benchmark.NetworkOnly(),
    benchmark_sim)]

experiment = Experiment(
  "grid-world", GameSpec(), params, Network, netparams, benchmark)


text:
using DiffEqOperators
# using Plots

n = 10000
@show n
x = range(0.0; length = n, stop = 2π)
dx = x[2]-x[1]
y = exp.(π*x)
y_ = y[2:(end-1)]

for dor in 1:4, aor in 2:6

    D1 = CenteredDifference(dor,aor,dx,length(x))

    #test result
    @show dor
    @show aor
    #take derivatives
    err_abs = abs.(D1*y .- (π^dor)*y_)
    err_percent = 100*err_abs./abs.(((π^dor)*y_))
    max_error = maximum(err_percent) # test operator with known derivative of exp(kx)
    avg_error = sum(err_percent)/length(err_percent)

    @show max_error
    @show avg_error
    #plot(x[2:(end-1)], err_percent, title="Percentage Error, n=$n aor = $aor, dor = $dor")

    #test result
end

ext: jl


text:
## main file for Markov chain related structs and functions

"Check if all the probabilities are greater equal than 0."
function check_prob(mat_prob::Union{Array{Float64,2}, Array{Int64,1}})
    if all(mat_prob .<= 0)
        throw("All probabilities must be greater equal than 0.")
    end
end

"Check if all the summations of rows in stpm are 1."
function check_prob_sum(stpm::Array{Float64,2})
    n_row = size(stpm)[1]
    if all(stpm * ones(n_row, 1) .== 1)
        throw("The sum of probabilities in $i-th row doesn't equal to 1.")
    end
end

"Calculate limiting probability distribution of a Markov chain family."
function cal_vec_lpd(stpm::Array{Float64,2})

    hcat([1; 1; 1])
    return vec_lpd
end

"Define the non-mutable struct for Markov chain family."
struct MarkovChainFamily
    stpm::Array{Float64,2}  # [Stationary Transition Probability Matrix]
    vec_lpd::Array{Float64,2}  # [Limiting Probability Distribution]

    function MarkovChainFamily(stpm::Array{Float64,2})
        check_prob(stpm)
        check_prob_sum(stpm)

        vec_lpd = cal_vec_lpd(stpm)

        new(stpm, vec_lpd)
    end
end


text:
###
### Indexing
###
###
### Indexing methods for mutable biological sequences.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

#=
Indexing into BioSequences works as follows:

A `Base.getindex` method is provided for `BioSequence`.
It's job is to, do boundschecking. And to then call `inbounds_getindex`.

It is the task of `inbounds_getindex` to extract from a sequence, the
appropriate bits, and pass them to an appropriate `decode` method.
=#

@inline function symbols_per_data_element(seq::BioSequence)
    return div(8 * sizeof(encoded_data_eltype(seq)), bits_per_symbol(seq))
end

@inline function bitindex(seq::BioSequence, i::Integer)
    return bitindex(BitsPerSymbol(seq), encoded_data_eltype(seq), i)
end

@inline function firstbitindex(seq::BioSequence)
    return bitindex(seq, 1)
end

@inline function lastbitindex(seq::BioSequence)
    return bitindex(seq, lastindex(seq))
end

@inline function bindata_mask(seq::BioSequence)
    return bitmask(Alphabet(seq))
end

function eachbitindex(seq::BioSequence, from = firstindex(seq), to = lastindex(seq))
    return bitindex(seq, from):bits_per_symbol(seq):bitindex(seq, to)
end

"""
    inbounds_getindex(seq::BioSequence, i::BitIndex)

An inbounds_getindex method defined by default for all BioSequences.

It is the job of this method to take a sequence and a bitindex, extract the
encoded bits and decode them into a symbol.

!!! note
    This method might be overloaded for specialised subtypes of BioSequence.
"""
@inline function inbounds_getindex(seq::BioSequence, bidx::BitIndex)
    encoded_symbol = extract_encoded_element(bidx, encoded_data(seq))
    return decode(Alphabet(seq), encoded_symbol)
end

"""
    inbounds_getindex(seq::BioSequence, i::Integer)

An inbounds_getindex method defined by default for all BioSequences.

It is the job of this method to convert the integer index to a `BitIndex`, and
call the method of inbounds_getindex, for the sequence and `BitIndex`
combination.

!!! note
    This method might be overloaded for specialised subtypes of BioSequence.
"""
@inline function inbounds_getindex(seq::BioSequence, i::Integer)
    bidx = bitindex(seq, i)
    return inbounds_getindex(seq::BioSequence, bidx)
end

function checkdimension(from::Integer, to::Integer)
    if from == to
        return true
    end
    throw(
        DimensionMismatch(
            string("attempt to assign ", from, " elements to ", to, " elements")
        )
    )
end

function checkdimension(seq::BioSequence, locs::AbstractVector)
    return checkdimension(length(seq), length(locs))
end

function checkdimension(seq::BioSequence, locs::AbstractVector{Bool})
    return checkdimension(length(seq), sum(locs))
end

@inline function unsafe_setindex!(s::BioSequence, v, i::Integer)
    return unsafe_setindex!(s, v, bitindex(s, i))
end

@inline function unsafe_setindex!(s::BioSequence, v, i::BitIndex)
    v_ = convert(eltype(s), v)
    encoded = encode(Alphabet(s), v_)
    return encoded_setindex!(s, encoded, i)
end

###
### Base methods
###

# Provided Base methods for indexing any BioSequence type. May be overloaded if
# needed but it is unlikely.

Base.eltype(::Type{T}) where T <: BioSequence = eltype(Alphabet(T))
Base.eltype(seq::BioSequence) = eltype(Alphabet(seq))
Base.size(seq::BioSequence) = (length(seq),)

Base.firstindex(seq::BioSequence) = 1
Base.lastindex(seq::BioSequence) = length(seq)
Base.eachindex(seq::BioSequence) = Base.OneTo(lastindex(seq))

# Bounds checking...
@inline function Base.checkbounds(seq::BioSequence, i::Integer)
    if 1 ≤ i ≤ lastindex(seq)
        return true
    end
    throw(BoundsError(seq, i))
end

@inline function Base.checkbounds(seq::BioSequence, locs::AbstractVector{Bool})
    if length(seq) == length(locs)
        return true
    end
    throw(BoundsError(seq, locs))
end

@inline function Base.checkbounds(seq::BioSequence, locs::AbstractVector)
    for i in locs
        checkbounds(seq, i)
    end
    return true
end

@inline function Base.checkbounds(seq::BioSequence, range::UnitRange)
    if 1 ≤ range.start && range.stop ≤ lastindex(seq)
        return true
    end
    throw(BoundsError(seq, range))
end

@inline function Base.getindex(seq::BioSequence, i::Integer)
    @boundscheck checkbounds(seq, i)
    return inbounds_getindex(seq, i)
end

@inline function Base.setindex!(s::BioSequence, v, i::Integer)
    @boundscheck checkbounds(s, i)
    return unsafe_setindex!(s, v, i)
end

function Base.setindex!(seq::BioSequence, x, locs::AbstractVector{<:Integer})
    @boundscheck checkbounds(seq, locs)
    return unsafe_setindex!(seq, x, locs)
end

@inline function Base.iterate(seq::BioSequence, i::Int = firstindex(seq))
    if i > lastindex(seq)
        return nothing
    else
        return inbounds_getindex(seq, i), i + 1
    end
end


text:
#Grid setup

type data1d
    x::Array{Float64,1} #cell centers
    xi::Array{Float64,1} #cell LEFT interfaces

    rho::Array{Float64,1}
    rhop::Array{Float64,1}
    rhom::Array{Float64,1}

    velx::Array{Float64,1}
    velxp::Array{Float64,1}
    velxm::Array{Float64,1}

    vely::Array{Float64,1}
    velyp::Array{Float64,1}
    velym::Array{Float64,1}

    eps::Array{Float64,1}
    epsp::Array{Float64,1}
    epsm::Array{Float64,1}

    press::Array{Float64,1}
    pressp::Array{Float64,1}
    pressm::Array{Float64,1}

    q::Array{Float64,2} #conserved quantities
    qp::Array{Float64,2}
    qm::Array{Float64,2}

    n::Int64
    g::Int64 #ghost cells

    function data1d(nzones::Int64)

        x = zeros(nzones)
        xi = zeros(nzones)

        rho = zeros(nzones)
        rhop = zeros(nzones)
        rhom = zeros(nzones)

        velx = zeros(nzones)
        velxp = zeros(nzones)
        velxm = zeros(nzones)

        vely = zeros(nzones)
        velyp = zeros(nzones)
        velym = zeros(nzones)

        eps = zeros(nzones)
        epsp = zeros(nzones)
        epsm = zeros(nzones)

        press = zeros(nzones)
        pressp = zeros(nzones)
        pressm = zeros(nzones)

        q = zeros(nzones, 4)
        qp = zeros(nzones, 4)
        qm = zeros(nzones, 4)
        n = nzones
        g = 3

        new(x,
            xi,
            rho,
            rhop,
            rhom,
            velx,
            velxp,
            velxm,
            vely,
            velyp,
            velym,
            eps,
            epsp,
            epsm,
            press,
            pressp,
            pressm,
            q,
            qp,
            qm,
            n,
            g)
    end
end

type data2d
    x::Array{Float64, 1} #cell centers for X
    y::Array{Float64, 1} #cell centers for Y
    xi::Array{Float64, 1} #cell LEFT interfaces
    yi::Array{Float64, 1} #cell LOWER interfaces

    rho::Array{Float64, 2}
    rhop::Array{Float64, 3}
    rhom::Array{Float64, 3}

    velx::Array{Float64, 2}
    velxp::Array{Float64, 3}
    velxm::Array{Float64, 3}
    vely::Array{Float64, 2}
    velyp::Array{Float64, 3}
    velym::Array{Float64, 3}

    eps::Array{Float64, 2}
    epsp::Array{Float64, 3}
    epsm::Array{Float64, 3}

    press::Array{Float64, 2}
    pressp::Array{Float64, 3}
    pressm::Array{Float64, 3}

    q::Array{Float64, 3} #conserved quantities
    qp::Array{Float64, 4}
    qm::Array{Float64, 4}

    nx::Int64
    ny::Int64
    g::Int64 #ghost cells

    function data2d(nx::Int64, ny::Int64)

        x = zeros(nx)
        y = zeros(ny)
        xi = zeros(nx)
        yi = zeros(ny)

        rho = zeros(ny, nx)
        rhop = zeros(ny, nx, 2)
        rhom = zeros(ny, nx, 2)

        velx = zeros(ny, nx)
        velxp = zeros(ny, nx, 2)
        velxm = zeros(ny, nx, 2)
        vely = zeros(ny, nx)
        velyp = zeros(ny, nx, 2)
        velym = zeros(ny, nx, 2)

        eps = zeros(ny, nx)
        epsp = zeros(ny, nx, 2)
        epsm = zeros(ny, nx, 2)

        press = zeros(ny, nx)
        pressp = zeros(ny, nx, 2)
        pressm = zeros(ny, nx, 2)

        q = zeros(ny, nx, 4)
        qp = zeros(ny, nx, 4, 2)
        qm = zeros(ny, nx, 4, 2)
        g = 3

        new(x,
            y,
            xi,
            yi,
            rho,
            rhop,
            rhom,
            velx,
            velxp,
            velxm,
            vely,
            velyp,
            velym,
            eps,
            epsp,
            epsm,
            press,
            pressp,
            pressm,
            q,
            qp,
            qm,
            nx,
            ny,
            g)
    end
end

#2-dim grid
function grid_setup(self::data2d, xmin, xmax, ymin, ymax)
    dx = (xmax - xmin) / (self.nx -self.g*2 - 1)
    dy = (ymax - ymin) / (self.ny -self.g*2 - 1)

    xmin = xmin - self.g*dx
    xmax = xmax + self.g*dx
    ymin = ymin - self.g*dy
    ymax = ymax + self.g*dy

    #x setup
    for i = 0:self.nx-1
        self.x[i+1] = xmin + i*dx
    end

    for i = 0:self.nx-1
        self.xi[i+1] = self.x[i+1] - 0.5dx
    end

    #y setup
    for i = 0:self.ny-1
        self.y[i+1] = ymin + i*dy
    end

    for i = 0:self.ny-1
        self.yi[i+1] = self.y[i+1] - 0.5dy
    end

    return self
end

#splice 2dim data into x-axis pencils
#(to save time we copy only the edge values)
function splice(hyd::data2d, j::Int64, dim::Int64)

    if dim == 1 #x-dim
        #println("x j=$j $(hyd.ny)")
        @assert j <= hyd.ny
        nn = hyd.nx
        xs = 1:nn
        ys = j
        hyd1 = data1d(nn)
        hyd1.x = hyd.x
        hyd1.xi = hyd.xi

        hyd1.velxp[:] = vec(hyd.velxp[ys, xs, dim])
        hyd1.velxm[:] = vec(hyd.velxm[ys, xs, dim])

        hyd1.velyp[:] = vec(hyd.velyp[ys, xs, dim])
        hyd1.velym[:] = vec(hyd.velym[ys, xs, dim])

        hyd1.qm[:, 2] = vec(hyd.qm[ys, xs, 2, dim])
        hyd1.qp[:, 2] = vec(hyd.qp[ys, xs, 2, dim])

        hyd1.qm[:, 3] = vec(hyd.qm[ys, xs, 3, dim])
        hyd1.qp[:, 3] = vec(hyd.qp[ys, xs, 3, dim])

    elseif dim == 2 #y-dim
        #println("y j=$j $(hyd.nx)")
        @assert j <= hyd.nx
        nn = hyd.ny
        xs = j
        ys = 1:nn
        hyd1 = data1d(nn)
        hyd1.x = hyd.y
        hyd1.xi = hyd.yi

        hyd1.velxp[:] = vec(hyd.velyp[ys, xs, dim])
        hyd1.velxm[:] = vec(hyd.velym[ys, xs, dim])

        hyd1.velyp[:] = vec(hyd.velxp[ys, xs, dim])
        hyd1.velym[:] = vec(hyd.velxm[ys, xs, dim])

        hyd1.qm[:, 2] = vec(hyd.qm[ys, xs, 3, dim])
        hyd1.qp[:, 2] = vec(hyd.qp[ys, xs, 3, dim])

        hyd1.qm[:, 3] = vec(hyd.qm[ys, xs, 2, dim])
        hyd1.qp[:, 3] = vec(hyd.qp[ys, xs, 2, dim])

    else
        error("dim = $dim is not valid")
    end

    hyd1.g = hyd.g

    hyd1.rhom[:] = vec(hyd.rhom[ys, xs, dim])
    hyd1.rhop[:] = vec(hyd.rhop[ys, xs, dim])

    hyd1.epsm[:] = vec(hyd.epsm[ys, xs, dim])
    hyd1.epsp[:] = vec(hyd.epsp[ys, xs, dim])

    hyd1.pressm[:] = vec(hyd.pressm[ys, xs, dim])
    hyd1.pressp[:] = vec(hyd.pressp[ys, xs, dim])

    hyd1.qm[:, 1] = vec(hyd.qm[ys, xs, 1, dim])
    hyd1.qp[:, 1] = vec(hyd.qp[ys, xs, 1, dim])

    hyd1.qm[:, 4] = vec(hyd.qm[ys, xs, 4, dim])
    hyd1.qp[:, 4] = vec(hyd.qp[ys, xs, 4, dim])

    return hyd1
end

max_stars_repo_name: natj/hydro

max_issues_repo_name: natj/hydro

max_forks_repo_name: natj/hydro


text:
using Documenter, Modia, Modia.ModiaResult, ModiaPlot_PyPlot

makedocs(
  #modules  = [Modia],
  sitename = "Modia",
  authors  = "Hilding Elmqvist (Mogram) and Martin Otter (DLR-SR)",
  format = Documenter.HTML(prettyurls = false),
  pages    = [
     "Home"      => "index.md",
	 "Tutorial"  => [
       "tutorial/Tutorial.md"
       "tutorial/GettingStarted.md"
       "tutorial/Modeling.md"
       "tutorial/Simulation.md"
       "tutorial/FloatingPointTypes.md"
       "tutorial/Appendix.md"
      ],
  	 "Functions" => "Functions.md",
     "Internal"  => "Internal.md"
  ]
)


text:
"""
    InvertedList{I<:Unsigned, U<:Unsigned}

Basic structure which corresponds to the points found within
a Voronoi cell. The fields `idxs` contains the indices of the
points while `codes` contains quantized vector data.
"""
struct InvertedList{I<:Unsigned, U<:Unsigned}
    idxs::Vector{I}
    codes::Vector{Vector{U}}
end

Base.show(io::IO, ivlist::InvertedList{I,U}) where {I,U} = begin
    n = length(ivlist.idxs)
    print(io, "InvertedList{$I,$U}, $n vectors")
end

"""
Simple alias for `Vector{InvertedList{I<:Unsigned, U<:Unsigned}}`.
"""
const InvertedIndex{I,U} = Vector{InvertedList{I,U}}

"""
    IVFADCIndex{U<:Unsigned, I<:Unsigned, Dc<:Distances.PreMetric, Dr<:Distances.PreMetric, T<:AbstractFloat, Q<:AbstractCoarseQuantizer{Dc,T}}

The inverse file system object. It allows for approximate nearest
neighbor search into the contained vectors.

# Fields
  * `coarse_quantizer::AbstractCoarseQuantizer{Dc,T}` contains the coarse vectors
  * `residual_quantizer::QuantizedArrays.OrthogonalQuantizer{U,Dr,T,2}`
is employed to quantize vectors when adding to the index
  * `inverse_index::InvertedIndex{I,U}` is the actual inverse index employed
to perform the search.
"""
struct IVFADCIndex{U<:Unsigned,
                   I<:Unsigned,
                   Dc<:Distances.PreMetric,
                   Dr<:Distances.PreMetric,
                   T<:AbstractFloat,
                   Q<:AbstractCoarseQuantizer{Dc,T}}
    coarse_quantizer::Q
    residual_quantizer::QuantizedArrays.OrthogonalQuantizer{U,Dr,T,2}
    inverse_index::InvertedIndex{I,U}
end

"""
    length(ivfadc::IVFADCIndex)

Returns the number of vectors indexed by `ivfadc`.
"""
Base.length(ivfadc::IVFADCIndex) =
    mapreduce(ivlist->length(ivlist.idxs), +, ivfadc.inverse_index)

"""
    size(ivfadc::IVFADCIndex)

Returns a tuple with the dimensionality and number of the vectors indexed by `ivfadc`.
"""
Base.size(ivfadc::IVFADCIndex) = (size(ivfadc.coarse_quantizer, 1), length(ivfadc))
Base.size(ivfadc::IVFADCIndex, i::Int) = size(ivfadc)[i]

Base.show(io::IO, ivfadc::IVFADCIndex{U,I,Dc,Dr,T,Q}) where {U,I,Dc,Dr,T,Q} = begin
    nvars, nvectors = size(ivfadc)
    idxsize = sizeof(I)
    m = length(ivfadc.residual_quantizer.codebooks)
    compsize = sizeof(U)
    codesize = m * compsize + idxsize
    cqstr = ifelse(Q <: HNSWQuantizer, "HNSW", "naive")
    print(io, "IVFADCIndex, $cqstr coarse quantizer, $codesize-byte encoding ($idxsize + $compsize×$m), $nvectors $T vectors")
end

"""
    IVFADCIndex(data [;kwargs])

Main constructor for building an inverse file system for billion-scale ANN search.

# Arguments
  * `Matrix{T<:AbstractFloat}` input data

# Keyword arguments
  * `kc::Int=DEFAULT_COARSE_K` number of clusters (Voronoi cells) to employ
in the coarse quantization step
  * `k::Int=DEFAULT_QUANTIZATION_K` number of residual quantization levels to use
  * `m::Int=DEFAULT_QUANTIZATION_M` number of residual quantizers to use
  * `coarse_quantizer::Symbol=DEFAULT_COARSE_QUANTIZER` coarse quantizer
  * `coarse_distance=DEFAULT_COARSE_DISTANCE` coarse quantization distance
  * `quantization_distance=DEFAULT_QUANTIZATION_DISTANCE` residual quantization distance
  * `quantization_method=DEFAULT_QUANTIZATION_METHOD` residual quantization method
  * `coarse_maxiter=DEFAULT_COARSE_MAXITER` number of clustering iterations for obtaining
the coarse vectors
  * `quantization_maxiter=DEFAULT_QUANTIZATION_MAXITER` number of clustering iterations for
residual quantization
  * `index_type=UInt32` type for the indexes of the vectors in the inverted list
"""
function IVFADCIndex(data::Matrix{T};
                     kc::Int=DEFAULT_COARSE_K,
                     k::Int=DEFAULT_QUANTIZATION_K,
                     m::Int=DEFAULT_QUANTIZATION_M,
                     coarse_quantizer::Symbol=DEFAULT_COARSE_QUANTIZER,
                     coarse_distance::Distances.PreMetric=DEFAULT_COARSE_DISTANCE,
                     quantization_distance::Distances.PreMetric=DEFAULT_QUANTIZATION_DISTANCE,
                     quantization_method::Symbol=DEFAULT_QUANTIZATION_METHOD,
                     coarse_maxiter::Int=DEFAULT_COARSE_MAXITER,
                     quantization_maxiter::Int=DEFAULT_QUANTIZATION_MAXITER,
                     index_type::Type{I}=UInt32
                    ) where {I<:Unsigned, T<:AbstractFloat}
    # Checks
    nrows, nvectors = size(data)
    bits_required = ceil(Int, log2(nvectors))
    @assert kc >= 2 "Number of coarse clusters has to be >= 2"
    @assert k <= nvectors "Number of quantization levels  has to be <= $nvectors"
    @assert m >= 1 && m <= nrows "Number of codebooks has to be between 1 and $nrows"
    @assert coarse_quantizer in [:naive, :hnsw] "Coarse quantizer can be :naive or :hnsw only"
    @assert coarse_maxiter > 0 "Number of clustering iterations has to be > 0"
    @assert quantization_maxiter > 0 "Number of clustering iterations has to be > 0"
    @assert QuantizedArrays.TYPE_TO_BITS[index_type] >=
        bits_required "$nvectors vectors require at least $bits_required index bits"

    # Run kmeans, build coarse quantizer
    @debug "Clustering..."
    cmodel = kmeans(data,
                    kc,
                    maxiter=coarse_maxiter,
                    distance=coarse_distance,
                    init=:kmpp,
                    display=:none)

    # Calculate residuals
    @debug "Residual calculation..."
    residuals = _build_residuals(cmodel, data)

    # Build residuals quantizer
    @debug "Building quantizer..."
    rq = build_quantizer(residuals,
                         k=k,
                         m=m,
                         method=quantization_method,
                         distance=quantization_distance,
                         maxiter=quantization_maxiter)

    # Quantize residuals for each cluster
    @debug "Building inverted index..."
    ii = _build_inverted_index(rq, cmodel, residuals, index_type=index_type)

    # Return
    @debug "Building coarse quantizer..."
    if coarse_quantizer == :naive
        cq = NaiveQuantizer{typeof(coarse_distance),
                eltype(cmodel.centers)}(cmodel.centers)
    else
        hnsw = HierarchicalNSW([cmodel.centers[:,i] for i in 1:kc],
                                metric=coarse_distance)
        add_to_graph!(hnsw)
        cq = HNSWQuantizer(hnsw)
    end
    return IVFADCIndex(cq, rq, ii)
end

function _build_residuals(km::KmeansResult, data::Matrix{T}) where {T}
    residuals = similar(data)
    for cluster in 1:nclusters(km)
        bidxs = km.assignments .== cluster
        @. residuals[:, bidxs] = data[:, bidxs] - km.centers[:, cluster]
    end
    return residuals
end

function _build_inverted_index(rq::QuantizedArrays.OrthogonalQuantizer{U,D,T,2},
                               km::KmeansResult,
                               data::Matrix{T};
                               index_type::Type{I}=UInt32
                              ) where {U,I<:Unsigned,D,T}
    n = nclusters(km)
    invindex = InvertedIndex{I,U}(undef, n)
    for cluster in 1:n
        idxs = findall(isequal(cluster), km.assignments)
        qdata = QuantizedArrays.quantize_data(rq, data[:, idxs])
        ivlist = InvertedList{I,U}(
                    idxs .- one(I),
                    [qdata[:, j] for j in 1:length(idxs)])
        invindex[cluster] = ivlist
    end
    return invindex
end

"""
    knn_search(ivfadc, point, k[; w=1])

Searches at most `k` closest neighbors of `point` in the index `ivfadc`;
the neighbors will be searched for in the points contained in the closest
`w` clusters.
"""
function knn_search(ivfadc::IVFADCIndex{U,I,Dc,Dr,T},
                    point::Vector{T},
                    k::Int;
                    w::Int=1
                   ) where {U,I,Dc,Dr,T}
    # Checks and initializations
    @assert k >= 1 "Number of neighbors must be k >= 1"
    @assert w >= 1 "Number of clusters to search in must be w >= 1"
    cq = ivfadc.coarse_quantizer
    nclusters = size(cq, 2)
    rq_cbooks = ivfadc.residual_quantizer.codebooks
    m, n = length(rq_cbooks), length(point)
    w = min(w, nclusters)

    # Find the 'w' closest coarse vectors and calculate residuals
    closest_clusters, coarse_distances = coarse_search(cq, point, w)
    residuals = _closest_cluster_residuals(cq, point, closest_clusters)

    # Calculate distances between point and the vectors
    # from the clusters using the residual distances.
    distances = Vector{T}()
    neighbors = SortedMultiDict{T,I}()
    maxdist = zero(T)
    difftables = Vector{LittleDict{U,T}}(undef, m)
    @inbounds for j = eachindex(closest_clusters)
        dc = coarse_distances[j]
        # Calculate all residual distances
        # (between the vector and the codebooks of the residual quantizer)
        for i = eachindex(difftables)
            rr = QuantizedArrays.rowrange(n, m, i)
            diffs = colwise(Dc(), rq_cbooks[i].vectors, residuals[rr, j])
            difftables[i] = LittleDict{U,T}(rq_cbooks[i].codes, diffs)
        end
        # Loop through the iverted list and calculate
        # the actual (quantized) distances through lookup;
        # use SortedMultiDict as a maxheap for `k` neighbors
        ivlist = ivfadc.inverse_index[closest_clusters[j]]
        for i = eachindex(ivlist.idxs)
            d = dc
            @simd for ii = eachindex(ivlist.codes[i])
                # add difference vector for the ii'th code of the i'th subquantizer
                d += difftables[ii][ivlist.codes[i][ii]]
            end
            if length(neighbors) < k
                push!(neighbors, d => ivlist.idxs[i])
                maxdist, _ = last(neighbors)
            elseif maxdist > d
                delete!((neighbors, lastindex(neighbors)))  # delete the max key
                push!(neighbors, d => ivlist.idxs[i])
                maxdist, _ = last(neighbors)
            end
        end
    end
    return collect(values(neighbors)), collect(keys(neighbors))
end

function knn_search(ivfadc::IVFADCIndex{U,I,Dc,Dr,T},
                    points::Vector{Vector{T}},
                    k::Int;
                    w::Int=1
                   ) where {U,I,Dc,Dr,T}
    n = length(points)
    idxs = Vector{Vector{I}}(undef, n)
    dists = Vector{Vector{T}}(undef, n)
    for i = eachindex(points)  # Threads.@threads for i = ... for multi-threading
        idxs[i], dists[i] = knn_search(ivfadc, points[i], k, w=w)
    end
    return idxs, dists
end


text:
using ImplicitGlobalGrid

@views d_xa(A) = A[2:end  , :     , :     ] .- A[1:end-1, :     , :     ];
@views d_xi(A) = A[2:end  ,2:end-1,2:end-1] .- A[1:end-1,2:end-1,2:end-1];
@views d_ya(A) = A[ :     ,2:end  , :     ] .- A[ :     ,1:end-1, :     ];
@views d_yi(A) = A[2:end-1,2:end  ,2:end-1] .- A[2:end-1,1:end-1,2:end-1];
@views d_za(A) = A[ :     , :     ,2:end  ] .- A[ :     , :     ,1:end-1];
@views d_zi(A) = A[2:end-1,2:end-1,2:end  ] .- A[2:end-1,2:end-1,1:end-1];
@views  inn(A) = A[2:end-1,2:end-1,2:end-1]

@views function diffusion3D()
    # Physics
    lam        = 1.0;                                       # Thermal conductivity
    cp_min     = 1.0;                                       # Minimal heat capacity
    lx, ly, lz = 10.0, 10.0, 10.0;                          # Length of computational domain in dimension x, y and z

    # Numerics
    nx, ny, nz = 128, 128, 128;                             # Number of gridpoints in dimensions x, y and z
    nt         = 10000;                                     # Number of time steps
    me, dims   = init_global_grid(nx, ny, nz);              # Initialize the implicit global grid
    dx         = lx/(nx_g()-1);                             # Space step in dimension x
    dy         = ly/(ny_g()-1);                             # ...        in dimension y
    dz         = lz/(nz_g()-1);                             # ...        in dimension z

    # Array initializations
    T     = zeros(nx,   ny,   nz  );
    Cp    = zeros(nx,   ny,   nz  );
    dTedt = zeros(nx-2, ny-2, nz-2);
    qx    = zeros(nx-1, ny-2, nz-2);
    qy    = zeros(nx-2, ny-1, nz-2);
    qz    = zeros(nx-2, ny-2, nz-1);

    # Initial conditions (heat capacity and temperature with two Gaussian anomalies each)
    Cp .= cp_min .+ [5*exp(-((x_g(ix,dx,Cp)-lx/1.5))^2-((y_g(iy,dy,Cp)-ly/2))^2-((z_g(iz,dz,Cp)-lz/1.5))^2) +
                     5*exp(-((x_g(ix,dx,Cp)-lx/3.0))^2-((y_g(iy,dy,Cp)-ly/2))^2-((z_g(iz,dz,Cp)-lz/1.5))^2) for ix=1:size(T,1), iy=1:size(T,2), iz=1:size(T,3)]
    T  .= [100*exp(-((x_g(ix,dx,T)-lx/2)/2)^2-((y_g(iy,dy,T)-ly/2)/2)^2-((z_g(iz,dz,T)-lz/3.0)/2)^2) +
            50*exp(-((x_g(ix,dx,T)-lx/2)/2)^2-((y_g(iy,dy,T)-ly/2)/2)^2-((z_g(iz,dz,T)-lz/1.5)/2)^2) for ix=1:size(T,1), iy=1:size(T,2), iz=1:size(T,3)]

    # Time loop
    dt = min(dx*dx,dy*dy,dz*dz)*cp_min/lam/8.1;                                               # Time step for the 3D Heat diffusion
    for it = 1:nt
        qx    .= -lam.*d_xi(T)./dx;                                                           # Fourier's law of heat conduction: q_x   = -λ δT/δx
        qy    .= -lam.*d_yi(T)./dy;                                                           # ...                               q_y   = -λ δT/δy
        qz    .= -lam.*d_zi(T)./dz;                                                           # ...                               q_z   = -λ δT/δz
        dTedt .= 1.0./inn(Cp).*(-d_xa(qx)./dx .- d_ya(qy)./dy .- d_za(qz)./dz);               # Conservation of energy:           δT/δt = 1/cₚ (-δq_x/δx - δq_y/dy - δq_z/dz)
        T[2:end-1,2:end-1,2:end-1] .= inn(T) .+ dt.*dTedt;                                    # Update of temperature             T_new = T_old + δT/δt
        update_halo!(T);                                                                      # Update the halo of T
    end

    finalize_global_grid();                                                                   # Finalize the implicit global grid
end

diffusion3D()


text:
# SOCP Lasso Testproblem

using COSMO, Test, LinearAlgebra, SparseArrays, Random

# generate problem data
rng = Random.MersenneTwister(12345)
n = 8
m = 50 * n
F = rand(rng, m, n)

vtrue = sprand(rng, n, 1, 0.1)
noise = 0.1 * rand(rng, m, 1)
b = F * vtrue + noise
μMax = norm(F' * b, Inf)
μ = 0.1 * μMax

# define lasso problem as SOCP

A1 = -sparse([1 zeros(1, 2 * n + 1) 1 zeros(1, m);
  -1 zeros(1, 2 * n) 1 zeros(1, m + 1);
  zeros(m, 1) -2 * F zeros(m, n + 2) Matrix(1.0I, m, m)])

A2 = -sparse([zeros(n, 1) Matrix(1.0I, n, n) -Matrix(1.0I, n, n) zeros(n, m + 2);
  zeros(n, 1) -Matrix(1.0I, n, n) -Matrix(1.0I, n, n) zeros(n, m + 2)])
A3 = -sparse([zeros(1, 2 * n + 1) -1 zeros(1, m + 1);
 zeros(1, 2 * n + 2) -1 zeros(1, m);
 zeros(m, 2 * n + 3) -Matrix(1.0I, m, m)])
b1 = [1; 1; -2 * b]
b2 = zeros(2 * n)
b3 = zeros(m + 2)

q = vec([1; zeros(n); μ * ones(n, 1); zeros(m + 2, 1)])
P = spzeros(length(q), length(q))

cs1 = COSMO.Constraint(A1 ,b1, COSMO.ZeroSet)
cs2 = COSMO.Constraint(A2 ,b2, COSMO.Nonnegatives)
cs3 = COSMO.Constraint(A3 ,b3, COSMO.SecondOrderCone)

# Solve with OSSDP
model = COSMO.Model()
assemble!(model, P, q, [cs1; cs2; cs3])
res = COSMO.optimize!(model);

@testset "SOCP - Lasso" begin
@test res.status == :Solved
@test isapprox(res.obj_val,0.4422849814458825,atol=1e-2)
end
nothing


text:
"""
The base type for all algorithms representing various distances metrics.
"""
abstract type AbstractDistance end

"""
    struct SnodeDistanceDijkstra(source_node::Int, target_node::Int, s::Int) <: AbstractDistance

Represent a distance between two nodes of the hypergraph `h`, which is
the minimum `s`-walk lenght between the two nodes. An `s`-walk between nodes
is a sequence of nodes that pairwise share at least `s` edges.
"""
struct SnodeDistanceDijkstra <: AbstractDistance
    source_node::Int
    target_node::Int
    s::Int
end

"""
    struct SedgeDistanceDijkstra(source_edge::Int, target_edge::Int, s::Int) <: AbstractDistance

Represent a distance between two hyperedges of the hypergraph `h`, which is
the minimum `s`-walk lenght between the two hyperedge. An `s`-walk between edges is a sequence
of edges such that consecutive pairwise edges intersect in at least `s` nodes.
"""
struct SedgeDistanceDijkstra <: AbstractDistance
    source_edge::Int
    target_edge::Int
    s::Int
end

"""
    distance(h::Hypergraph, distance_method::SnodeDistanceDijkstra)

Return the shortest `distance_method.s`-walk distance between the `distance_method.source_node` and
the node `distance_method.target_node` in the hypergraph `h`.

NOTE
The concepts of `s`-distance and `s`-walk have been defined in the
Python library [HyperNetX](https://github.com/pnnl/HyperNetX)

From [HyperNetX](https://pnnl.github.io/HyperNetX/build/_modules/classes/hypergraph.html#Hypergraph.distance)
The `s`-distance is the shortest `s`-walk length between two nodes.
An `s`-walk between nodes is a sequence of nodes that pairwise share
at least `s` edges. The length of the shortest `s`-walk is 1 less than
the number of nodes in the path sequence. If no such path exists returns typemax(T).
"""
function distance(h::Hypergraph, distance_method::SnodeDistanceDijkstra)
    checkbounds(h.v2he, distance_method.source_node)
    checkbounds(h.v2he, distance_method.target_node)
    A = adjacency_matrix(h; s=distance_method.s)
    g = LightGraphs.Graph(A)
    dj = LightGraphs.dijkstra_shortest_paths(g, distance_method.source_node)
    dj.dists[distance_method.target_node]
end

"""
    distance(h::Hypergraph, distance_method::SedgeDistanceDijkstra)

Return the shortest `distance_method.s`-walk distance between the `distance_method.source_edge` and
the node `distance_method.target_edge` in the hypergraph `h`.

NOTE
The concepts of `s`-distance and `s`-walk have been defined in the
Python library [HyperNetX](https://github.com/pnnl/HyperNetX)

From [HyperNetX](https://pnnl.github.io/HyperNetX/build/_modules/classes/hypergraph.html#Hypergraph.edge_distance)
The `s`-distance is the shortest `s`-walk length between the edges.
An `s`-walk between edges is a sequence of edges such that consecutive pairwise
edges intersect in at least `s` nodes. The length of the shortest `s`-walk is 1 less than
the number of edges in the path sequence. If no such path exists returns typemax(T).
"""
function distance(h::Hypergraph, distance_method::SedgeDistanceDijkstra)
    checkbounds(h.he2v, distance_method.source_edge)
    checkbounds(h.he2v, distance_method.target_edge)
    A = edge_adjacency_matrix(h; s=distance_method.s)
    g = LightGraphs.Graph(A)
    dj = LightGraphs.dijkstra_shortest_paths(g, distance_method.source_edge)
    dj.dists[distance_method.target_edge]
end


text:
push!(LOAD_PATH,"../src/")
using Firebase
using Documenter

DocMeta.setdocmeta!(Firebase, :DocTestSetup, :(using Firebase); recursive=true)

makedocs(;
    modules=[Firebase],
    authors="Ashwani Rathee",
    repo="github.com/ashwani-rathee/Firebase.jl/blob/{commit}{path}#{line}",
    sitename="Firebase.jl",
    format=Documenter.HTML(;
        prettyurls=Base.get(ENV, "CI", "false") == "true",
        canonical="https://ashwani-rathee.github.io/Firebase.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Cloud Firestore" => "firestore.md",
        "Cloud Authentication" => "authentication.md",
        "Realtime Database" => "realtime.md",
        "Cloud Storage" => "storage.md",
        "Utility Functions" => "utility.md",
    ],
)

deploydocs(;
    repo="github.com/ashwani-rathee/Firebase.jl",
    devbranch="main",
    push_preview = true
)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 24734878-7d42-11eb-10e7-35421889faaf
using GR

# ╔═╡ 167d66f2-7d42-11eb-1967-3fff5fa301dd
histogram(randn(10000))

# ╔═╡ Cell order:
# ╠═24734878-7d42-11eb-10e7-35421889faaf
# ╠═167d66f2-7d42-11eb-1967-3fff5fa301dd


text:
import ViscousFlow: Fields, SaddlePointSystems

using LinearAlgebra

@testset "Saddle-Point Systems" begin

nx = 130; ny = 130
Lx = 2.0
dx = Lx/(nx-2)
w = Nodes(Dual,(nx,ny))

L = plan_laplacian(size(w),with_inverse=true)
Linv(x) = L\x
#L⁻¹(w::T) where {T} = L\w

n = 128
θ = range(0,stop=2π,length=n+1)
R = 0.5
xb = 1.0 .+ R*cos.(θ[1:n])
yb = 1.0 .+ R*sin.(θ[1:n])
ds = (2π/n)*R
X = VectorData(xb,yb)
f = ScalarData(X)

E = Regularize(X,dx;ddftype=Fields.Roma,issymmetric=true)
Hmat,Emat = RegularizationMatrix(E,f,w)

PS = SaddleSystem((w,f),(Linv,Hmat,Emat),issymmetric=true,isposdef=true)

ψb = ScalarData(X)
w = Nodes(Dual,(nx,ny))
ψb .= -(xb .- 1)
f .= ones(Float64,n)*ds
ψ = Nodes(Dual,w)
ψ,f = PS\(w,ψb)

fex = -2*cos.(θ[1:n])
@test norm(f./ds-fex,Inf) ≈ 0.035288024
@test ψ[nx,65] ≈ -ψ[1,65]
@test ψ[65,ny] ≈ ψ[65,1]

ru = ones(Float64,2)
rf = Vector{Float64}()
A⁻¹(u::Vector{Float64}) = u
B₁ᵀ(f::Vector{Float64}) = zeros(Float64,2)
B₂(u::Vector{Float64}) = Vector{Float64}()
sys = SaddleSystem((ru,rf),(A⁻¹,B₁ᵀ,B₂))

u,fnull = sys\(ru,rf)
@test u ≈ [1.0,1.0]
@test length(fnull) == 0

sysys = (PS,sys)
rhs = ((w,ψb),(ru,rf))
(ψ,f),(u,fnull) = sysys\rhs

@test ψ[nx,65] ≈ -ψ[1,65]
@test ψ[65,ny] ≈ ψ[65,1]
@test u ≈ [1.0,1.0]
@test length(fnull) == 0

end


text:
for (A,F) in ((:gamma, :arb_gamma),
              (:lgamma, :arb_lgamma),(:rgamma, :arb_rgamma), (:digamma, :arb_digamma),
              (:zeta, :arb_zeta),

              (:erf, :arb_hypgeom_erf), (:erfc, :arb_hypgeom_erfc), (:erfi, :arb_hypgeom_erfi),
              (:ei, :arb_hypgeom_ei), (:si, :arb_hypgeom_si), (:ci, :arb_hypgeom_ci),
              (:shi, :arb_hypgeom_shi), (:chi, :arb_hypgeom_chi),
             )
    @eval begin
        function ($A)(x::ArbReal{P}, prec::Int=P) where {P}
            z = ArbReal{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, prec)
            return z
         end
    end
end

const Cint0 = zero(Cint)
const FmpzZero = Fmpz(0)

function lambertw(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_lambertw), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Cint, Clong), z, x, Cint0, P)
    return z
end

lambertw(x::ArbFloat{P}) where {P} = ArbFloat{P}(lambertw(ArbReal{P}(x)))

function lambertw(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_lambertw), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{Fmpz}, Cint, Clong), z, x, FmpzZero, 0, 128)
    return z
end

for (A,F) in ((:gamma, :acb_gamma),
              (:lgamma, :acb_lgamma),(:rgamma, :acb_rgamma), (:digamma, :acb_digamma),

              (:logsinpi, :acb_log_sin_pi), (:barnesg, :acb_barnes_g), (:logbarnesg, :acb_log_barnes_g),
              (:agm1, :acb_agm1),

              (:erf, :acb_hypgeom_erf), (:erfc, :acb_hypgeom_erfc), (:erfi, :acb_hypgeom_erfi),
              (:ei, :acb_hypgeom_ei), (:si, :acb_hypgeom_si), (:ci, :acb_hypgeom_ci),
              (:shi, :acb_hypgeom_shi), (:chi, :acb_hypgeom_chi),
             )
    @eval begin
        function ($A)(x::ArbComplex{P}, prec::Int=P) where P
            z = ArbComplex{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, prec)
            return z
         end
    end
end

for (A,F) in (
              (:eta, :acb_dirichlet_eta), (:xi, :acb_dirichlet_xi)
             )
    @eval begin
        function ($A)(x::ArbReal{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            xc = ArbComplex{P}(x)
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, xc, prec)
            rea = real(z)
            return rea
         end
    end
end
#=
for (A,F) in ((:elliptick, :elliptic_k), (:elliptice, :elliptic_e))
    @eval begin
        function ($A)(k::ArbReal{P}, prec::Int=P) where {P}
           m = square(k)
           return ($F)(m, prec)
        end
    end
end
=#
for (A,F) in ((:agm, :arb_agm), )
    @eval begin
        function ($A)(x::ArbReal{P}, y::ArbReal{P}, prec::Int=P) where {P}
            z = ArbReal{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, y, prec)
            return z
         end
    end
end

function polylog(s::ArbComplex{P}, z::ArbComplex{P}, prec::Int=P) where {P}
    w = ArbComplex{P}()
    ccall(@libarb(acb_polylog), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Cint), w, s, z, P)
    return w
end

function polylog(s::ArbReal{P}, z::ArbReal{P}, prec::Int=P) where {P}
    sc = ArbComplex(s)
    zc = ArbComplex(z)
    wc = polylog(sc, zc, prec)
    return wc
end

function polylog(s::ArbFloat{P}, z::ArbFloat{P}, prec::Int=P) where {P}
    sc = ArbComplex(s)
    zc = ArbComplex(z)
    wc = polylog(sc, zc, prec)
    return wc
end

function polylog(s::Int, z::ArbComplex{P}, prec::Int=P) where {P}
    w = ArbComplex{P}()
    ccall(@libarb(acb_polylog_si), Cvoid, (Ref{ArbComplex}, Cint, Ref{ArbComplex}, Cint), w, s, z, P)
    return w
end

function polylog(s::Int, z::ArbReal{P}, prec::Int=P) where {P}
    zc = ArbComplex(z)
    wc = polylog(s, zc, prec)
    return wc
end

function polylog(s::Int, z::ArbFloat{P}, prec::Int=P) where {P}
    zc = ArbComplex(z)
    wc = polylog(s, zc, prec)
    return wc
end

#=
for (A,F) in ((:ellipticp, :acb_elliptic_p), (:ellipticpi, :acb_elliptic_pi),
              (:ellipticzeta, :acb_elliptic_zeta), (:ellipticsigma, :acb_elliptic_sigma),
             )
    @eval begin
        function ($A)(x::ArbComplex{P}, y::ArbComplex{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, y, prec)
            return z
         end
    end
end

for (A,F) in ((:ellipticpi, :acb_elliptic_pi), (:ellipticp, :acb_elliptic_p),
              (:ellipticzeta, :acb_elliptic_zeta), (:ellipticsigma, :acb_elliptic_sigma),
             )
    @eval begin
        function ($A)(x::ArbReal{P}, y::ArbReal{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            xc = ArbComplex{P}(x)
            yc = ArbComplex{P}(y)
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, xc, yc, prec)
            rea = real(z)
            return rea
         end
    end
end
=#
#=
function ellipticpi(n::ArbComplex{P}, k::ArbComplex{P}, prec::Int=P) where {P}
    m = square(k)
    return elliptic_pi(n, m, prec)
end
=#

for (A,F) in ((:gamma, :arb_gamma),
              (:lgamma, :arb_lgamma),(:rgamma, :arb_rgamma), (:digamma, :arb_digamma),
              (:zeta, :arb_zeta),

              (:erf, :arb_hypgeom_erf), (:erfc, :arb_hypgeom_erfc), (:erfi, :arb_hypgeom_erfi),
              (:ei, :arb_hypgeom_ei), (:si, :arb_hypgeom_si), (:ci, :arb_hypgeom_ci),
              (:shi, :arb_hypgeom_shi), (:chi, :arb_hypgeom_chi),
             )
    @eval begin
        function ($A)(x::ArbFloat{P}, prec::Int=P) where {P}
            z = ArbReal{P}()
            xb = ArbReal{P}(x)
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, xb, prec)
            return midpoint_byref(z)
         end
    end
end

for (A,F) in ((:eta, :acb_dirichlet_eta), (:xi, :acb_dirichlet_xi))
    @eval begin
        function ($A)(x::ArbFloat{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            xc = ArbComplex{P}(x)
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, xc, prec)
            return midpoint_byref(real(z))
         end
    end
end

function agm1(x::ArbReal{P}) where {P}
    x1 = ArbComplex{P}(x)
    w  = agm1(x1)
    z = real(w)
    return z
end

function agm1(x::ArbFloat{P}) where {P}
    x1 = ArbReal{P}(x)
    w  = agm1(x1)
    z = midpoint_byref(w)
    return z
end

function agm(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    w  = agm(x1, y1)
    z = midpoint_byref(w)
    return z
end

#=
function erfcx(z::ArbComplex{P}) where {P}
    hiprec = round(Int, P*1.25)
    setprecision(ArbComplex, hiprec)
    w = ArbComplex(real(z), imag(z), bits=hiprec)
    ww = w*w
    a  = exp(ww)
    b  = erfc(w)
    res  = a * b
    setprecision(ArbComplex, P)
    return ArbComplex(real(res), imag(res), bits=P-extrabits())
end

function erfcx(x::ArbReal{P}) where {P}
    hiprec = round(Int, P*1.25)
    setprecision(ArbReal, hiprec)
    w = ArbReal(x, bits=hiprec)
    ww = w*w
    a  = exp(ww)
    b  = erfc(w)
    res  = a * b
    setprecision(ArbReal, P)
    return ArbReal(res, bits=P-extrabits())
end

function erfcx(x::ArbFloat{P}) where {P}
    res = erfcx(ArbReal(x))
    return ArbFloat(res)
end
=#

function erfcx(z::ArbComplex{P}; magnify=2.125) where {P}
    external_precision = P
    external_extrabits = extrabits()
    internal_precision = round(Int, P*magnify)
    internal_extrabits = 0
    setextrabits(internal_extrabits)
    setworkingprecision(ArbReal, internal_precision)

    re = round(real(z), sigdigits=internal_precision, base=2)
    im = round(imag(z), sigdigits=internal_precision, base=2)
    w  = ArbComplex(re, im)
    ww = w*w
    a  = exp(ww)
    b  = erfc(w)
    ab = a * b
    abre = round(real(ab), sigdigits=external_precision, base=2)
    abim = round(imag(ab), sigdigits=external_precision, base=2)
    re = ArbReal{P}(abre)
    im = ArbReal{P}(abim)
    z  = ArbComplex{P}(re, im)
    setextrabits(external_extrabits)
    setworkingprecision(ArbReal, external_precision)
    return z
end

function erfcx(x::ArbReal{P}; magnify=2.125) where {P}
    external_precision = P
    external_extrabits = extrabits()
    internal_precision = round(Int, P*magnify)
    internal_extrabits = 0
    setextrabits(internal_extrabits)
    setworkingprecision(ArbReal, internal_precision)

    w = round(x, sigdigits=internal_precision, base=2)
    ww = w*w
    a  = exp(ww)
    b  = erfc(w)
    ab = a * b
    abre = round(real(ab), sigdigits=external_precision, base=2)
    z = ArbReal{P}(abre)
    setextrabits(external_extrabits)
    setworkingprecision(ArbReal, external_precision)
    return z
end

function erfcx(x::ArbFloat{P}; magnify=2.125) where {P}
    res = erfcx(ArbReal(x); magnify=magnify)
    return ArbFloat(res)
end


text:
include("file_blacklist.jl")

# Formats the file names to title case
function apply_formatting(name)
    title = name
    if !(lowercase(name) in format_blacklist) && !(lowercase(name) in uppercase_list)
        title = titlecase(replace(name, "-" => " "))
    elseif lowercase(name) in uppercase_list
        title = uppercase(title)
    end

    return title
end

# Checks wether a path leads to a folder and a non blacklisted file
function folder_and_not_blacklisted(file)

    is_folder = isdir(file)
    black_listed = folder_blacklisted(file)

    return is_folder && !black_listed
end

# Checks wether a folder is blacklisted.
function folder_blacklisted(file)
    blacklisted = false

    for exclude in folder_black_list
        if contains(file, exclude)
            blacklisted = true
            break
        end
    end

    return blacklisted
end

# Checks wether a file is blacklisted.
function file_blacklisted(file)
    blacklisted = false

    for exclude in file_black_list

        if contains(file, exclude)
            blacklisted = true
            break
        end
    end

    return blacklisted
end


text:
export Descent, InvDecay, backend, set_backend!

using Flux
using Knet

import Flux.Optimise: apply!, Descent, InvDecay, Dense
import Flux: params!
import Base: show

function apply!(o::Descent, x, δ::Number)
    o.eta * δ
end

function apply!(o::InvDecay, x, δ::Number)
    γ = o.gamma
    n = get!(o.state, x, 0)
    o.state[x] = n + 1
    δ / (1 + γ * n)
end

#####
# extend layers defined in Flux to support Knet
#####

# relief the AbstractArray constraint to support KnetArray
function (a::Dense)(x)
  W, b, σ = a.W, a.b, a.σ
  σ.(W*x .+ b)
end

function Dense(in::Integer, out::Integer, σ = identity; initW = Flux.glorot_uniform, initb = n -> zeros(Float32, n), backend=:Zygote)
    if backend == :Zygote
        Dense(initW(out, in), initb(out), σ)
    elseif backend == :Knet
        Dense(Knet.param(initW(out, in)), Knet.param(initb(out)), σ)
    else
        throw(ArgumentError("unknown backend $backend"))
    end
end

params!(p::Flux.Params, x::Knet.Param, seen = IdSet()) = push!(p, x)

const BACKENDS_CACHE = IdDict()
const SUPPORTED_BACKENDS = Set([:Zygote, :Knet, :Any])

function set_backend!(m, x)
    if x in SUPPORTED_BACKENDS
        BACKENDS_CACHE[m] = x
    else
        @error "unknown backend $x, supported backends are $SUPPORTED_BACKENDS"
    end
end

backend(m) = haskey(BACKENDS_CACHE, m) ? BACKENDS_CACHE[m] : :Any
backend(::Dense) = :Zygote
backend(::Dense{<:Any, <:Knet.Param, <:Knet.Param}) = :Knet
backend(m::Chain) = backend(m.layers)
backend(m::Conv) = :Zygote
backend(m::Tuple{}) = :Any

function backend(m::Tuple)
    if haskey(BACKENDS_CACHE, m)
        BACKENDS_CACHE[m]
    else
        res = :Any
        res_first = backend(m[1])
        res_rest = backend(m[2:end])
        if res_first === :Any
            res = res_rest
        elseif res_rest === :Any
            res = res_first
        elseif res_first === res_rest
            res = res_first
        else
            @error "inconsistent backend detected for $m"
        end
        BACKENDS_CACHE[m] = res
        res
    end
end

function Base.show(io::IO, l::Dense)
  print(io, "Dense(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, "; backend=", backend(l))
  print(io, ")")
end


text:
using Binning
using Flux
using Plots

# binning function stuff -------------------------------

x = range(-10, stop=10, length=200)
plot(x, SmoothStep.(x, 0, 4), label="binning function")

# demonstrate the action of the bin function
y = sin.(x)
b_out = [sin(xᵢ)*SmoothStep(xᵢ, 0, 4) for xᵢ ∈ x]

plot!(x, y, color=:red, label="signal")
plot!(x, b_out, color=:green, label="binned signal")

savefig("../figures/smoothstep.pdf")
savefig("../figures/smoothstep.svg")

#---- test binning layer action --------------------------
x = Array(0:.1:10)
f = sin.(x)
testBinner = DomainBinner(x, 20)
binit = testBinner.b

fnew = testBinner(f)
xnew = getBinCenters(testBinner)

p = plot()
b = sort(vcat(testBinner.b₁,testBinner.b, testBinner.bₑ))
println(size(fnew), size(b))

for i ∈ 2:length(b)
    println(i)
    plot!(p, b[i-1:i], [fnew[i-1], fnew[i-1]], color=:black, fill=(0, 0.5, :black), label="")
end
xlabel!("bin(s)")
ylabel!("count")

savefig("../figures/binPlot.pdf")

#---- test training via integration problem ----------

# set up loss function
function loss()
    bins = sort(vcat(testBinner.b₁,testBinner.b, testBinner.bₑ))
    fnew = testBinner(f)
    Δx = [bins[i]-bins[i-1] for i ∈ 2:length(bins)]
    Î = sum([fnew[i]*Δx[i] for i ∈ 1:length(fnew)])
    I = sum(0.1.*f)
    return (Î-I)^2 # square to make sure loss > 0
end

function binSpreadLoss()
    bins = sort(vcat(testBinner.b₁,testBinner.b, testBinner.bₑ))
    bmin = bins[1]
    bmax = bins[end]
    E = abs(bmax-bmin)-(testBinner.bₑ-testBinner.b₁)
end

λ=1
function composedLoss()
    return loss() + λ*binSpreadLoss()
end

println(loss())
println(binSpreadLoss())
println(composedLoss())

opt = Descent(0.0005)
data = Iterators.repeated((), 200)
ps = params(testBinner)
println(ps)

# train the model!
Flux.train!(composedLoss, ps, data, opt, cb=()->println(composedLoss()))

println(composedLoss())

#p1 = plot()
b = sort(vcat(testBinner.b₁,testBinner.b, testBinner.bₑ))

for i ∈ 2:length(b)
    plot!(p, b[i-1:i], [fnew[i-1], fnew[i-1]], color=:blue, fill=(0, 0.5, :blue), label="")
end
xlabel!("bin(s)")
ylabel!("count")
savefig("../figures/trained_bin_comparison.pdf")
savefig("../figures/trained_bin_comparison.png")


text:
import Dagger: indexes, project

@testset "UnitDomain" begin
    @test domain(1) == UnitDomain()
end

@testset "ArrayDomain" begin
    @test domain(rand(10, 10)) == ArrayDomain((1:10, 1:10))

    d1 = domain(Array{Int}(10, 10))
    @test indexes(d1) == (1:10, 1:10)
    @test size(d1) == (10, 10,)
    @test length(d1) == 100

    d2 = ArrayDomain(5:10, 5:10)
    d3 = ArrayDomain((8:12, 1:8))
    @testset "Domain methods for ArrayDomain" begin
        @test intersect(d2,d2) == d2
        @test intersect(d2,d1) == d2
        @test intersect(d1,d2) == d2
        @test intersect(d2,d3) == ArrayDomain((8:10, 5:8))

        @test project(ArrayDomain((11:25, 21:100)), ArrayDomain((15:20, 30:40))) == ArrayDomain((5:10,10:20))
        # should this cause a BoundsError? vvv
        @test project(ArrayDomain((11:25, 21:35)), ArrayDomain((15:20, 30:40))) == ArrayDomain((5:10,10:20))

        @test ArrayDomain((11:20, 11:20))[ArrayDomain((5:8, 5:10))] == ArrayDomain((15:18, 15:20))
        @test_throws BoundsError ArrayDomain((11:20, 11:20))[ArrayDomain((5:8, 5:11))]

        @test alignfirst(ArrayDomain(11:25, 21:100)) == ArrayDomain((1:15), (1:80))
    end
end


text:
# An eigenproblem described in the LaTeX file eigen.tex.

using FinElt
using FinElt.PlanarPoisson
using Arpack: eigs

const h = 0.05
run(`gmsh -2 -clmax $h -o roundL.msh -format msh22 roundL.geo`)

mesh = read_msh_file("roundL.msh")
essential_bc = [ "Dirichlet" ]
ep = EigenProblem(mesh, essential_bc)
add_bilin_form!(ep, "Omega", grad_dot_grad!,   :LHS)
add_bilin_form!(ep, "Omega", func_times_func!, :RHS)

A, B = assembled_eigenproblem_matrices(ep)

d, vfree, nconv, niter, nmult, resid = eigs(A, B, which=:SM)

v = complete_soln(vfree, ep)

write_pos_file("roundL.pos") do fid
    for k = 1:size(v,2)
        save_warp_nodal_scalar_field(v[:,k], "v$k", fid)
    end
end

run(`gmsh eigen.script`)


text:
using Statistics
import Base.+, Base.-, Base.*, Base./

+(x::PValue, y::PValue) = PValue(x.contribs+y.contribs,
                                 x.expected+y.expected,
                                 x.value+y.value,
                                 max(x.updated,y.updated))

-(x::PValue, y::PValue) = PValue(x.contribs+y.contribs,
                                 x.expected+y.expected,
                                 x.value-y.value,
                                 max(x.updated,y.updated))

*(x::PValue, y::PValue) = PValue(x.contribs+y.contribs,
                                 x.expected+y.expected,
                                 x.value*y.value,
                                 max(x.updated,y.updated))

/(x::PValue, y::PValue) = PValue(x.contribs+y.contribs,
                                 x.expected+y.expected,
                                 x.value/y.value,
                                 max(x.updated,y.updated))

"""

"""
function add_formulas(abaco, df)
    for (tag, name, expression) in eachrow(df)
        if !haskey(abaco.cfg, tag)
            abaco.cfg[tag] = SnapsSetting(nothing, false, nothing)
        end
        setting = abaco.cfg[tag]
        f = formula(setting, name, expression)

        # create a formula state for each node with el.tag==tag
        for el in values(abaco.node)
            if el.tag == tag
                for snap in values(el.snap)
                    snap.outputs[f.output] = FormulaState(false, f.output)
                end
            end
        end
    end
end

function formula(abaco::Context, formula_def)
    setting = abaco.cfg[DEFAULT_TYPE]
    formula(setting, formula_def)
end

function formula(abaco::Context, name, expression, tag)
    if !haskey(abaco.cfg, tag)
        abaco.cfg[tag] = SnapsSetting(nothing, false, abaco.oncompletedefault)
    end
    setting = abaco.cfg[tag]
    f = formula(setting, name, expression)

    # create a formula state for each node with tag==row.tag
    for el in values(abaco.node)
        if el.tag == tag
            for snap in values(el.snap)
                snap.outputs[f.output] = FormulaState(false, f.output)
            end
        end
    end
    f
end

formula(abaco, name, expression) = formula(abaco, name, expression, "")

"""
    formula(setting::SnapsSetting, name, expression)

Add the formula `name` defined by `expression`:
a mathematical expression like `x + y*w`.
"""
formula(setting::SnapsSetting, name, expression) = formula(setting::SnapsSetting, "$name=$expression")

"""
    formula(setting::SnapsSetting, formula_def::String)

Add a formula, with `formula_def` formatted as `"formula_name = expression"`,
where expression is a mathematical expression, like `x + y*w`.
"""
function formula(setting::SnapsSetting, formula_def)
    formula = extractor(formula_def)
    setting.formula[formula.output] = formula
    setup_formula(setting, formula)
end

function setup_formula(setting::SnapsSetting, formula::Formula)
    setting.formula[formula.output] = formula

    # delete already present dependents, otherwise formula redefinition
    # could not work as expected
    for dependent in values(setting.dependents)
        delete!(dependent, formula.output)
    end

    # update the dependents
    for invar in formula.inputs
        if haskey(setting.dependents, invar)
            push!(setting.dependents[invar], formula.output)
        else
            setting.dependents[invar] = Set([formula.output])
        end
    end

    formula
end

function delete_formula(abaco::Context, name::String, tag::String)
    if !haskey(abaco.cfg, tag)
        abaco.cfg[tag] = SnapsSetting(nothing, false, nothing)
    end
    setting = abaco.cfg[tag]
    formula = setting.formula[name]

    for invar in formula.inputs
        delete!(setting.dependents[invar], name)
    end
    delete!(setting.formula, name)
end

"""
    dependents(abaco::Context, tag::String, var::String)

Returns the list of expressions that depends on `var`.
"""
function dependents(abaco::Context, tag::String, var::String)
    result = String[]
    if haskey(abaco.cfg, tag)
        deps = abaco.cfg[tag].dependents
        if haskey(deps, var)
            append!(result, deps[var])
        end
    end

    # TODO: merge glob dependents

    result
end

function dependents(abaco::Context, tag::String, vars)
    result = String[]
    if haskey(abaco.cfg, tag)
        deps = abaco.cfg[tag].dependents
        append!(result, union([haskey(deps, var) ? deps[var] : [] for var in vars]...))
    end

    # TODO?: merge glob dependents

    result
end

#
#     eval(formula::Formula, map::Dict{String,SValue})
#
# Compute the `formula`, looking up the values of
# inputs variables in `map`, and returns the result.
#
# In case of error throws [`Abaco.EvalError`](@ref).
#
function eval(formula::Formula, map::Dict)
    try
        return advance(formula.expr, map)
    catch ex
        # internal error inspection
        # showerror(stdout, ex, catch_backtrace())
        throw(EvalError(string(formula.expr), ex))
    end
end

function advance(s::Symbol, map::Dict)
    var = String(s)
    if haskey(map, var)
        return map[var].value
    else
        throw(UndefVarError(s))
    end
end

function advance(x::Number, map::Dict)
    return x
end

function advance(e::Expr, map::Dict)
    return advance(Val(e.head), e.args, map)
end

function advance(::Val{:call}, args, map::Dict)
    return advance(Val(args[1]), args[2:end], map)
end

function advance(::Val{:ref}, args, map::Dict)
    var = String(args[1])
    idx = args[2]
    if haskey(map, var)
        if idx <= length(map[var].value)
            return map[var].value[idx]
        else
            throw(ValueNotFound(var, idx))
        end
    else
        throw(UndefVarError(var))
    end
end

getsymbol(::Val{x}) where x = x

function advance(fsym::Val, args, map::Dict)
    sym = getsymbol(fsym)
    if hasproperty(Statistics, sym) || sym === :sum
        var = string(args[1])
        if haskey(map, var)
            return PValue(length(map[var].value),
                          map[var].contribs,
                          eval(sym)(map[var].value),
                          nowts())
        else
            throw(UndefVarError(var))
        end
    else
        return eval(sym)([advance(arg, map) for arg in args]...)
    end
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using DelimitedFiles
Random.seed!(233)

function mpi_tensor_transpose(row,col,ncol,val,n,rank,nt)
    require_mpi()
    mpi_tensor_transpose_ = load_op_and_grad("./build/libMPITensor.so","mpi_tensor_transpose", multiple=true)
    row,col,ncol,val,n,rank,nt = convert_to_tensor(Any[row,col,ncol,val,n,rank,nt], [Int32,Int32,Int32,Float64,Int64,Int64,Int64])
    indices, vals = mpi_tensor_transpose_(row,col,ncol,val,n,rank,nt)
end

mpi_init()
row = [0;1] .+ mpi_rank()*2
col = [0;1;2;3;4;5;0;1;2;3;4;5]
ncol = [6;6]
val = Array(1:12) .+ 12 * mpi_rank()

n = 2
rank = mpi_rank()
nt = 6

# TODO: specify your input parameters
u = mpi_tensor_transpose(row,col,ncol,val,n,rank,nt)
sess = Session(); init(sess)

indices, vals = run(sess, u)
sleep(mpi_rank())
println("===================================================")
writedlm(stdout, indices)
writedlm(stdout, vals)

mpi_finalize()


text:
using ComplexRegions
CR = ComplexRegions

using Test
@testset "Utilities" begin
	 @test(CR.scaleto(1im,3im,[0.5,0.75]) ≈ [2.0im,2.5im])
	 @test(CR.scalefrom(1im,3im,[2im,1im,1.5im]) ≈ [0.5,0,0.25])
	 x = CR.realroots(3,16,1)
	 @test( all(@. abs(3x^2 + 16x + 1)<1e-12 ) )
	 z = 2im .+ 3*exp.(2im*pi*[0.8,0.1,0.25])
	 @test( CR.isccw(z...) )
	 @test( CR.intadapt(exp,0,4,1e-13) ≈ (exp(4)-1)  )
	 z = t -> exp(1im*t)
	 @test( CR.fdtangent(z,0) ≈ 1im )
	 @test( CR.fdtangent(z,1) ≈ 1im*exp(1im) )
	 @test( CR.fdtangent(z,0.2) ≈ 1im*exp(0.2im) )
end

@testset "Curves" begin
	f = t -> 2*cos(t) + 1im*sin(t)
	@test( Curve(f,-1,1) isa Curve)
	@test( point(Curve(f,-1,1),0.5) ≈ f(0) )
	f = t -> 2*cos(t) + 3im*sin(t)
	c = ClosedCurve(f,0,2π)
	@test( point(5-3im*c,0.125) ≈ 5-3im*f(pi/4) )
	@test( angle(normal(c,0.75)) ≈ 0.5π )
end

@testset "Circles" begin
	c = Circle(Spherical(1-1im),sqrt(2))
	@test( arclength(c) ≈ 2*sqrt(2)*pi )
	@test( dist(-1+1im,c) ≈ sqrt(2) )
	@test( closest(1+4im,c) ≈ 1+1im*(sqrt(2)-1) )
	@test( isinside(1.5-1im,c) && isoutside(1.5+1im,reverse(c)) )
	@test( isinf(reflect(c.center,c)) )
	@test( reflect(reflect(-1+2im,c),c) ≈ -1+2im )
	@test( all( abs(arg(c,c(t))-t) < 1e-11 for t in 0.1:0.1:1 ) )
	@test( angle(unittangent(c,0.125)) ≈ 0.75π )
	@test( abs(tangent(c,0.125)) ≈ 2π*sqrt(2) )

	c = Circle(1-1im,sqrt(2))
	@test( point(c,.25) ≈ complex(1,sqrt(2)-1) )
	@test( 2/c isa Line )
	c = Circle(1f0,-1im,0)
	@test( c.radius ≈ 1/sqrt(2f0) )
	c = Circle(1,-1im,false)
	@test( all( abs(arg(c,c(t))-t) < 1e-11 for t in 0.1:0.1:1 ) )
	zz = point(5im - c/3im,.23)
	@test( abs(zz-(5im-c.center/3im)) ≈ c.radius/3 )
	@test( tangent(c,0.7) ≈ CR.fdtangent(c,0.7) )
	@test( Circle(1+3im,Polar(1-1im),1.0) isa Line )
end

@testset "Arcs" begin
	a = Arc(exp.(1im*[pi/2,pi/5,0])...)
	zz = 1/sqrt(2)*(1+1im)
	@test( point(a,0.5) ≈ zz )
	@test( all( abs(arg(a,a(t))-t) < 1e-11 for t in 0:0.1:1 ) )
	@test( dist(3im+.5*2im*exp(1im*pi/5),3im+2im*a) ≈ 2*0.5 )
	@test( tangent(a,0.2) ≈ CR.fdtangent(a,0.2) )
	@test( angle(unittangent(a,0.5)) ≈ -0.25π )
	a = (a-2)/1im
	@test( all( abs(arg(a,a(t))-t) < 1e-11 for t in 0:0.1:1 ) )

	b = Arc(-1im,1im,-1)
	@test( point(b,2/3) ≈ 1im )
	@test( point(0.1 - 3im*b,2/3) ≈ 0.1 - 3im*1im )
	@test( closest(5im,b) ≈ 1im )
	@test( closest(2-5im,b+2) ≈ 2-1im )
	@test( all( abs(arg(b,b(t))-t) < 1e-11 for t in 0:0.1:1 ) )

	b = reverse(b)
	@test( point(b,1/3) ≈ 1im )
	@test( all( abs(arg(b,b(t))-t) < 1e-11 for t in 0:0.1:1 ) )
	@test( angle(tangent(b,2/3)) ≈ -0.5π )
end

@testset "Lines" begin
	@test( Line(1,5) isa Line )
	l = Line(1im,direction=1+2im)
	@test( isleft(2im,l) && !isleft(0,l) )
	dz = point(l,.6)-point(l,.1)
	@test( angle(dz) ≈ angle(1+2im) )
	@test( 1/l isa Circle )
	zz = point(5im - l/3im,.23)
	z0 = 5im - l.base/3im
	@test( angle(zz-z0) ≈ angle(tangent(l,0.5)/3im) )
	@test( tangent(l,0.2) ≈ CR.fdtangent(l,0.2) )
	z = l(0.3) + 1im*sign(l.direction)
	@test( dist(z,l) ≈ 1 )
	@test( closest(z,l) ≈ l(0.3) )
	@test( reflect(z,l) ≈ l(0.3) - 1im*sign(l.direction))
	@test( all( abs(arg(l,l(t))-t) < 1e-11 for t in 0:0.1:0.9 ) )
end

@testset "Segments" begin
	s = Segment(1,3.0+5.0im)
	@test( isleft(-1,s) && !isleft(2,s) )
	zz = 2 + 2.5im
	@test( point(2 - 3im*s,0.5) ≈ 2 - 3im*zz)
	@test( closest(4+6im,s) ≈ 3+5im )
	@test( dist(-1,s) ≈ 2 )
	@test( angle(tangent(s,2/3)) ≈ angle(s(0.6)-s(0.1)) )
	@test( tangent(s,0.75) ≈ CR.fdtangent(s,0.75) )
	z = s(0.7) + 1im*sign(s(0.9)-s(0.7))
	@test( closest(z,s) ≈ s(0.7) )
	@test( reflect(z,s) ≈ s(0.7) - (z-s(0.7)) )
	@test( all( abs(arg(s,s(t))-t) < 1e-11 for t in 0:0.1:1 ) )
end

@testset "Rays" begin
	s = Ray(Polar(2,0),pi/2)
	@test( isinf(arclength(s)) )
	@test( isleft(-1im,s) && !isleft(-1im,reverse(s)) )
	@test( real(s(0.23)) ≈ 2 )
	@test( imag(s(0.9)) > imag(s(0.7)) )
	@test( closest(5im,s) ≈ 2+5im )
	@test( all( abs(arg(s,s(t))-t) < 1e-11 for t in 0:0.1:1 ) )
	@test( angle(tangent(s,.1))≈π/2 )
	@test( tangent(s,0.1) ≈ CR.fdtangent(s,0.1) )
	@test( angle(tangent(reverse(s),1))≈-π/2  )
	s = Ray(Spherical(2im),pi,true)
	@test( imag(s(0.5)) ≈ 2 )
	@test( real(s(0.3)) < real(s(0.4)) )
	@test( tangent(s,0.6) ≈ CR.fdtangent(s,0.6) )
	@test( !isleft(4,s) && isleft(-1+3im,s) )
	@test( closest(-4+1im,s) ≈ -4+2im )
	@test( closest(6,s) ≈ 2im )
	@test( all( abs(arg(s,s(t))-t) < 1e-11 for t in 0:0.1:1 ) )
end

@testset "Intersections" begin
	z = intersect(Circle(0,1),Circle(.2+0.5im,1.5))
	@test( all( @. abs(z-0)≈1 ) )
	@test( all( @. abs(z-(.2+0.5im))≈1.5 ) )
	z = intersect(Circle(0,1),Circle(1+1im,1.5))
	@test( all( @. abs(z-0)≈1 ) )
	@test( all( @. abs(z-(1+1im))≈1.5 ) )
	@test( isempty(intersect(Circle(0,1),Circle(.2+0.5im,.1))) )
	@test( isempty(intersect(Circle(0,1),Circle(.2+0.5im,6))) )

	z = intersect(Line(1,direction=1im),Line(-1,direction=1+1im))
	@test( z[1]≈1+2im )
	@test( isempty(intersect(Line(1,direction=1im),Line(-2,direction=1im))) )
	l = Line(2,direction=3+1im)
	@test( intersect(l,l+1e-15) ≈ l )

	z = intersect(Segment(0,1),Segment(.4-1im,.7+2im))
	@test( z[1]≈0.5 )
	z = intersect(Segment(2+3im,3+3im),Segment(-1+3im,4+3im))
	@test( z≈Segment(2+3im,3+3im) )
	z = intersect(2-Segment(2+3im,3+3im),2-Segment(3im,2.4+3im))
	@test( z≈2-Segment(2+3im,2.4+3im) )
	z = intersect(1im*Segment(2+3im,3+3im),1im*Segment(2.7+3im,3.4+3im))
	@test( z≈1im*Segment(2.7+3im,3+3im) )
	z = intersect(Segment(2+3im,3+3im),Segment(2.7+4im,3.4+6im))
	@test( isempty(z) )

	z = intersect(Ray(0,pi/4),Segment(.5-1im,.5+2im))
	@test( z[1]≈0.5+0.5im )
	@test( isempty(intersect(Ray(0,-3pi/4),Segment(.5-1im,.5+2im))) )
	z = intersect(Ray(0,pi/4),Ray(2+2im,pi/4))
	@test( z≈Ray(2+2im,pi/4) )
	z = intersect(Ray(0,pi/6),Ray(-1,-pi/2))
	@test( isempty(z) )

	z = intersect(3+Circle(0,1),3+Line(0.5,0.5+3im))
	@test( z[1]≈(3.5+sqrt(3)/2*1im) || z[1]≈(3.5-sqrt(3)/2*1im) )
	z = intersect(Circle(-2im,2),Line(3im,3))
	@test( isempty(z) )

	z = intersect(3+Ray(0.5+0.1im,pi/2,true),3+Circle(0,1))
	@test( z[1]≈(3.5+sqrt(3)/2*1im) && length(z)==1 )

	z = intersect(3+Circle(0,1),3+Segment(0.5,0.5+3im))
	@test( z[1]≈(3.5+sqrt(3)/2*1im) || z[1]≈(3.5-sqrt(3)/2*1im) )
	z = intersect(Circle(-2im,2),Segment(-1im,.5-1im))
	@test( isempty(z) )

end

@testset "Paths" begin
	S = Segment(1,1im)
	A = Arc(1im,-1+0.5im,-1)
	P = Path([S,A,-S])
	@test( all( point(P,[0,1,1.5,2.5,3]) .≈ [S(0),S(1),A(0.5),-S(0.5),-S(1)] ) )
	Q = 1 - 3im*P
	@test( Q(1.5) ≈ 1 - 3im*A(0.5) )

	P = ClosedPath([S,1im*S,-S,-1im*S])
	@test( vertex(P,3) ≈ -1 )
	@test( arclength(P) ≈ 4*sqrt(2) )
	@test( isa(reverse(P),ClosedPath) )
	@test( all( point(P,[0,1,1.25,2.5,3,4]) .≈ [S(0),S(1),1im*S(.25),-S(0.5),-S(1),S(0)] ) )

end

@testset "Polygons" begin
	s = Segment(2,2im)
	p = Polygon([s,1im*s,-s,-1im*s])
	@test( arclength(p) ≈ 8*sqrt(2) )
	@test( angle(normal(p,1.1+length(p))) ≈ -π/4 )
	@test( winding(p,-0.4+0.5im) == 1 )
	@test( winding(reverse(p),-0.4+0.5im) == -1 )
	@test( winding(p,-4-0.5im) == 0 )
	@test( all( angles(p) .≈ 0.5*pi ) )
	@test( ispositive(p) )
	@test( !ispositive(reverse(p)) )

	p = Polygon([4,4+3im,3im,-2im,6-2im,6])
	@test( arclength(p) ≈ (3+4+5+6+2+2) )
	@test( tangent(p,2.3-length(p)) ≈ -5im )
	@test( winding(p,5-im) == 1 )
	@test( winding(p,-1) == 0 )
	@test( sum(angles(p)) ≈ 4*pi )

	p = CircularPolygon([Arc(1,2+1im,1im),Segment(1im,-1),Arc(-1,-0.5im,-1im),Segment(-1im,1)])
	@test( all(winding(p,z)==1 for z in [1+0.5im,1.7+1im,0,-1+0.05*exp(1im*pi/5),-1im+0.05*exp(1im*0.3*pi)]) )
	@test( all(winding(p,z)==0 for z in [-.999im,0.001-1im,-.999,-1.001,1.001,1.999im]) )
end

@testset "Unbounded polygons" begin
	p = Polygon([5,4+3im,3im,-2im,6-2im,(-pi/2,0)])
	a = angles(p)/pi
	@test( a[6] ≈ -0.5  )
	@test( sum(a.-1) ≈ -2 )

	p = Polygon([(pi/2,pi/2),5,4+3im,3im,-2im,6-2im])
	a = angles(p)/pi
	@test( abs(a[1]) < 1e-10 )
	@test( sum(a.-1) ≈ -2 )
	@test( all(winding(p,z)==1 for z in [1+2im,5-1im,5.5+6im]) )
	@test( all(winding(p,z)==0 for z in [-3im,3+5im,5.5-6im]) )

	p = Polygon([(-pi/2,pi/2),7,4+3im,3im,-2im,6-2im])
	a = angles(p)/pi
	@test( a[1] ≈ -1 )
	@test( sum(a.-1) ≈ -2 )
	@test( all(winding(p,z)==1 for z in [4,7-2im,9]) )
	@test( all(winding(p,z)==0 for z in [4+4im,6+2im,4-3im]) )

	p = Polygon([4+3im,7,(0,0),6-2im,-2im,3im])
	a = angles(p)/pi
	@test( a[3] ≈ -2 )
	@test( sum(a.-1) ≈ -2 )
end

@testset "Möbius" begin
	z = [0,1,2+2im]
	w = [Inf,-1im,-1]
	f = Möbius(z,w)
	@test( all(f.(z).≈w) )
	g = inv(f)
	@test( all(g.(w).≈z) )
	c = Circle(z...)
	l = Line(w[2],w[3])
	@test( f(c)≈l && g(l)≈c )
	d = exterior(c)
	h = halfplane(l)
	@test( f(d)≈!h && g(h)≈!d )
	@test( f(!d)≈h && g(!h)≈d )
	f = Möbius([2,3+im,-5],z)
	g = Möbius(w,[2,3+im,-5])
	h = f ∘ g
	@test( all(h.(w).≈z) )
end


text:
@testset "attach" begin

# tests the attach macro
@attach a = 2
@attach include("attach_include.jl")

# test of basic funtionality
sc = SparkContext(master="local")
txt = parallelize(sc, ["hello", "world"])
rdd = map_partitions(txt, it -> map(func_from_include, it))

@test reduce(rdd, +) == 14

close(sc)

end


text:
# This file was generated, do not modify it. # hide
DataFrames.describe(features)


text:
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

# D, Δx = 1.0, 0.35
# Γ = meshsphere(D, Δx)

Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
X = raviartthomas(Γ)
#Δt, Nt = 0.08, 400
Δt, Nt = 0.6, 200
T = timebasisc0d1(Δt, Nt)
U = timebasiscxd0(Δt, Nt)

V = X ⊗ T
W = X ⊗ U

width, delay, scaling = 8.0, 12.0, 1.0
gaussian = creategaussian(width, delay, scaling)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, derive(gaussian), 1.0)

@hilbertspace j; @hilbertspace j′
T = MWSingleLayerTDIO(1.0,-1.0,-1.0,2,0)
tdefie = @discretise T[j′,j] == -1E[j′]   j∈V  j′∈W
xefie = solve(tdefie)

Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
ω = collect(ω0 .+ (0:Nt-1)*Δω)
_, i1 = findmin(abs.(ω.-1.0))

ω1 = ω[i1]
ue = Xefie[:,i1]/ fouriertransform(gaussian)(ω1)


text:
# ellipse to rotated and shifted ellipse, initial and final landmark positions observed
using Revise
using Polynomials.PolyCompat
using BridgeLandmarks
const BL = BridgeLandmarks
using Distributions
using LinearAlgebra
using Random
using FileIO
using JLD2
using DelimitedFiles
using Distances

Random.seed!(9)
workdir = @__DIR__
cd(workdir)
outdir = workdir

d1 = readdlm("dataset1.csv",',')
d2 = readdlm("dataset2.csv",',')
n = 12

center = PointF(40.0, 20.0)
xobs0 = map(i->PointF(d1[i,:]) - center, 1:n)
xobsT = map(i->PointF(d2[i,:]) - center, 1:n)

# read data
# dat = load("data_exp1.jld2")
# xobs0 = dat["xobs0"]
# xobsT = dat["xobsT"]

# set pars
p_ms = Pars_ms(δmom=0.1, cinit=0.2, γinit=2.0, σobs = 0.01)

# for ahs adjust domain bounds
#p_ahs = Pars_ahs(δmom=0.1, cinit=0.02, γinit=0.2, db=[30.0,25.0],stdev=5.0, σobs = 0.01)
ainit = 0.5*mean(norm.(diff(xobs0)))  #pairwise(dist, X, Y, dims=2)
ups = [:innov, :mala_mom, :parameter]#, :matching]
ups = [:innov, :rmrw_mom, :parameter]#, :matching]

# run algorithm
@time landmarksmatching(xobs0,xobsT; ITER=200,pars=p_ms, outdir=outdir, ainit=ainit, updatescheme=ups)
#@time landmarksmatching(xobs0,xobsT; ITER=25,pars=p_ahs, outdir=outdir, ainit=ainit, updatescheme=ups)

plotlandmarksmatching(outdir)


text:
"""
Usage:

using HOI4CombatSimulator
HOI4CombatSimulator.quick.init()
using HOI4CombatSimulator.quick: Infantry, Artillery, Tank, Motor, SPA,
    inf7art2, inf10, tank5mot2spa2, div72, div10, div522
"""
module quick

using ..HOI4CombatSimulator: get_unit_list, DivisionTemplate, Division

function init()
    global unit_list = get_unit_list()
    # there're ("Support", "Artillery") and ("Infantry", "Artillery")
    global unit_dict = Dict((unit.CategoryHead, unit.Unit) => unit for unit in unit_list)

    # unit type (base: Light)
    global Infantry = unit_dict[("Infantry", "Infantry")]
    global Artillery = unit_dict[("Infantry", "Artillery")]
    global Tank = unit_dict[("Tanks", "Light TD")]
    global Motor = unit_dict[("Mobile Battalions", "Motorized Infantry")]
    global SPA = unit_dict[("Tanks", "Light SP Artillery")]

    # common 20 width template
    global inf7art2 = DivisionTemplate((Infantry, 7), (Artillery, 2))
    global inf10 = DivisionTemplate((Infantry, 10))
    global tank5mot2spa2 = DivisionTemplate((Tank, 5), (Motor, 2), (SPA, 2))

    # quick Division example
    global div72 = Division(inf7art2)
    global div10 = Division(inf10)
    global div522 = Division(tank5mot2spa2)

end

end


text:
### MultipleLines

using Statistics: mean, std, var

struct ErrorBarLines
    x
    ys
end

@recipe function f(lines::ErrorBarLines; nstd=1, std_warn=nothing)
    @unpack x, ys = lines

    y = mean(ys)
    y_s = std(ys)
    dy = nstd * y_s

    # Add a series for an error band
    @series begin
        seriestype := :path
        # Ignore series in legend and color cycling
        primary := false
        linecolor := nothing
        fillcolor := :lightgray
        fillalpha := 0.5
        fillrange := y .- dy
        # ensure no markers are shown for the error band
        markershape := :none
        # Return series data
        x, y + dy
    end

    # Highlight `s` larger than `σ_warn` via :x as markers
    # NOTE: this is NOT properly supported by PGFPlotsX
    if !isnothing(std_warn)
        s = get(plotattributes, :markershape, :none)
        markershape := ifelse.(y_s .> std_warn, :xcross, s)
    end

    x, y
end

### TwoYAxesLines

using Plots: mm
import Plots: plot

struct TwoYAxesLines
    x
    y1
    y2
end

function plot(
    lines::TwoYAxesLines;
    color1=1,
    color2=2,
    xlabel=nothing,
    ylabel1=nothing,
    ylabel2=nothing,
    kwargs...
)
    @unpack x, y1, y2 = lines
    right_margin = isnothing(ylabel2) ? :match : 13mm
    p = plot(x, y1; color=color1, xlabel=xlabel, ylabel=ylabel1, right_margin=right_margin, kwargs...)
    plot!(twinx(p), x, y2; color=color2, ylabel=ylabel2, kwargs...)
    return p
end


text:
"""
    bubble(x, y)

Creates an `EChart` scatterplot, with additional dimension represented by circle size.

## Methods
```julia
bubble(x::AbstractVector{<:Union{Missing, Real}}, y::AbstractVector{<:Union{Missing, Real}},
    size::AbstractVector{<:Union{Missing, Real}})
bubble(df::AbstractDataFrame, x::Symbol, y::Symbol, size::Symbol)
bubble(df::AbstractDataFrame, x::Symbol, y::Symbol, size::Symbol, group::Symbol)
```

## Arguments
* `legend::Bool` : display legend?
* `scale::Bool = false` : show full Y-axis or truncated
* `large::Bool = true` : minimize overplotting
* `largeThreshold::Int = 2000` : number of points before overplotting optimization occurs
* `bubblesize::Real = 50` : maximum size of bubbles (see notes)
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

The displayed size of the bubble is normalized within the function, as the square root of the size values. Using the square root of the value maintains the proper visual perception of the difference. The overall circle size can be modified by the `bubblesize` keyword argument, which is roughly the area in pixels of the largest circle.

## Examples
```julia
Random.seed!(13579)
xval = rand() .* (0:1000:30000)
yval = 60 .+ (rand(31) * 20)
sizeval = shuffle!(6000 .* (rand(31) * 50))
sp = bubble(xval, yval, sizeval, scale = true)
```
"""
function bubble(x::AbstractVector{<:Union{Missing, Real}},
				y::AbstractVector{<:Union{Missing, Real}},
				size::AbstractVector{<:Union{Missing, Real}};
		legend::Bool = false,
		scale::Bool = false,
		large::Bool = true,
		largeThreshold::Int = 2000,
		bubblesize::Real = 50,
		kwargs...)

		#Build base via scatter
		ec = scatter(x, y;
					legend = legend,
					scale = scale,
					large = large,
					largeThreshold = largeThreshold,
					kwargs...)

		#Append size into existing data structure
		[push!(x, y) for (x,y) in zip(ec.series[1].data, size)]

		normalized = maximum(sqrt.(size))

		[x.symbolSize =  JSFunction("""function (data) {return $bubblesize * Math.sqrt(data[2]) / ($normalized);}""") for x in ec.series]

		return ec

end

#dataframe, single series
function bubble(df::AbstractDataFrame, x::Symbol, y::Symbol, size::Symbol;
		legend::Bool = false,
		scale::Bool = false,
		large::Bool = true,
		largeThreshold::Int = 2000,
		bubblesize::Real = 50,
		kwargs...)

		#Not currently handling NAs
		ec = bubble(df[x], df[y], df[size],
					legend = legend,
					scale = scale,
					large = large,
					largeThreshold = largeThreshold,
					bubblesize = bubblesize,
					kwargs...)

		#Name axes since we know them
		yaxis!(ec, name = string(y))
		xaxis!(ec, name = string(x))

end

#dataframe with group: uses single series df method
function bubble(df::AbstractDataFrame, x::Symbol, y::Symbol, size::Symbol, group::Symbol;
			legend::Bool = true,
			scale::Bool = false,
			large::Bool = true,
			largeThreshold::Int = 2000,
			kwargs...)

			#Create grouped df
			subdf = groupby(df, group)

			#Get number of groups
			numgroups = length(subdf)

			#Initialize with single series
			ec = bubble(subdf[1], x, y, size; legend = legend, scale = scale, kwargs...)

			#Append remaining series, put symbolSize
			for i in 2:numgroups
				push!(ec.series, XYSeries(_type = "scatter", data = arrayofarray(subdf[i][x], subdf[i][y], subdf[i][size])))

				ec.series[i].symbolSize = ec.series[1].symbolSize
			end

			#Add series names based on levels of grouped df
			seriesnames!(ec, [unique(subdf[x][group])[1] for x in 1:numgroups])

			#Add legend if desired
			legend == true ? legend!(ec) : nothing

			#Enable optimization for lots of data
			[x.large = large for x in ec.series]
			[x.largeThreshold = largeThreshold for x in ec.series]

			return ec

end


text:
push!(LOAD_PATH, "./")

using Rando
using Stiffness
using BenchmarkTools

function prSubDiag(m)

    for i=0:5
        for j=0:i
            print(m[div(i*(i+1), 2) + j+1], " ")
        end

        println()
    end
end

function RandomTriangle!(R::RandoData, x::Array{Float64, 1}, y::Array{Float64, 1})
    for i in 1:3
        x[i] = Rando.fv!(R, 10.)
    end

    for i in 1:3
        y[i] = Rando.fv!(R, 10.)
    end
end

const S = StiffnessData()
const x = Float64[0., 1., 0.]
const y = Float64[0., 0., 1.]
const m = zeros(21)
const R = RandoData()
const ntri = 1_000_000

print("\nVerify that, on the reference element, we are coherent with sage ")
println("(see ../sage/):\n")
Stiffness.op!(S, x, y, m)
prSubDiag(m)

print("\nWe must get the same result if we dilate the triangle:\n")
x[:] *= 2.0
y[:] *= 2.0
Stiffness.op!(S, x, y, m)
prSubDiag(m)

println("\nNow, start the benchmark:")
println(ntri, " triangles.")

io = IOContext(stdout, :compact => false)
bench = @benchmarkable Stiffness.op!($S, $x, $y, $m) setup = (RandomTriangle!($R, $x, $y)) samples = ntri
bench_res = run(bench)

show(io, bench_res)

println("\n")
timeByTr = BenchmarkTools.median(bench_res).time * 1e-9

open("RunningOn" * gethostname(), "w") do f
    write(f, string(timeByTr), "\n")
end

println("end.")

# to profile the code, uncomment:
# ntri1 = 100000
# @profile for t=1:ntri1
#     RandomTriangle!(R, x, y)
#     Stiffness.op!(S, x, y, m)
# end
# Profile.print(format=:flat)

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
global MINSTEP  = -0.6;
global MAXVALUE =  0.6;
global AGCGAIN  =  1.0;

"""
    setagc(;gain=1.0, maxvalue=0.6, minstep=-0.6)

Set the automatic-gain-control module's parameters.
# Optional Arguments
- `gain    `: inital speech signal's gain value
- `maxvalue`: specified maximum value for speech signal, maxvalue ∈ (0, +Inf)
- `minstep `: minimum value to change the gain, minstep ∈ (-1, 0)
"""
function setagc(;gain=1.0, maxvalue=0.6, minstep=-0.6)
    global MINSTEP;
    global MAXVALUE;
    global AGCGAIN;
    MINSTEP  = (-1.0 < minstep < 0.0) ? minstep : -0.6
    MAXVALUE = ( 0.0 < maxvalue ) ? maxvalue : 0.6
    MAXVALUE = 1/MAXVALUE
    AGCGAIN  = gain
    return nothing
end

"""
    onlineagc!(wav::Array)

Streaming automatic-gain-control module for speech signal.
# Arguments
- `wav`: one frame of speech signal
"""
function onlineagc!(wav::Array)
    global MINSTEP;
    global MAXVALUE;
    global AGCGAIN;
    fmax = maximum(abs.(wav));
    step = 1.0 - fmax*MAXVALUE;
    step = step * abs(step);
    step = max(MINSTEP, min(step, 0.0));
    AGCGAIN = 0.8*AGCGAIN + 0.2*(1 + step);
    wav .*= AGCGAIN;
    return nothing
end

"""
    onlineagc(wav::Array) -> Array

Streaming automatic-gain-control module for speech signal.
# Arguments
- `wav`: one frame of speech signal
"""
function onlineagc(wav::Array)
    data = deepcopy(wav)
    onlineagc!(data)
    return data
end

max_stars_repo_name: sonosole/AGC

max_issues_repo_name: sonosole/AutoGainCtrl

max_forks_repo_name: sonosole/AutoGainCtrl


text:
import Base.in

export SetConstraint
export in

### Set constraint
mutable struct SetConstraint{T <: ContiSet} <: Constraint
  head::Symbol
  id_hash::UInt64
  lhs::AbstractExpr
  rhs::T
  size::Tuple{Int, Int}
  dual::ValOrNothing
end

function SetConstraint(lhs::Variable, rhs::ContiSet)
  id_hash = hash((lhs, rhs, :(==)))
  return SetConstraint(:(set), id_hash, lhs, rhs, (1, 1), nothing)
end

function get_cv(constraint::SetConstraint)
  return get_cv(constraint.lhs)
end

function get_sv(constraint::SetConstraint)
  return get_sv(constraint.lhs)
end

in(lhs::Variable, rhs::ContiSet) = SetConstraint(lhs, rhs)

# the fenchel conjugate of a set indicator function
fenchel(setconstraint::SetConstraint, w::AbstractArray) = fenchel(setconstraint.rhs, w)

# the prox of a set indicator function
prox(setconstraint::SetConstraint, w::AbstractArray) = prox(setconstraint.rhs, w)


text:
using QuasiArrays, Test

@testset "AbstractQuasiArray" begin
    A = QuasiArray(rand(5,4,3), (range(0;stop=1,length=5), Base.OneTo(4), [2,3,6]))

    @testset "Convert" begin
        @test convert(typeof(A),A) === convert(AbstractQuasiArray{Float64},A) ===
                convert(AbstractQuasiArray{Float64,3},A) === A
        @test Array(A) == Array{Float64}(A) == Array{Float64,3}(A) ==
                convert(Array,A) == convert(Array{Float64},A) ==
                convert(Array{Float64,3},A) ==
                convert(AbstractArray,A) == convert(AbstractArray{Float64},A) ==
                convert(AbstractArray{Float64,3},A) == parent(A)
        @test QuasiArray(A) == QuasiArray{Float64}(A) == QuasiArray{Float64,3}(A) ==
                QuasiArray{Float64,3,typeof(A.axes)}(A) == copy(A) == A
        v = QuasiArray(1:5, (range(0;stop=1,length=5),))
        @test Vector(v) == convert(AbstractVector,v) == Vector{Int}(v) ==
                    convert(AbstractVector{Int},v) == parent(v)
        @test QuasiVector(v) == v
        M = QuasiArray(rand(5,4), (range(0;stop=1,length=5), Base.OneTo(4)))
        @test Matrix(M) == convert(AbstractMatrix, M) == Matrix{Float64}(M) ==
                convert(AbstractMatrix{Float64}, M) == parent(M)
        @test QuasiMatrix(M) == M
    end

    @testset "QuasiArray basics" begin
        @test_throws ArgumentError QuasiArray(rand(5,4,3), (range(0;stop=1,length=6), Base.OneTo(4), [1,3,4]))
        @test_throws MethodError QuasiArray(rand(5,4,3), Inclusion.((range(0;stop=1,length=6), Base.OneTo(4), [1,3,4])))
        @test axes(A) == (Inclusion(range(0;stop=1,length=5)), Base.OneTo(4), Inclusion([2,3,6]))
        @test A[0.25,2,6] == A.parent[2,2,3]
        @test_throws BoundsError A[0.1,2,6]
    end

    @testset "Bounds checking" begin
        @test QuasiArrays.checkbounds(Bool, A, 0, 1, 2) == QuasiArrays.checkbounds(Bool, A, 0.0, 1, 2) == true
        @test QuasiArrays.checkbounds(Bool, A, 1, 4, 3) == true
        @test QuasiArrays.checkbounds(Bool, A, -1, 1, 1) == false
        @test QuasiArrays.checkbounds(Bool, A, 1, -1, 1) == false
        @test QuasiArrays.checkbounds(Bool, A, 1, 1, -1) == false
        @test QuasiArrays.checkbounds(Bool, A, 6, 4, 3) == false
        @test QuasiArrays.checkbounds(Bool, A, 5, 5, 3) == false
        @test QuasiArrays.checkbounds(Bool, A, 5, 4, 4) == false
        @test QuasiArrays.checkbounds(Bool, A, 0, 2, 2, 1) == true  # extra indices
        @test QuasiArrays.checkbounds(Bool, A, 2, 2, 2, 2) == false
        @test QuasiArrays.checkbounds(Bool, A, 1, 1)  == false
        @test QuasiArrays.checkbounds(Bool, A, 1, 12) == false
        @test QuasiArrays.checkbounds(Bool, A, 5, 12) == false
        @test QuasiArrays.checkbounds(Bool, A, 1, 13) == false
        @test QuasiArrays.checkbounds(Bool, A, 6, 12) == false
    end

    @testset "vector indices" begin
        @test QuasiArrays.checkbounds(Bool, A, 0:0.25:1, 1:4, [2,3,6]) == true
        @test QuasiArrays.checkbounds(Bool, A, 0:5, 1:4, 1:3) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:5, 0:4, 1:3) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:5, 1:4, 0:3) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:6, 1:4, 1:3) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:5, 1:5, 1:3) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:5, 1:4, 1:4) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:60) == true
        @test QuasiArrays.checkbounds(Bool, A, 1:61) == false
        @test QuasiArrays.checkbounds(Bool, A, 0.25, 2, 2, 1:1) == true  # extra indices
        @test QuasiArrays.checkbounds(Bool, A, 2, 2, 2, 1:2) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:5, 1:4) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:5, 1:12) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:5, 1:13) == false
        @test QuasiArrays.checkbounds(Bool, A, 1:6, 1:12) == false
    end

    @testset "inferred axes/size" begin
        @inferred axes(A)
        @test @inferred(size(A)) == (5,4,3)
        @test @inferred(size(A, 2)) == 4
    end

    @testset "isinteger and isreal" begin
        v = QuasiArray(1:5, (range(0;stop=1,length=5),))
        @test all(isinteger, QuasiDiagonal(v))
        @test isreal(QuasiDiagonal(v))
    end

    @testset "unary ops" begin
        A = QuasiDiagonal(QuasiArray(1:5, (range(0;stop=1,length=5),)))
        @test +(A) == A
        @test *(A) == A
    end

    @testset "ndims and friends" begin
        A = QuasiDiagonal(QuasiArray(1:5, (range(0;stop=1,length=5),)))
        @test ndims(A) == 2
        @test ndims(QuasiDiagonal{Float64}) == 2
    end

    @testset "empty" begin
        v = QuasiArray([1, 2, 3],(2:4,))
        v2 = empty(v)
        v3 = empty(v, Float64)
        @test !isempty(v)
    end

    @testset "QuasiArray indexing" begin
        v = QuasiArray([1, 2, 3],(0:0.5:1,))
        @test axes(v) == axes(v[:]) == axes(v[Inclusion(0:0.5:1)]) == (Inclusion(0:0.5:1),)
        @test v[0.5] == v[:][0.5] == v[Inclusion(0:0.5:1)][0.5] == 2
    end
end