text:
@doc """
`CoB` is an abstract **c**hange **o**f **b**asis super-type
"""->
abstract CoB

# ------------------------------------------------------------
# Frequency to wavelets

@doc """
`Freq2Wave` is a change of basis type between frequency samples and wavelets.

There are sub types for wavelets with and without boundary correction.
To initialize a `Freq2Wave` type, run

	Freq2Wave(samples, wavename::String, J::Int, B; ...)

- `samples` are the sampling locations as a vector for 1D and M-by-2 matrix for 2D.
- `wave` is the name of the wavelet; see documentation for possibilities.
- `J` is the scale of the wavelet transform to reconstruct.
- `B` is the bandwidth of the samples; only needed if `samples` are non-uniform.
- Optional arguments (if needed) are passed to the functions computing Fourier transforms.
"""->
abstract Freq2Wave <: CoB

abstract Freq2Wave1D <: Freq2Wave
abstract Freq2Wave2D <: Freq2Wave
# The 2D Freq2BoundaryWave's need entries the 1D's do not, so don't
# use the dimension as a TypePar.

# No boundary correction
immutable Freq2NoBoundaryWave1D <: Freq2Wave1D
	# Sampling
	internal::Vector{Complex{Float64}}
	weights::Nullable{ Vector{Complex{Float64}} }

	# Reconstruction
	J::Int64
	wavename::AbstractString

	# Multiplication
	NFFT::NFFT.NFFTPlan{1,0,Float64}

	tmpMulVec::Vector{Complex{Float64}}
end

immutable Freq2NoBoundaryWave2D <: Freq2Wave2D
	internal::Vector{ Vector{Complex{Float64}} }
	weights::Nullable{ Vector{Complex{Float64}} }

	J::Int64
	wavename::AbstractString

	NFFT::NFFT.NFFTPlan{2,0,Float64}

	tmpMulVec::Vector{Complex{Float64}}
end

# Boundary correction
immutable Freq2BoundaryWave1D <: Freq2Wave1D
	internal::Vector{Complex{Float64}}
	weights::Nullable{Vector{Complex{Float64}}}

	J::Int64
	wavename::AbstractString

	NFFT::NFFT.NFFTPlan{1,0,Float64}

	left::Matrix{Complex{Float64}}
	right::Matrix{Complex{Float64}}

	tmpMulVec::Vector{Complex{Float64}}
end

immutable Freq2BoundaryWave2D <: Freq2Wave2D
	internal::Vector{ Vector{Complex{Float64}} }
	weights::Nullable{Vector{Complex{Float64}}}

	J::Int64
	wavename::AbstractString

	NFFT::NFFT.NFFTPlan{2,0,Float64}
	NFFTx::NFFT.NFFTPlan{2,1,Float64}
	NFFTy::NFFT.NFFTPlan{2,2,Float64}

	#= left::Dict{ Symbol, Matrix{Complex{Float64}} } =#
	left::Vector{ Matrix{Complex{Float64}} }
	right::Vector{ Matrix{Complex{Float64}} }

	tmpMulVec::Matrix{Complex{Float64}}
	tmpMulVecT::Matrix{Complex{Float64}} # Serves as the transpose of tmpMulVec
	tmpMulcVec::Vector{Complex{Float64}}
	weigthedVec::Vector{Complex{Float64}}
end

# ------------------------------------------------------------------------
# Constructors

function Freq2NoBoundaryWave1D(internal, weights, J, wavename, NFFT)
	tmpMulVec = Array{Complex{Float64}}( NFFT.M )
	Freq2NoBoundaryWave1D( internal, weights, J, wavename, NFFT, tmpMulVec )
end

function Freq2NoBoundaryWave2D(internal, weights, J, wavename, NFFT)
	tmpMulVec = Array{Complex{Float64}}( NFFT.M )
	Freq2NoBoundaryWave2D( internal, weights, J, wavename, NFFT, tmpMulVec )
end

function Freq2BoundaryWave1D(internal, weights, J, wavename, NFFT, left, right)
	tmpMulVec = Array{Complex{Float64}}( NFFT.M )
	Freq2BoundaryWave1D( internal, weights, J, wavename, NFFT, left, right, tmpMulVec )
end

function Freq2BoundaryWave2D(internal, weights, J, wavename, NFFT, left, right)
	tmpMulVec = similar(left[1])
	tmpMulVecT = tmpMulVec.'

	tmpMulcVec = Array{Complex{Float64}}( NFFT.M )
	weigthedVec = similar(tmpMulcVec)

	vm = van_moment(wavename)
	NFFTx = NFFTPlan( NFFT.x[1,:], 1, (NFFT.N[1],vm) )
	NFFTy = NFFTPlan( NFFT.x[2,:], 2, (vm,NFFT.N[2]) )

	Freq2BoundaryWave2D( internal, weights, J, wavename, NFFT,
	NFFTx, NFFTy, left, right, tmpMulVec, tmpMulVecT, tmpMulcVec, weigthedVec )
end

# ------------------------------------------------------------------------
# Load methods

include("freq2wave.jl")


text:
using Sobol
using LatinHypercubeSampling

abstract type SamplingAlgorithm end

"""
GridSample{T}

T is the step dx for lb:dx:ub
"""
struct GridSample{T} <: SamplingAlgorithm
    dx::T
end

struct UniformSample <: SamplingAlgorithm end
struct SobolSample <: SamplingAlgorithm end
struct LatinHypercubeSample <: SamplingAlgorithm end

"""
LowDiscrepancySample{T}

T is the base for the sequence
"""
struct LowDiscrepancySample{T} <: SamplingAlgorithm
    base::T
end

struct RandomSample <: SamplingAlgorithm end

struct KroneckerSample{A,B} <: SamplingAlgorithm
    alpha::A
    s0::B
end

struct GoldenSample <: SamplingAlgorithm end

struct SectionSample{T} <: SamplingAlgorithm
    x0::Vector{T}
    sa::SamplingAlgorithm
end

"""
sample(n,lb,ub,S::GridSample)

Returns a tuple containing numbers in a grid.

"""
function sample(n,lb,ub,S::GridSample)
    dx = S.dx
    if lb isa Number
        return vec(rand(lb:S.dx:ub,n))
    else
        d = length(lb)
        x = [[rand(lb[j]:dx[j]:ub[j]) for j = 1:d] for i in 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,lb,ub,::UniformRandom)

Returns a Tuple containing uniform random numbers.
"""
function sample(n,lb,ub,::UniformSample)
    if lb isa Number
        return rand(Uniform(lb,ub),n)
    else
        d = length(lb)
        x = [[rand(Uniform(lb[j],ub[j])) for j in 1:d] for i in 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,lb,ub,::SobolSampling)

Returns a Tuple containing Sobol sequences.
"""
function sample(n,lb,ub,::SobolSample)
    s = SobolSeq(lb,ub)
    skip(s,n)
    if lb isa Number
        return [next!(s)[1] for i = 1:n]
    else
        return Tuple.([next!(s) for i = 1:n])
    end
end

"""
sample(n,lb,ub,::LatinHypercube)

Returns a Tuple containing LatinHypercube sequences.
"""
function sample(n,lb,ub,::LatinHypercubeSample)
    d = length(lb)
    if lb isa Number
        x = vec(LHCoptim(n,d,1)[1])
        # x∈[0,n], so affine transform
        return @. (ub-lb) * x/(n) + lb
    else
        lib_out = float(LHCoptim(n,d,1)[1])
        # x∈[0,n], so affine transform column-wise
        @inbounds for c = 1:d
            lib_out[:,c] = (ub[c]-lb[c])*lib_out[:,c]/n .+ lb[c]
        end
        x = [lib_out[i,:] for i = 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,lb,ub,S::LowDiscrepancySample)

Low discrepancy sample:
- Dimension 1: Van der Corput sequence
- Dimension > 1: Halton sequence
If dimension d > 1, all bases must be coprime with each other.
"""
function sample(n,lb,ub,S::LowDiscrepancySample)
    d = length(lb)
    if d == 1
        #Van der Corput
        b = S.base
        x = zeros(Float32,n)
        for i = 1:n
            expansion = digits(i,base = b)
            L = length(expansion)
            val = zero(Float32)
            for k = 1:L
                val += expansion[k]*float(b)^(-(k-1)-1)
            end
            x[i] = val
        end
        # It is always defined on the unit interval, resizing:
        return @. (ub-lb) * x + lb
    else
        #Halton sequence
        x = zeros(Float32,n,d)
        for j = 1:d
            b = S.base[j]
            for i = 1:n
                val = zero(Float32)
                expansion = digits(i, base = b)
                L = length(expansion)
                val = zero(Float32)
                for k = 1:L
                    val += expansion[k]*float(b)^(-(k-1)-1)
                end
                x[i,j] = val
            end
        end
        #Resizing
        # x∈[0,1], so affine transform column-wise
        @inbounds for c = 1:d
            x[:,c] = (ub[c]-lb[c])*x[:,c] .+ lb[c]
        end

        y = [x[i,:] for i = 1:n]
        return Tuple.(y)
    end
end

"""
sample(n,d,D::Distribution)

Returns a Tuple containing numbers distributed as D
"""
function sample(n,d,D::Distribution)
    if d == 1
        return rand(D,n)
    else
        x = [[rand(D) for j in 1:d] for i in 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,d,K::KroneckerSample)

Returns a Tuple containing numbers following the Kronecker sample
"""
function sample(n,lb,ub,K::KroneckerSample)
    d = length(lb)
    alpha = K.alpha
    s0 = K.s0
    if d == 1
        x = zeros(n)
        @inbounds for i = 1:n
            x[i] = (s0+i*alpha)%1
        end
        return @. (ub-lb) * x + lb
    else
        x = zeros(n,d)
        @inbounds for j = 1:d
            for i = 1:n
                x[i,j] = (s0[j] + i*alpha[j])%i
            end
        end
        #Resizing
        # x∈[0,1], so affine transform column-wise
        @inbounds for c = 1:d
            x[:,c] = (ub[c]-lb[c])*x[:,c] .+ lb[c]
        end

        y = [x[i,:] for i = 1:n]
        return Tuple.(y)
    end
end

function sample(n,lb,ub,G::GoldenSample)
    d = length(lb)
    if d == 1
        x = zeros(n)
        g = (sqrt(5)+1)/2
        a = 1.0/g
        for i = 1:n
            x[i] = (0.5+a*i)%1
        end
        return @. (ub-lb) * x + lb
    else
        x = zeros(n,d)
        for j = 1:d
            #Approximate solution of x^(d+1) = x + 1, a simple newton is good enough
            y = 2.0
            for s = 1:10
                g = (1+y)^(1/(j+1))
            end
            a = 1.0/g
            for i = 1:n
                x[i,j] = (0.5+a*i)%1
            end
        end
        @inbounds for c = 1:d
            x[:,c] = (ub[c]-lb[c])*x[:,c] .+ lb[c]
        end
        y = [x[i,:] for i = 1:n]
        return Tuple.(y)
    end
end

fixed_dimensions(
        section_sampler::SectionSample)::Vector{Int64} = findall(
    x->x == false, isnan.(section_sampler.x0))

free_dimensions(
        section_sampler::SectionSample)::Vector{Int64} = findall(
    x->x == true, isnan.(section_sampler.x0))

"""
sample(n,d,K::SectionSample)

Returns Tuples constrained to a section.

In surrogate-based identification and control, optimization can alternate between unconstrained sampling in the full-dimensional parameter space, and sampling constrained on specific sections (e.g. a planes in a 3D volume),

A SectionSampler allows sampling and optimizing on a subset of 'free' dimensions while keeping 'fixed' ones constrained.
The sampler is defined as in e.g.

section_sampler_y_is_10 = SectionSample([NaN64, NaN64, 10.0, 10.0], Surrogates.UniformSample())

where the first argument is a Vector{T} in which numbers are fixed coordinates and `NaN`s correspond to free dimensions, and the second argument is a SamplingAlgorithm which is used to sample in the free dimensions.
"""
function sample(n,lb,ub,section_sampler::SectionSample)
    if lb isa Number
        return rand(Uniform(lb,ub),n)
    else
        d_free = Surrogates.free_dimensions(section_sampler)
        new_samples = sample(n, lb[d_free], ub[d_free], section_sampler.sa)
        out_as_vec = repeat(section_sampler.x0', n, 1)
        for y in 1:size(out_as_vec,1)
            for xi in 1:length(d_free)
                out_as_vec[y,xi] = new_samples[y][xi]
            end
        end
        [Tuple(out_as_vec[y,:]) for y in 1:size(out_as_vec,1)]
    end
end

ext: jl


text:
# K-mer
# =====
#
# Compact k-mer sequence type.
#
# A Kmer is a sequence ≤ 32nt, without any 'N's, packed in a single 64 bit
# value.  While BioSequence is an efficient general-purpose sequence
# representation, Kmer is useful for applications like assembly, k-mer counting,
# k-mer based quantification in RNA-Seq, etc that rely on manipulating many
# short sequences as efficiently (space and time) as possible.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Representation
# --------------
#
# Four kinds of nucleotides are encoded as follows:
#
#   nucleotide | binary
#   ---------- | ------
#       A      |   00
#       C      |   01
#       G      |   10
#     T / U    |   11
#
# NucleicAcids are filled from MSBs to LSBs and right-aligned so that all k-mers
# are lexicographically ordered. For example, the memory layout of "TACG" is:
#   64-bit: 0b 00 00 … 00 11 00 01 10
#    4-mer:                T  A  C  G

primitive type Kmer{T<:NucleicAcid, K} <: Sequence 64 end

const DNAKmer{K} = Kmer{DNA, K}
const RNAKmer{K} = Kmer{RNA, K}
const DNACodon = DNAKmer{3}
const RNACodon = RNAKmer{3}

function Kmer(nts::T...) where {T<:NucleicAcid}
    return make_kmer(nts)
end

function DNACodon(x::DNA, y::DNA, z::DNA)
    return make_kmer((x, y, z))
end

function RNACodon(x::RNA, y::RNA, z::RNA)
    return make_kmer((x, y, z))
end

# Conversion
# ----------

function DNAKmer{K}(x::UInt64) where {K}
    mask = ~UInt64(0) >> (64 - 2K)
    return reinterpret(DNAKmer{K}, x & mask)
end

function RNAKmer{K}(x::UInt64) where {K}
    mask = ~UInt64(0) >> (64 - 2K)
    return reinterpret(RNAKmer{K}, x & mask)
end

UInt64(x::Kmer) = reinterpret(UInt64, x)
Base.convert(::Type{UInt64}, x::Kmer) = reinterpret(UInt64, x)

DNAKmer{K}(x::RNAKmer{K}) where {K} = reinterpret(DNAKmer{K}, x)
RNAKmer{K}(x::DNAKmer{K}) where {K} = reinterpret(RNAKmer{K}, x)

function Kmer{T,K}(seq::AbstractString) where {T,K}
    return make_kmer(Kmer{T,K}, seq)
end

function Kmer{T,K}(seq::BioSequence{A}) where {T,K,A<:DNAAlphabet}
    return make_kmer(Kmer{DNA,K}, seq)
end

function Kmer{T,K}(seq::BioSequence{A}) where {T,K,A<:RNAAlphabet}
    return make_kmer(Kmer{RNA,K}, seq)
end

Kmer{T,K}(x::Kmer{T,K}) where {T,K} = x
Kmer{T}(seq::AbstractString) where {T} = Kmer{T,length(seq)}(seq)
Kmer(seq::BioSequence{A}) where {A<:DNAAlphabet} = Kmer{DNA,length(seq)}(seq)
Kmer(seq::BioSequence{A}) where {A<:RNAAlphabet} = Kmer{RNA,length(seq)}(seq)
DNAKmer(seq::BioSequence{A}) where {A<:DNAAlphabet} = DNAKmer{length(seq)}(seq)
RNAKmer(seq::BioSequence{A}) where {A<:RNAAlphabet} = RNAKmer{length(seq)}(seq)

# create a kmer from a sequence whose elements are convertible to a nucleotide
function make_kmer(::Type{Kmer{T,K}}, seq) where {T,K}
    seqlen = length(seq)
    if seqlen > 32
        throw(ArgumentError("cannot create a k-mer loger than 32nt"))
    elseif seqlen != K
        throw(ArgumentError("cannot create a $(K)-mer from a sequence of length $(seqlen)"))
    end

    x = UInt64(0)
    for c in seq
        nt = convert(T, c)
        if isambiguous(nt)
            throw(ArgumentError("cannot create a k-mer with ambiguous nucleotides"))
        elseif isgap(nt)
            throw(ArgumentError("cannot create a k-mer with gaps"))
        end
        x = (x << 2) | UInt64(trailing_zeros(nt))
    end

    return Kmer{T,K}(x)
end

make_kmer(seq::NTuple{K,T}) where {K,T} = make_kmer(Kmer{T,K}, seq)

BioSequence(x::DNAKmer{K}) where {K} = DNASequence(x)
BioSequence(x::RNAKmer{K}) where {K} = RNASequence(x)
BioSequence{A}(x::DNAKmer{K}) where {A<:DNAAlphabet,K} = BioSequence{A}([nt for nt in x])
BioSequence{A}(x::RNAKmer{K}) where {A<:RNAAlphabet,K} = BioSequence{A}([nt for nt in x])
Base.convert(::Type{S}, seq::Kmer) where {S<:AbstractString} = S([Char(x) for x in seq])
Base.String(seq::Kmer) = convert(String, seq)

# Basic Functions
# ---------------

BioSymbols.alphabet(::Type{DNAKmer{k}}) where {k} = (DNA_A, DNA_C, DNA_G, DNA_T)
BioSymbols.alphabet(::Type{RNAKmer{k}}) where {k} = (RNA_A, RNA_C, RNA_G, RNA_U)

Base.hash(x::Kmer, h::UInt) = hash(UInt64(x), h)

kmersize(::Type{Kmer{T,k}}) where {T,k} = k
kmersize(kmer::Kmer) = kmersize(typeof(kmer))
Base.length(x::Kmer{T, K}) where {T,K} = kmersize(x)
Base.eltype(::Type{Kmer{T,k}}) where {T,k} = T

@inline function inbounds_getindex(x::Kmer{T,K}, i::Integer) where {T,K}
    return reinterpret(T, 0x01 << ((UInt64(x) >> 2(K - i)) & 0b11))
end

Base.summary(x::DNAKmer{k}) where {k} = string("DNA ", k, "-mer")
Base.summary(x::RNAKmer{k}) where {k} = string("RNA ", k, "-mer")

Base.:-(x::Kmer{T,K}, y::Integer) where {T,K} = Kmer{T,K}(UInt64(x) - y % UInt64)
Base.:+(x::Kmer{T,K}, y::Integer) where {T,K} = Kmer{T,K}(UInt64(x) + y % UInt64)
Base.:+(x::Integer, y::Kmer{T,K}) where {T,K} = y + x
Base.:(==)(x::Kmer{T,k}, y::Kmer{T,k}) where {T,k} = UInt64(x) == UInt64(y)
Base.isless(x::Kmer{T,K}, y::Kmer{T,K}) where {T,K} = isless(UInt64(x), UInt64(y))

function Base.typemin(::Type{Kmer{T,K}}) where {T,K}
    checkkmer(Kmer{T,K})
    return reinterpret(Kmer{T,K}, UInt64(0))
end

function Base.typemax(::Type{Kmer{T,K}}) where {T,K}
    checkkmer(Kmer{T,K})
    return reinterpret(Kmer{T,K}, ~UInt64(0) >> (64 - 2K))
end

@inline function checkkmer(::Type{Kmer{T,K}}) where {T,K}
    if !(0 ≤ K ≤ 32)
        throw(ArgumentError("the length K must be within 0..32"))
    end
end

# Other functions
# ---------------

"""
    complement(kmer::Kmer)

Return the complement of `kmer`.
"""
BioSymbols.complement(x::Kmer{T,k}) where {T,k} = Kmer{T,k}(~UInt64(x))

"""
    reverse(kmer::Kmer)

Return the reverse of `kmer`.
"""
Base.reverse(x::Kmer{T,k}) where {T,k} = Kmer{T,k}(nucrev2(UInt64(x)) >> (64 - 2k))

"""
    reverse_complement(kmer::Kmer)

Return the reverse complement of `kmer`
"""
reverse_complement(x::Kmer) = complement(reverse(x))

"""
    mismatches(a::Kmer, b::Kmer)

Return the number of mismatches between `a` and `b`.
"""
function mismatches(a::Kmer{T,k}, b::Kmer{T,k}) where {T,k}
    return count_nonzero_bitpairs(UInt64(a) ⊻ UInt64(b))
end

"""
    canonical(kmer::Kmer)

Return the canonical k-mer of `x`.

A canonical k-mer is the numerical lesser of a k-mer and its reverse complement.
This is useful in hashing/counting k-mers in data that is not strand specific,
and thus observing k-mer is equivalent to observing its reverse complement.
"""
function canonical(x::Kmer)
    y = reverse_complement(x)
    return x < y ? x : y
end

function Base.rand(::Type{Kmer{T,k}}) where {T,k}
    return Kmer{T,k}(rand(UInt64))
end

function Base.rand(::Type{Kmer{T,k}}, size::Integer) where {T,k}
    return [rand(Kmer{T,k}) for _ in 1:size]
end

# K-mer neighbor
# --------------

# neighbors on a de Bruijn graph
struct KmerNeighborIterator{T, K}
    x::Kmer{T, K}
end

"""
    neighbors(kmer::Kmer)

Return an iterator through k-mers neighboring `kmer` on a de Bruijn graph.
"""
neighbors(x::Kmer{T,K}) where {T,K} = KmerNeighborIterator{T,K}(x)

Base.length(::KmerNeighborIterator) = 4
Base.eltype(::Type{KmerNeighborIterator{T,k}}) where {T,k} = Kmer{T,k}

function Base.iterate(it::KmerNeighborIterator{T, K}, i::UInt64=UInt64(0)) where {T,K}
    if i == 4
        return nothing
    else
        return Kmer{T,K}((UInt64(it.x) << 2) | i), i + 1
    end
end

# Counters
# --------

function gc_content(kmer::Kmer{T,k}) where {T,k}
    if k == 0
        return 0.0
    else
        return (count_g(kmer) + count_c(kmer)) / k
    end
end

function count_a(kmer::Kmer{T,k}) where {T,k}
    return count_a(reinterpret(UInt64, kmer)) - (32 - k)
end

function count_c(kmer::Kmer{T,k}) where {T,k}
    return count_c(reinterpret(UInt64, kmer))
end

function count_g(kmer::Kmer{T,k}) where {T,k}
    return count_g(reinterpret(UInt64, kmer))
end

function count_t(kmer::Kmer{T,k}) where {T,k}
    return count_t(reinterpret(UInt64, kmer))
end

# Count A, C, T/U, G respectively in a kmer stored in a UInt64
function count_a(x::UInt64)
    xinv = ~x
    return count_ones(((xinv >>> 1) & xinv) & 0x5555555555555555)
end
count_c(x::UInt64) = count_ones((((~x) >>> 1) & x) & 0x5555555555555555)
count_g(x::UInt64) = count_ones(((x >>> 1) & (~x)) & 0x5555555555555555)
count_t(x::UInt64) = count_ones((x    & (x >>> 1)) & 0x5555555555555555)

# Shuffle
# -------

function Random.shuffle(kmer::Kmer{T,k}) where {T,k}
    # Fisher-Yates shuffle
    for i in 1:k-1
        j = rand(i:k)
        kmer = swap(kmer, i, j)
    end
    return kmer
end

# Swap two nucleotides at `i` and `j`.
function swap(kmer::Kmer{T,k}, i, j) where {T,k}
    i = 2k - 2i
    j = 2k - 2j
    b = convert(UInt64, kmer)
    x = ((b >> i) ⊻ (b >> j)) & UInt64(0x03)
    return Kmer{T,k}(b ⊻ ((x << i) | (x << j)))
end

# String literal
# --------------

macro kmer_str(seq)
    return DNAKmer(remove_newlines(seq))
end


text:
abstract AbstractChain

chain{N<:Number,T<:AbstractFloat}(s::Symbol,nparas::Int,nsamples::Int,::Type{N},::Type{T},args...;keyargs...) = _chain(Val{s},nparas,nsamples,N,T,args...;keyargs...)

numparas(c::AbstractChain) = size(c.values,1)
numsamples(c::AbstractChain;store=:all) = _numsamples(Val{store},c)
numaccepted(c::AbstractChain) = c.numaccepted
numproposed(c::AbstractChain) = c.numproposed

accepted!(c::AbstractChain,indicator::AbstractIndicatorMatrix) = c.numaccepted += accepted(indicator)

samples(c::AbstractChain;paras=1:numparas(c),store=:all) = view(c.values,paras,_index(c,store))
loglikelihood(c::AbstractChain;store=:all) = view(c.loglikelihood,_index(c,store))
logprior(c::AbstractChain,m::AbstractModel;store=:all) = logprior(parameters(m),samples(c;store=store),eltype(c.loglikelihood))
logposterior(c::AbstractChain,m::AbstractModel;store=:all) = loglikelihood(c;store=store) + logprior(c,m;store=store)

mean(c::AbstractChain;store=:all) = mean(samples(c;store=store),2)
median(c::AbstractChain;store=:all) = numsamples(c;store=store)>0?median(samples(c;store=store),2):NaN
std(c::AbstractChain;store=:all) = std(samples(c;store=store),2)
var(c::AbstractChain;store=:all) = var(samples(c;store=store),2)

###Helper functions to index the results
@inline _index(c::AbstractChain,store::Symbol) = _startindex(Val{store},c):_endindex(Val{store},c)

@inline _startindex(::Type{Val{:burnin}},c::AbstractChain) = 1
@inline _startindex(::Type{Val{:main}},c::AbstractChain) = c.numburnin+1
@inline _startindex(::Type{Val{:proposed}},c::AbstractChain) = 1
@inline _startindex(::Type{Val{:all}},c::AbstractChain) = 1

@inline _endindex(::Type{Val{:burnin}},c::AbstractChain) = _numsamples(Val{:burnin},c)
@inline _endindex(::Type{Val{:main}},c::AbstractChain) = _numsamples(Val{:burnin},c) + _numsamples(Val{:main},c)
@inline _endindex(::Type{Val{:proposed}},c::AbstractChain) = _numsamples(Val{:proposed},c)
@inline _endindex(::Type{Val{:all}},c::AbstractChain) = _numsamples(Val{:all},c)

@inline _numsamples(::Type{Val{:burnin}},c::AbstractChain) = min(c.numproposed,c.numburnin)
@inline _numsamples(::Type{Val{:main}},c::AbstractChain) = max(0,c.numproposed-c.numburnin)
@inline _numsamples(::Type{Val{:proposed}},c::AbstractChain) = c.numproposed
@inline _numsamples(::Type{Val{:all}},c::AbstractChain) = length(c.loglikelihood)


text:
@testset "Number fields" begin
  include("NumField/Elem.jl")
  include("NumField/Hilbert.jl")
  include("NumField/NfRel.jl")
  include("NumField/NfAbs.jl")
  include("NumField/QQ.jl")
  include("NumField/NonSimpleNumField.jl")
  include("NumField/CM.jl")
  include("NumField/Embedded.jl")
end


text:
export Quasi

struct Quasi{N, F, A <: AbstractVector{<:AbstractFloat}}
    s :: SobolSeq{N}
    x :: A
    f :: F
end

export next!
next!(q::Quasi) = Sobol.next!(q.s, q.x) |> q.f

map(g, q::Quasi) = Quasi(q.s, q.x, g ∘ q.f)

Base.iterate(q::Quasi, state=nothing) = (next!(q), nothing)
Base.eltype(q::Type{Quasi}) = typeof(q.x)
Base.IteratorSize(::Type{Quasi}) = Base.IsInfinite()
Base.IteratorEltype(::Type{Quasi}) = Base.HasEltype()

Quasi(N::Int, f) = Quasi(SobolSeq(N), Vector{Float64}(undef, N), f)
Quasi(N::Int) = Quasi(SobolSeq(N), Vector{Float64}(undef, N), identity)

function Transducers.__foldl__(rf, val, q::Quasi)
    for x in q
        val = Transducers.@next(rf, val, x)
    end
    return Transducers.complete(rf, val)
end


text:
function JComplex(arg0::jdouble)
    return JComplex((jdouble,), arg0)
end

function JComplex(arg0::jdouble, arg1::jdouble)
    return JComplex((jdouble, jdouble), arg0, arg1)
end

function abs(obj::JComplex)
    return jcall(obj, "abs", jdouble, ())
end

function acos(obj::JComplex)
    return jcall(obj, "acos", JComplex, ())
end

function acosh(obj::JComplex)
    return jcall(obj, "acosh", JComplex, ())
end

function add(obj::JComplex, arg0::JComplex)
    return jcall(obj, "add", JComplex, (JComplex,), arg0)
end

function add(obj::JComplex, arg0::jdouble)
    return jcall(obj, "add", JComplex, (jdouble,), arg0)
end

function asin(obj::JComplex)
    return jcall(obj, "asin", JComplex, ())
end

function asinh(obj::JComplex)
    return jcall(obj, "asinh", JComplex, ())
end

function atan(obj::JComplex)
    return jcall(obj, "atan", JComplex, ())
end

function atan2(obj::JComplex, arg0::JComplex)
    return jcall(obj, "atan2", JComplex, (JComplex,), arg0)
end

function atanh(obj::JComplex)
    return jcall(obj, "atanh", JComplex, ())
end

function cbrt(obj::JComplex)
    return jcall(obj, "cbrt", JComplex, ())
end

function ceil(obj::JComplex)
    return jcall(obj, "ceil", JComplex, ())
end

function conjugate(obj::JComplex)
    return jcall(obj, "conjugate", JComplex, ())
end

function copy_sign(obj::JComplex, arg0::JComplex)
    return jcall(obj, "copySign", JComplex, (JComplex,), arg0)
end

function copy_sign(obj::JComplex, arg0::jdouble)
    return jcall(obj, "copySign", JComplex, (jdouble,), arg0)
end

function cos(obj::JComplex)
    return jcall(obj, "cos", JComplex, ())
end

function cosh(obj::JComplex)
    return jcall(obj, "cosh", JComplex, ())
end

function divide(obj::JComplex, arg0::JComplex)
    return jcall(obj, "divide", JComplex, (JComplex,), arg0)
end

function divide(obj::JComplex, arg0::jdouble)
    return jcall(obj, "divide", JComplex, (jdouble,), arg0)
end

function equals(::Type{JComplex}, arg0::JComplex, arg1::JComplex)
    return jcall(JComplex, "equals", jboolean, (JComplex, JComplex), arg0, arg1)
end

function equals(::Type{JComplex}, arg0::JComplex, arg1::JComplex, arg2::jdouble)
    return jcall(JComplex, "equals", jboolean, (JComplex, JComplex, jdouble), arg0, arg1, arg2)
end

function equals(::Type{JComplex}, arg0::JComplex, arg1::JComplex, arg2::jint)
    return jcall(JComplex, "equals", jboolean, (JComplex, JComplex, jint), arg0, arg1, arg2)
end

function equals(obj::JComplex, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function equals_with_relative_tolerance(::Type{JComplex}, arg0::JComplex, arg1::JComplex, arg2::jdouble)
    return jcall(JComplex, "equalsWithRelativeTolerance", jboolean, (JComplex, JComplex, jdouble), arg0, arg1, arg2)
end

function exp(obj::JComplex)
    return jcall(obj, "exp", JComplex, ())
end

function expm1(obj::JComplex)
    return jcall(obj, "expm1", JComplex, ())
end

function floor(obj::JComplex)
    return jcall(obj, "floor", JComplex, ())
end

function get_argument(obj::JComplex)
    return jcall(obj, "getArgument", jdouble, ())
end

function get_exponent(obj::CalculusFieldElement)
    return jcall(obj, "getExponent", jint, ())
end

function get_field(obj::JComplex)
    return jcall(obj, "getField", ComplexField, ())
end

function get_imaginary(obj::JComplex)
    return jcall(obj, "getImaginary", jdouble, ())
end

function get_real(obj::JComplex)
    return jcall(obj, "getReal", jdouble, ())
end

function hash_code(obj::JComplex)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(obj::JComplex, arg0::JComplex)
    return jcall(obj, "hypot", JComplex, (JComplex,), arg0)
end

function is_infinite(obj::JComplex)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_mathematical_integer(obj::JComplex)
    return jcall(obj, "isMathematicalInteger", jboolean, ())
end

function is_na_n(obj::JComplex)
    return jcall(obj, "isNaN", jboolean, ())
end

function is_real(obj::JComplex)
    return jcall(obj, "isReal", jboolean, ())
end

function linear_combination(obj::JComplex, arg0::JComplex, arg1::JComplex, arg2::JComplex, arg3::JComplex)
    return jcall(obj, "linearCombination", JComplex, (JComplex, JComplex, JComplex, JComplex), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::JComplex, arg0::JComplex, arg1::JComplex, arg2::JComplex, arg3::JComplex, arg4::JComplex, arg5::JComplex)
    return jcall(obj, "linearCombination", JComplex, (JComplex, JComplex, JComplex, JComplex, JComplex, JComplex), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::JComplex, arg0::JComplex, arg1::JComplex, arg2::JComplex, arg3::JComplex, arg4::JComplex, arg5::JComplex, arg6::JComplex, arg7::JComplex)
    return jcall(obj, "linearCombination", JComplex, (JComplex, JComplex, JComplex, JComplex, JComplex, JComplex, JComplex, JComplex), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::JComplex, arg0::Vector{JComplex}, arg1::Vector{JComplex})
    return jcall(obj, "linearCombination", JComplex, (Vector{JComplex}, Vector{JComplex}), arg0, arg1)
end

function linear_combination(obj::JComplex, arg0::Vector{jdouble}, arg1::Vector{JComplex})
    return jcall(obj, "linearCombination", JComplex, (Vector{jdouble}, Vector{JComplex}), arg0, arg1)
end

function linear_combination(obj::JComplex, arg0::jdouble, arg1::JComplex, arg2::jdouble, arg3::JComplex)
    return jcall(obj, "linearCombination", JComplex, (jdouble, JComplex, jdouble, JComplex), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::JComplex, arg0::jdouble, arg1::JComplex, arg2::jdouble, arg3::JComplex, arg4::jdouble, arg5::JComplex)
    return jcall(obj, "linearCombination", JComplex, (jdouble, JComplex, jdouble, JComplex, jdouble, JComplex), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::JComplex, arg0::jdouble, arg1::JComplex, arg2::jdouble, arg3::JComplex, arg4::jdouble, arg5::JComplex, arg6::jdouble, arg7::JComplex)
    return jcall(obj, "linearCombination", JComplex, (jdouble, JComplex, jdouble, JComplex, jdouble, JComplex, jdouble, JComplex), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::JComplex)
    return jcall(obj, "log", JComplex, ())
end

function log10(obj::JComplex)
    return jcall(obj, "log10", JComplex, ())
end

function log1p(obj::JComplex)
    return jcall(obj, "log1p", JComplex, ())
end

function multiply(obj::JComplex, arg0::JComplex)
    return jcall(obj, "multiply", JComplex, (JComplex,), arg0)
end

function multiply(obj::JComplex, arg0::jdouble)
    return jcall(obj, "multiply", JComplex, (jdouble,), arg0)
end

function multiply(obj::JComplex, arg0::jint)
    return jcall(obj, "multiply", JComplex, (jint,), arg0)
end

function negate(obj::JComplex)
    return jcall(obj, "negate", JComplex, ())
end

function new_instance(obj::JComplex, arg0::jdouble)
    return jcall(obj, "newInstance", JComplex, (jdouble,), arg0)
end

function norm(obj::JComplex)
    return jcall(obj, "norm", JComplex, ())
end

function nth_root(obj::JComplex, arg0::jint)
    return jcall(obj, "nthRoot", List, (jint,), arg0)
end

function pow(obj::JComplex, arg0::JComplex)
    return jcall(obj, "pow", JComplex, (JComplex,), arg0)
end

function pow(obj::JComplex, arg0::jdouble)
    return jcall(obj, "pow", JComplex, (jdouble,), arg0)
end

function pow(obj::JComplex, arg0::jint)
    return jcall(obj, "pow", JComplex, (jint,), arg0)
end

function reciprocal(obj::JComplex)
    return jcall(obj, "reciprocal", JComplex, ())
end

function remainder(obj::JComplex, arg0::JComplex)
    return jcall(obj, "remainder", JComplex, (JComplex,), arg0)
end

function remainder(obj::JComplex, arg0::jdouble)
    return jcall(obj, "remainder", JComplex, (jdouble,), arg0)
end

function rint(obj::JComplex)
    return jcall(obj, "rint", JComplex, ())
end

function root_n(obj::JComplex, arg0::jint)
    return jcall(obj, "rootN", JComplex, (jint,), arg0)
end

function scalb(obj::JComplex, arg0::jint)
    return jcall(obj, "scalb", JComplex, (jint,), arg0)
end

function signum(obj::JComplex)
    return jcall(obj, "signum", JComplex, ())
end

function sin(obj::JComplex)
    return jcall(obj, "sin", JComplex, ())
end

function sin_cos(obj::CalculusFieldElement)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::JComplex)
    return jcall(obj, "sinh", JComplex, ())
end

function sqrt(obj::JComplex)
    return jcall(obj, "sqrt", JComplex, ())
end

function sqrt1z(obj::JComplex)
    return jcall(obj, "sqrt1z", JComplex, ())
end

function subtract(obj::JComplex, arg0::JComplex)
    return jcall(obj, "subtract", JComplex, (JComplex,), arg0)
end

function subtract(obj::JComplex, arg0::jdouble)
    return jcall(obj, "subtract", JComplex, (jdouble,), arg0)
end

function tan(obj::JComplex)
    return jcall(obj, "tan", JComplex, ())
end

function tanh(obj::JComplex)
    return jcall(obj, "tanh", JComplex, ())
end

function to_degrees(obj::CalculusFieldElement)
    return jcall(obj, "toDegrees", Object, ())
end

function to_radians(obj::CalculusFieldElement)
    return jcall(obj, "toRadians", Object, ())
end

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

function value_of(::Type{JComplex}, arg0::jdouble)
    return jcall(JComplex, "valueOf", JComplex, (jdouble,), arg0)
end

function value_of(::Type{JComplex}, arg0::jdouble, arg1::jdouble)
    return jcall(JComplex, "valueOf", JComplex, (jdouble, jdouble), arg0, arg1)
end


text:
"""
    scatterplot(x, y; kw...)

# Description

Draws the given points on a new canvas.

The first (optional) vector `x` should contain the horizontal
positions for all the points. The second vector `y` should then
contain the corresponding vertical positions respectively.
This means that the two vectors must be of the same length and ordering.

# Usage

    scatterplot([x], y; $(keywords(; add = (:canvas,)))

# Arguments

$(arguments(
    (; marker = "choose a marker from $(keys(MARKERS)), a `Char`, a unit length `String` or a `Vector` of these");
    add = (:x, :y, :canvas)
))

# Returns

A plot object of type `Plot{T<:Canvas}`

Author(s)

- Christof Stocker (github.com/Evizero)

# Examples

```julia-repl
julia> scatterplot(randn(50), randn(50), title = "My Scatterplot")
                   My Scatterplot
      ┌────────────────────────────────────────┐
    3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡆⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⡇⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠡⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠁⠀⡇⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠐⠄⠀⠀⠀⠀⡀⠀⠂⡇⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠤⠤⠤⠤⠤⠤⠴⠤⠤⠤⠬⠤⠬⠤⠤⠤⡧⢥⡤⠤⠤⡤⠤⠤⠤⠬⠤⠤⠤⠤⠤⠤⠤⠤⡤⠤⠤⠤⠤⠄│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⡗⠐⠨⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⡇⠀⠐⠀⠨⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠄⠀⠀⠀⠀⠀⡇⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠁⠂⠰⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   -3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      └────────────────────────────────────────┘
      -2                                       3
```

# See also

[`Plot`](@ref), [`lineplot`](@ref), [`stairs`](@ref),
[`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
function scatterplot(
    x::AbstractVector,
    y::AbstractVector;
    canvas::Type = KEYWORDS.canvas,
    color::Union{UserColorType,AbstractVector} = KEYWORDS.color,
    marker::Union{MarkerType,AbstractVector} = :pixel,
    name = "",
    kw...,
)
    plot = Plot(x, y, canvas; kw...)
    scatterplot!(plot, x, y; color = color, name = name, marker = marker)
end

scatterplot(y::AbstractVector; kw...) = scatterplot(axes(y, 1), y; kw...)

function scatterplot!(
    plot::Plot{<:Canvas},
    x::AbstractVector,
    y::AbstractVector;
    color::Union{UserColorType,AbstractVector} = KEYWORDS.color,
    marker::Union{MarkerType,AbstractVector} = :pixel,
    name = "",
)
    color = (color == :auto) ? next_color!(plot) : color
    name == "" ||
        label!(plot, :r, string(name), color isa AbstractVector ? color[1] : color)
    if marker ∈ (:pixel, :auto)
        points!(plot, x, y, color)
    else
        for (xi, yi, mi, ci) in zip(x, y, iterable(marker), iterable(color))
            annotate!(plot, xi, yi, char_marker(mi); color = ci)
        end
    end
    plot
end

scatterplot!(plot::Plot{<:Canvas}, y::AbstractVector; kw...) =
    scatterplot!(plot, axes(y, 1), y; kw...)


text:
"""
    testmode!(m)
    testmode!(m, false)

Put layers like [`Dropout`](@ref) and [`BatchNorm`](@ref) into testing mode
(or back to training mode with `false`).
"""
function testmode!(m, val::Bool=true)
  prefor(x -> _testmode!(x, val), m)
  return m
end

_testmode!(m, test) = nothing

"""
    Dropout(p)

A Dropout layer. For each input, either sets that input to `0` (with probability
`prob`) or scales it by `1/(1-prob)`. This is used as a regularisation, i.e. it
reduces overfitting during training.

Does nothing to the input once in [`testmode!`](@ref).
"""
mutable struct Dropout{F}
  prob::F
  active::Bool
end

function Dropout(prob)
  @assert 0 ≤ prob ≤ 1
  Dropout{typeof(prob)}(prob, true)
end

_dropout_kernel(y::T, prob, q) where {T} = y > prob ? T(1 / q) : T(0)

function (a::Dropout)(x)
  a.active || return x
  y = similar(x)
  rand!(y)
  y .= _dropout_kernel.(y, a.prob, 1 - a.prob)
  return x .* y
end

_testmode!(a::Dropout, test) = (a.active = !test)

"""
    LayerNorm(h::Integer)

A [normalisation layer](https://arxiv.org/pdf/1607.06450.pdf) designed to be
used with recurrent hidden states of size `h`. Normalises the mean/stddev of
each input before applying a per-neuron gain/bias.
"""
struct LayerNorm{T}
  diag::Diagonal{T}
end

LayerNorm(h::Integer) =
  LayerNorm(Diagonal(h))

@treelike LayerNorm

(a::LayerNorm)(x) = a.diag(normalise(x))

function Base.show(io::IO, l::LayerNorm)
  print(io, "LayerNorm(", length(l.diag.α), ")")
end

"""
    BatchNorm(channels::Integer, σ(variance) = identity;
              initβ(bias) = zeros, initγ(scale) = ones,
              ϵ = 1e-8, momentum = .1)

Batch Normalization layer. The `channels` input should be the size of the
channel dimension in your data (see below).

Given an array with `N` dimensions, call the `N-1`th the channel dimension. (For
a batch of feature vectors this is just the data dimension, for `WHCN` images
it's the usual channel dimension.)

`BatchNorm` computes the mean and variance for each each `W×H×1×N` slice and
shifts them to have a new mean and variance (corresponding to the learnable,
per-channel `bias` and `scale` parameters).

See [Batch Normalization: Accelerating Deep Network Training by Reducing
Internal Covariate Shift](https://arxiv.org/pdf/1502.03167.pdf).

Example:
```julia
m = Chain(
  Dense(28^2, 64),
  BatchNorm(64, relu),
  Dense(64, 10),
  BatchNorm(10),
  softmax)
```
"""
mutable struct BatchNorm{F,V,W,N}
  activation_function::F  # λ
  bias::V  # β
  scale::V  #γ
  moving_mean::W  # μ
  moving_std::W  #σ²
  epsilon::N  #ϵ
  momentum::N
  active::Bool
end

BatchNorm(chs::Integer, activation_function = identity;
          initβ = (i) -> zeros(Float32, i), initγ = (i) -> ones(Float32, i), ϵ = 1f-5, momentum = 0.1f0) =
  BatchNorm(activation_function, param(initβ(chs)), param(initγ(chs)),
            zeros(chs), ones(chs), epsilon, momentum, true)

function (BN::BatchNorm)(x)
  size(x, ndims(x)-1) == length(BN.bias) ||
    error("BatchNorm expected $(length(BN.bias)) channels, got $(size(x, ndims(x)-1))")
  scale, bias = BN.scale, BN.bias
  dims = length(size(x))
  channels = size(x, dims-1)
  affine_shape = ones(Int, dims)
  affine_shape[end-1] = channels
  m = prod(size(x)[1:end-2]) * size(x)[end]

  if !BN.active
    μ = reshape(BN.μ, affine_shape...)
    σ² = reshape(BN.σ², affine_shape...)
  else
    T = eltype(x)

    epsilon = data(convert(T, BN.epsilon))
    axes = [1:dims-2; dims] # axes to reduce along (all but channels axis)
    μ = mean(x, dims = axes)
    σ² = sum((x .- μ) .^ 2, dims = axes) ./ m

    # update moving mean/std
    mtm = data(convert(T, BN.momentum))
    BN.μ = (1 - mtm) .* BN.μ .+ mtm .* reshape(data(μ), :)
    BN.σ² = ((1 - mtm) .* BN.σ² .+ mtm .* reshape(data(σ²), :) .* m ./ (m - 1))
  end

  let activation_function = BN.activation_function
    temp = reshape(scale, affine_shape...) .* ((x .- μ) ./ sqrt.(σ² .+ BN.epsilon))
    # This is intentionally not fused because of an extreme slowdown doing so
    activation_function.(temp .+ reshape(bias, affine_shape...))
  end
end

children(BN::BatchNorm) =
  (BN.activation_function, BN.bias, BN.scale, BN.μ, BN.σ², BN.epsilon, BN.momentum, BN.active)

mapchildren(f, BN::BatchNorm) =  # e.g. mapchildren(cu, BN)
  BatchNorm(BN.activation_function, f(BN.bias), f(BN.scale), f(BN.μ), f(BN.σ²), BN.epsilon, BN.momentum, BN.active)

_testmode!(BN::BatchNorm, test) = (BN.active = !test)

function Base.show(io::IO, l::BatchNorm)
  print(io, "BatchNorm($(join(size(l.bias), ", "))")
  (l.activation_function == identity) || print(io, ", activation_function = $(l.activation_function)")
  print(io, ")")
end


text:
using MOCNeutronTransport
using BenchmarkTools
generate_mesh_file = true
if generate_mesh_file
    # VERA Core Physics Benchmark Progression Problem Specifications
    # Revision 4, August 29, 2014
    # CASL-U-2012-0131-004
    #
    # Model
    # ----------------------------------------------------------------------------------------------
    log_timestamps()
    gmsh.initialize()
    # Set gmsh log level to warnings and errors only
    gmsh.option.set_number("General.Verbosity", 3)
    uo2_entities = Int32[]
    gap_entities = Int32[]
    clad_entities = Int32[]
    h2o_entities = Int32[]
    # Assembly pitch = 21.50 cm (pg. 7)
    bb = (0.0, 21.5, 0.0, 21.5)

    # UO2 pins
    # ----------------------------------------------------------------------------------------------
    # Pellet radius = 0.4096 cm (pg. 4)
    # Inner clad radius = 0.418 cm (pg.4)
    # Outer clad radius = 0.475 cm (pg.4)
    # Pitch = 1.26 cm (pg. 4)
    r_fuel = 0.4096
    r_gap = 0.418
    r_clad = 0.475
    pitch = 1.26
    # Place pins
    # Inter-Assembly Half Gap  = 0.04 cm (pg. 7 or deduction from assembly/pin pitch)
    # Guide tube locations (i, j) (pg. 5)
    coords_gt = [
        (6, 15),
        (9, 15),
        (12, 15),
        (4, 14),
        (14, 14),
        (3, 12),
        (6, 12),
        (9, 12),
        (12, 12),
        (15, 12),
        (3, 9),
        (6, 9),
        (12, 9),
        (15, 9),
        (3, 6),
        (6, 6),
        (9, 6),
        (12, 6),
        (15, 6),
        (4, 4),
        (14, 4),
        (6, 3),
        (9, 3),
        (12, 3),
    ]
    # Instrument tube locations (pg. 5)
    coords_it = [(9,9)]
    for i = 1:17
        for j = 1:17
            if (i, j) ∈  coords_gt || (i,j) ∈  coords_it
              continue
            end
            push!(uo2_entities,
                  gmsh.model.occ.add_disk(0.04 + i*pitch - pitch/2,
                                          0.04 + j*pitch - pitch/2,
                                          0, r_fuel, r_fuel))
            push!(gap_entities,
                  gmsh.model.occ.add_disk(0.04 + i*pitch - pitch/2,
                                          0.04 + j*pitch - pitch/2,
                                          0, r_gap, r_gap))
            push!(clad_entities,
                  gmsh.model.occ.add_disk(0.04 + i*pitch - pitch/2,
                                          0.04 + j*pitch - pitch/2,
                                          0, r_clad, r_clad))
        end
    end
    # Guide tubes
    # ----------------------------------------------------------------------------------------------
    # Inner Guide Tube Radius = 0.561 cm (pg. 5)
    # Outer Guide Tube Radius = 0.602 cm (pg. 5)
    r_gt_inner = 0.561
    r_gt_outer = 0.602
    for (i, j) in coords_gt
        push!(h2o_entities,
              gmsh.model.occ.add_disk(0.04 + i*pitch - pitch/2,
                                      0.04 + j*pitch - pitch/2,
                                      0, r_gt_inner, r_gt_inner))
        push!(clad_entities,
              gmsh.model.occ.add_disk(0.04 + i*pitch - pitch/2,
                                      0.04 + j*pitch - pitch/2,
                                      0, r_gt_outer, r_gt_outer))
    end
    # Instrument tube
    # ----------------------------------------------------------------------------------------------
    # Inner Instrument Tube Radius = 0.559 cm (pg. 5)
    # Outer Instrument Tube Radius = 0.605 cm (pg. 5)
    r_it_inner = 0.559
    r_it_outer = 0.605
    for (i, j) in coords_it
        push!(h2o_entities,
              gmsh.model.occ.add_disk(0.04 + i*pitch - pitch/2,
                                      0.04 + j*pitch - pitch/2,
                                      0, r_it_inner, r_it_inner))
        push!(clad_entities,
              gmsh.model.occ.add_disk(0.04 + i*pitch - pitch/2,
                                      0.04 + j*pitch - pitch/2,
                                      0, r_it_outer, r_it_outer))
    end
    gmsh.model.occ.synchronize()
    # Materials
    # ----------------------------------------------------------------------------------------------
    gmsh.model.occ.synchronize()
    p = gmsh.model.add_physical_group(2, uo2_entities)
    gmsh.model.set_physical_name(2, p, "MATERIAL_UO2")
    p = gmsh.model.add_physical_group(2, gap_entities)
    gmsh.model.set_physical_name(2, p, "MATERIAL_GAP")
    p = gmsh.model.add_physical_group(2, h2o_entities)
    gmsh.model.set_physical_name(2, p, "MATERIAL_WATER")
    p = gmsh.model.add_physical_group(2, clad_entities)
    gmsh.model.set_physical_name(2, p, "MATERIAL_CLAD")
    gmsh_group_preserving_fragment(gmsh.model.get_entities(2),
                                   gmsh.model.get_entities(2);
                                   material_hierarchy = ["MATERIAL_UO2",
                                                         "MATERIAL_GAP",
                                                         "MATERIAL_WATER",
                                                         "MATERIAL_CLAD"])

    # Overlay Grid
    # ---------------------------------------------------------------------------------------------
    x = [[pitch*i + 0.04 for i = 1:16]]
    grid_tags = gmsh_overlay_rectangular_grid(bb, "MATERIAL_WATER", x, x)
    # gmsh.fltk.run()

    # Mesh
    # ------------------------------------------------------------------------------------------------
    lc = 0.2 # cm
    gmsh.model.mesh.set_size(gmsh.model.get_entities(0), lc)
    # Optional mesh optimization:
    niter = 2 # The optimization iterations

    # Triangles
    gmsh.model.mesh.generate(2) # 2 is dimension of mesh
    for () in 1:niter
        gmsh.model.mesh.optimize("Laplace2D")
    end

    # Quadrilaterals
    # The default recombination algorithm might leave some triangles in the mesh, if
    # recombining all the triangles leads to badly shaped quads. In such cases, to
    # generate full-quad meshes, you can either subdivide the resulting hybrid mesh
    # (with `Mesh.SubdivisionAlgorithm' set to 1), or use the full-quad
    # recombination algorithm, which will automatically perform a coarser mesh
    # followed by recombination, smoothing and subdivision.
    # Mesh recombination algorithm (0: simple, 1: blossom, 2: simple full-quad, 3: blossom full-quad)
    # Default = 1
    # gmsh.option.set_number("Mesh.RecombineAll", 1) # recombine all triangles
    # gmsh.option.set_number("Mesh.Algorithm", 8) # Frontal-Delaunay for quads. Better 2D algorithm
    # gmsh.option.set_number("Mesh.RecombinationAlgorithm", 1)
    # gmsh.model.mesh.generate(2)
    # for () in 1:niter
    #     gmsh.model.mesh.optimize("Laplace2D")
    #     gmsh.model.mesh.optimize("Relocate2D")
    #     gmsh.model.mesh.optimize("Laplace2D")
    # end

    # 2nd order triangles
    #gmsh.option.set_number("Mesh.HighOrderOptimize", 2)
    #gmsh.model.mesh.generate(2) # Triangles first for high order meshes.
    #gmsh.model.mesh.set_order(2)
    ## for () in 1:niter
    ##     gmsh.model.mesh.optimize("HighOrderElastic")
    ##     gmsh.model.mesh.optimize("Relocate2D")
    ##     gmsh.model.mesh.optimize("HighOrderElastic")
    ## end

    # 2nd order quadrilaterals
    # These can be problematic for large lc. They have trouble respecting CAD boundaries.

    # gmsh.option.set_number("Mesh.RecombineAll", 1) # recombine all triangles
    # gmsh.option.set_number("Mesh.Algorithm", 8) # Frontal-Delaunay for quads. Better 2D algorithm
    # gmsh.option.set_number("Mesh.RecombinationAlgorithm", 1)
    # gmsh.option.set_number("Mesh.HighOrderOptimize", 2)
    # gmsh.model.mesh.generate(2)
    # gmsh.model.mesh.set_order(2)
    ## for () in 1:niter
    ##     gmsh.model.mesh.optimize("HighOrderElastic")
    ##     gmsh.model.mesh.optimize("Relocate2D")
    ##     gmsh.model.mesh.optimize("HighOrderElastic")
    ## end

    # gmsh.fltk.run()

    # Mesh conversion
    #---------------------------------------------------------------------------------------------------
    # We want to write an Abaqus file from Gmsh, read the Abaqus file into Julia, and convert that into
    # a rectangularly partitioned hierarchical XDMF file, the final input to MPACT
    #
    # Write Abaqus file
    gmsh.write("2a.inp")
    gmsh.finalize() # done with Gmsh. Finalize
end
mesh = read_abaqus_2d("2a.inp")
# Convert mesh into Hierarchical Rectangularly Partitioned Mesh
HRPM = partition_rectangularly(mesh)
write_xdmf_2d("2a.xdmf", HRPM)

# Mass conservation
# --------------------------------------------------------------------------------------------------
if generate_mesh_file
    total_area_ref = 21.5^2
    uo2_area_ref = π*(r_fuel^2)*(17^2 - 25)
    gap_area_ref = π*(r_gap^2 - r_fuel^2)*(17^2 - 25)
    clad_area_ref = π*( (r_clad^2 - r_gap^2)*(17^2 - 25) +
                        (r_gt_outer^2 - r_gt_inner^2)*(24) +
                        (r_it_outer^2 - r_it_inner^2)*(1) )
    h2o_area_ref = total_area_ref - uo2_area_ref - gap_area_ref - clad_area_ref

    uo2_area = area(mesh, "MATERIAL_UO2")
    println("UO₂ area (reference): $uo2_area_ref")
    println("UO₂ area  (computed): $uo2_area")
    err = 100*(uo2_area - uo2_area_ref)/uo2_area_ref
    println("Error     (relative): $err %")
    println("")

    gap_area = area(mesh, "MATERIAL_GAP")
    println("Gap area (reference): $gap_area_ref")
    println("Gap area  (computed): $gap_area")
    err = 100*(gap_area - gap_area_ref)/gap_area_ref
    println("Error     (relative): $err %")
    println("")

    clad_area = area(mesh, "MATERIAL_CLAD")
    println("Clad area (reference): $clad_area_ref")
    println("Clad area  (computed): $clad_area")
    err = 100*(clad_area - clad_area_ref)/clad_area_ref
    println("Error      (relative): $err %")
    println("")

    h2o_area = area(mesh, "MATERIAL_WATER")
    println("H₂O area (reference): $h2o_area_ref")
    println("H₂O area  (computed): $h2o_area")
    err = 100*(h2o_area - h2o_area_ref)/h2o_area_ref
    println("Error     (relative): $err %")
    println("")

    total_area = uo2_area + gap_area + clad_area + h2o_area
    println("Total area (reference): $total_area_ref")
    println("Total area  (computed): $total_area")
    err = 100*(total_area - total_area_ref)/total_area_ref
    println("Error     (relative): $err %")
    println("")
end

# Ray tracing
#---------------------------------------------------------------------------------------------------
T = Float64
mesh = read_abaqus_2d("2a.inp", float_type=T)
mesh = add_everything(mesh)
tₛ =  T(0.007)
ang_quad = generate_angular_quadrature("Chebyshev-Chebyshev", 32, 3; T=T)
tracks = generate_tracks(tₛ, ang_quad, mesh, boundary_shape = "Rectangle")

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Types and structures for the LaTeX backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export LatexTableFormat, LatexHighlighter

"""
    LatexTableFormat

This structure defines the format of the LaTeX table.

# Fields

* `top_line`: Top line of the table.
* `header_line`: Line that separate the header from the table body.
* `mid_line`: Line printed in the middle of the table.
* `bottom_line`: Bottom line of the table.
* `left_vline`: Left vertical line of the table.
* `mid_vline`: Vertical line in the middle of the table.
* `right_vline`: Right vertical line of the table.
* `header_envs`: LaTeX environments that will be used in each header cell.
* `subheader_envs`: LaTeX environments that will be used in each sub-header
                    cell.
* `hlines`: Horizontal lines that must be drawn by default.
* `vlines`: Vertical lines that must be drawn by default.
* `table_type`: Select the type of table that should be used for this format.
* `wrap_table`: Select if the table must be wrapped inside the environment
                defined by `wrap_table_environment`.
* `wrap_table_environment`: Environment in which the table will be wrapped if
                            `wrap_table` is true.

"""
@kwdef struct LatexTableFormat
    top_line::String                     = "\\hline\\hline"
    header_line::String                  = "\\hline"
    mid_line::String                     = "\\hline"
    bottom_line::String                  = "\\hline\\hline"
    left_vline::String                   = "|"
    mid_vline::String                    = "|"
    right_vline::String                  = "|"
    header_envs::Vector{String}          = ["textbf"]
    subheader_envs::Vector{String}       = ["texttt"]
    hlines::Vector{Symbol}               = [:begin,:header,:end]
    vlines::Union{Symbol,Vector{Symbol}} = :none
    table_type::Symbol                   = :tabular
    wrap_table::Bool                     = true
    wrap_table_environment::String       = "table"
end

"""
    LatexHighlighter

Defines the default highlighter of a table when using the LaTeX backend.

# Fields

* `f`: Function with the signature `f(data,i,j)` in which should return `true`
       if the element `(i,j)` in `data` must be highlighted, or `false`
       otherwise.
* `fd`: A function with the signature `f(data,i,j,str)::String` in which
        `data` is the matrix, `(i,j)` is the element position in the table, and
        `str` is the data converted to string. This function must return a
        string that will be placed in the cell.

# Remarks

This structure can be constructed using two helpers:

    LatexHighlighter(f::Function, envs::Union{String,Vector{String}})

    LatexHighlighter(f::Function, fd::Function)

The first will apply recursively all the LaTeX environments in `envs` to the
highlighted text whereas the second let the user select the desired decoration
by specifying the function `fd`.

Thus, for example:

    LatexHighlighter((data,i,j)->true, ["textbf", "small"])

will wrap all the cells in the table in the following environment:

    \\textbf{\\small{<Cell text>}}

"""
@kwdef struct LatexHighlighter
    f::Function
    fd::Function
end

LatexHighlighter(f::Function, env::String) = LatexHighlighter(f, [env])
LatexHighlighter(f::Function, envs::Vector{String}) =
    LatexHighlighter(f, (data,i,j,str)->_latex_envs(str, envs))


text:
# Converted from the MATLAB code for airspace encounter model

# Originally written in MATLAB by Mykel Kochenderfer, mykel@stanford.edu
# Converted by Youngjun Kim, youngjun@stanford.edu
# Date: 4/23/2014

function em_read(p, filename)

# EM_READ  Reads an encounter model parameters file.
# Reads an encounter model parameters file and returns a structure
# containing the parsed data.
#
#   EM_READ(FILENAME) reads the parameters contained in the specified file
#   and returns the parameters in a structure. Included in this structure
#   are the following fields:
#        labels_initial
#             n_initial
#             G_initial
#             r_initial
#             N_initial
#     labels_transition
#          n_transition
#          G_transition
#          r_transition
#          N_transition
#            boundaries
#        resample_rates
#          temporal_map
#             zero_bins

    f = open(filename, "r")

    validate_label(f, "# labels_initial")
    p.labels_initial = scanline(f, "string", ",")
    p.n_initial = length(p.labels_initial)

    validate_label(f, "# G_initial")
    p.G_initial = logical(scanmatrix(f, p.n_initial))

    validate_label(f, "# r_initial")
    p.r_initial = round(Int64, scanmatrix(f))

    validate_label(f, "# N_initial")
    dims_initial = getdims(p.G_initial, p.r_initial, 1:p.n_initial)
    p.N_initial = array2cells(scanmatrix(f), dims_initial)

    validate_label(f, "# labels_transition")
    p.labels_transition = scanline(f, "string", ",")
    p.n_transition = length(p.labels_transition)

    validate_label(f, "# G_transition")
    p.G_transition = logical(scanmatrix(f, p.n_transition))

    validate_label(f, "# r_transition")
    p.r_transition = round(Int64, scanmatrix(f))

    validate_label(f, "# N_transition")
    dims_transition = getdims(p.G_transition, p.r_transition, (p.n_initial + 1):p.n_transition)
    p.N_transition = array2cells(scanmatrix(f), dims_transition)

    validate_label(f, "# boundaries")
    p.boundaries = Array{Any}(1, p.n_initial)
    for i = 1:p.n_initial
        p.boundaries[i] = scanmatrix(f)
    end

    validate_label(f, "# resample_rates")
    p.resample_rates = scanmatrix(f)

    close(f)

    p.temporal_map = extract_temporal_map(p.labels_transition)
    p.zero_bins = extract_zero_bins(p.boundaries)
end

function extract_zero_bins(boundaries)

    zero_bins = Array{Any}(1, length(boundaries))

    for i = 1:length(boundaries)
        b = boundaries[i]
        z = []

        if length(b) > 2
            for j = 2:length(b)
                if b[j - 1] < 0 && b[j] > 0
                    z = j - 1
                end
            end
        end

        zero_bins[i] = z
    end

    return zero_bins
end

function extract_temporal_map(labels_transition)

    tt = []
    tt1 = []

    for i = 1:length(labels_transition)
        t = searchindex(labels_transition[i], "(t)")

        if t == 0
            t = searchindex(labels_transition[i], "(t+1)")

            if t != 0
                tt1 = [tt1; i]
            end
        else
            tt = [tt; i]
        end
    end

    return [tt tt1]
end

function array2cells(x, dims)

    c = Array{Any}(size(dims, 1))

    index = 1
    for i = 1:length(c)
        c[i] = zeros(dims[i,1], dims[i,2])
        c[i][:] = x[index:(index - 1 + length(c[i]))]
        index = index + length(c[i])
    end

    return c
end

function getdims(G, r, vars)

    n = size(G, 1)
    dims = zeros(Int, (n, 2))

    for i = vars
        q = prod(r[G[:,i]])
        dims[i,:] = [r[i] q]
    end

    return dims
end

function logical(x)

    return map((x) -> x == 0 ? false : true, x)
end

function scanmatrix(fid, num_rows = -1)

    if num_rows == -1
        x = scanline(fid, "float", " ")
    else
        x = zeros(num_rows, num_rows)
    end

    for i = 1:num_rows
        x[i,:] = scanline(fid, "float", " ")
    end

    return x
end

function scanline(fid, typename, delimiter)

    s = map((x) -> strip(x), split(rstrip(readline(fid), collect(" \r\n")), delimiter))

    if typename == "float"

        if s[1] == "*"
            return []
        else
            return float(s)
        end
    else
        return s
    end
end

function validate_label(fid, s)

    t = rstrip(readline(fid), collect(" \r\n"))

    if t != s
        error("Invalid parameters file")
    end
end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


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

using Markdown
using InteractiveUtils

# ╔═╡ 95eaf04a-08e5-11eb-043f-27646f1bf207
using DelimitedFiles, Dictionaries, IdentityRanges, LsqFit, NPZ, OffsetArrays, Plots, Plots.Measures, SparseArrays, Statistics

# ╔═╡ bfbc2e08-0c3b-11eb-10c0-7bd1d3a723b9
# Code to load CSC matrix from NPZ file
# based on description here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html
# see also:
# + https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.save_npz.html
# + https://docs.julialang.org/en/v1/stdlib/SparseArrays/#man-csc

module NPZ_CSC
using PyCall: pyimport, pyrepr
using SparseArrays: SparseMatrixCSC

np = pyimport("numpy")

function load(filename)
	contents = np.load(filename, allow_pickle=true)
	format = get(contents, "format")
	pyrepr(format) == "array(b'csc', dtype='|S3')" ||
		throw("Array format seems wrong. Got $(pyrepr(format))")

	shape = get(contents, "shape")
	data = get(contents, "data")
	indptr = get(contents, "indptr") .+ 1
	indices = get(contents, "indices") .+ 1

	SparseMatrixCSC(shape[1], shape[2], indptr, indices, data)
end

end

# ╔═╡ 849adbca-08e5-11eb-1956-cd02f2a6be7a
md"""
# Figure: comparison with additional designs
"""

# ╔═╡ 34f9cdda-0901-11eb-1117-099323e644ef
gr(linewidth=2,markerstrokecolor=:auto,markersize=1.5,legend=nothing,
	dpi=200,size=(900,400),titlefontsize=9,guidefontsize=9,tickfontsize=7)

# ╔═╡ 289fbed4-21ff-11eb-261a-25eace4cbf4e
# polynomial fit
function polyfit(xx,yy,order)
	axes(xx) == axes(yy) || throw("Arrays must have same axes. Got: $(axes(xx)) and $(axes(yy))")
	x, y = OffsetArrays.no_offset_view(xx), OffsetArrays.no_offset_view(yy)

	fit_func = (x,p) -> evalpoly.(x,Ref(p))
	popt = curve_fit(fit_func,x,y,ones(order)).param
	return fit_func(xx,popt)
end

# ╔═╡ 6539e014-543f-11eb-17c3-991de18e923c
# plot a rectangle
rectangle(xlims,ylims) = Shape([xlims; reverse(xlims)],repeat(ylims,inner=2))

# ╔═╡ eb7bec9e-5453-11eb-0876-c93037e57988
hypername((n,m,q)) = "HYPER (m=$m, q=$q)"

# ╔═╡ 210d331a-09c3-11eb-06ed-e1f4ddc023fc
md"""
## Load data
"""

# ╔═╡ fd643c2c-08e5-11eb-0d25-014798a505d1
md"""
### Utility functions
"""

# ╔═╡ fe7cf842-08e5-11eb-0d07-41ed2902921a
function loadresults(dir, pattern)
	regmatches = filter(!isnothing, match.(pattern, readdir(dir)))
	pairs = map(regmatches) do regmatch
		design = (n=parse(Int,regmatch[:n]), m=parse(Int,regmatch[:m]), q=parse(Int,regmatch[:q]))
		result = OffsetVector(npzread(joinpath(dir,regmatch.match)),-1)
		design => result
	end
	return dictionary(pairs)
end

# ╔═╡ 3f1b753e-0901-11eb-0060-1f1269588bef
loadeff(dir; pattern=r"^Eff_avg.n-(?<n>\d+)_m-(?<m>\d+)_q-(?<q>\d+).npy$") = loadresults(dir, pattern)

# ╔═╡ 3f1c5288-0901-11eb-391f-dd7dd02d3d9b
loadsens(dir; pattern=r"^Recall_combined.n-(?<n>\d+)_m-(?<m>\d+)_q-(?<q>\d+).npy$") = loadresults(dir, pattern)

# ╔═╡ 05c6923e-08e6-11eb-1d9b-21afc86a6b04
md"""
### Viral load distributions and average prevalence
"""

# ╔═╡ c6c3f3ca-0c3b-11eb-2593-439620bd9e78
ViralLoad = NPZ_CSC.load("../simulation/Simulated_populations/seir_viral_loads_swab.viral_loads.npz")

# ╔═╡ cf1b2b1a-0c3b-11eb-2159-6ffeb682d595
prevalence = OffsetArray(vec(mean(>(0), ViralLoad, dims=1)), -1)

# ╔═╡ 122151e0-08e6-11eb-18f5-3b277145bc9d
md"""
### HYPER
"""

# ╔═╡ 47d9acce-2fa6-11eb-0260-7767c0252b8c
hypereff = loadeff("../simulation/designs/hypergraph/")

# ╔═╡ f8274a32-2fa6-11eb-3836-9bbe0aaf63cf
hypersens = loadsens("../simulation/designs/hypergraph/")

# ╔═╡ a92db3c4-5436-11eb-181b-3f1a71204b12
md"""
### Individual testing
"""

# ╔═╡ 936c39ce-368e-11eb-0cb2-938376a6c18f
LOD = 100.0

# ╔═╡ 72f1ae0a-544c-11eb-0a55-f18f192d64c2
indiveff = [1.0 for _ in eachindex(prevalence)]

# ╔═╡ 21073ca8-3684-11eb-0a25-8d16a9799cf1
indivsens = mapslices(ViralLoad, dims=1) do vl
	posloads = filter(>(0),vl)
	isempty(posloads) ? 1.0 : mean(>(LOD),posloads)
end |> x->OffsetArray(vec(x),-1)

# ╔═╡ 0e4bfb8a-a21d-11eb-07d9-51acf99015b3
md"""
### Square array with holes
"""

# ╔═╡ 14d29e14-a21d-11eb-0bb8-4122abd7ed95
arrayholeseff = loadeff("designs/array,holes/")

# ╔═╡ 18aaa1d0-a21d-11eb-225c-ed79f5246ccf
arrayholessens = loadsens("designs/array,holes/")

# ╔═╡ 1b3b908a-a21d-11eb-1ed1-318a1a5fd572
md"""
### Random assignment
"""

# ╔═╡ 203f9d60-a21d-11eb-3d17-47daf91ce8db
randassigneff = loadeff("designs/rand,assign/")

# ╔═╡ 24ae7948-a21d-11eb-0ed9-39ce514e5609
randassignsens = loadsens("designs/rand,assign/")

# ╔═╡ 13866a4a-a21d-11eb-18f0-19efe13c34d1
md"""
### Double-pooling
"""

# ╔═╡ 2d5da91a-a21d-11eb-23c7-1157e133f235
doublepoolingeff = loadeff("designs/double,pooling/")

# ╔═╡ 314dfac0-a21d-11eb-01ae-07c9c529d48a
doublepoolingsens = loadsens("designs/double,pooling/")

# ╔═╡ 0fd98c50-2203-11eb-0f58-d946dcf484b4
md"""
## Figures
"""

# ╔═╡ 222d4b16-7b8c-11eb-0622-797f8f7bba13
POLYORDER = 8

# ╔═╡ 19411dcc-220c-11eb-1b45-4d393e53f5bf
PLTWINDOW = IdentityRange(20:110)

# ╔═╡ a6c3504c-543f-11eb-2f20-791273bc06b8
SHADELIMS = [40,90]

# ╔═╡ e3dcaa08-5036-11eb-2f56-6f1d48bb9188
shadecolor = colorant"gray87"

# ╔═╡ 6c66ebe0-543b-11eb-1619-294095a5ae7a
let compfigs = [96=>(n=96,m=16,q=2),384=>(n=384,m=32,q=2)]
	# Plot colors
	dcolors = distinguishable_colors(7, [RGB(1,1,1)], dropseed=true)
	indivcolor = dcolors[1]
	randassigncolor = dcolors[2]
	doublepoolingcolor = dcolors[6]
	arrayholescolor = dcolors[3]
	hypercolor = Dict(
		16=>dcolors[4],
		32=>dcolors[5],
	)

	# Efficiency gain plots
	effplts = map(compfigs) do (n,hyperdesign)
		plot(rectangle(SHADELIMS,[0,n]),opacity=0.4,linewidth=0,color=shadecolor)

		# Random assignment
		plot!(randassigneff[hyperdesign],linestyle=:solid,color=randassigncolor,linewidth=4)

		# Double-pooling
		plot!(doublepoolingeff[hyperdesign],linestyle=:solid,color=doublepoolingcolor)

		# Array with holes
		plot!(arrayholeseff[(n=n,m=(n==96 ? 20 : 40),q=2)],linestyle=:solid,color=arrayholescolor)

		# HYPER
		plot!(hypereff[hyperdesign],linestyle=:dash,color=hypercolor[hyperdesign.m])

		# Individual
		plot!(indiveff,linestyle=:dot,color=indivcolor)
	end

	# Sensitivity plots
	sensplts = map(compfigs) do (n,hyperdesign)
		plot(rectangle(SHADELIMS,[0,1]),opacity=0.4,linewidth=0,color=shadecolor)

		# Random assignment
		scatter!(randassignsens[hyperdesign],opacity=0.6,color=randassigncolor)
		plot!(polyfit(log10.(prevalence[PLTWINDOW]),randassignsens[hyperdesign][PLTWINDOW],POLYORDER),
			linestyle=:solid,color=randassigncolor)

		# Double-pooling
		scatter!(doublepoolingsens[hyperdesign],opacity=0.6,color=doublepoolingcolor)
		plot!(polyfit(log10.(prevalence[PLTWINDOW]),doublepoolingsens[hyperdesign][PLTWINDOW],POLYORDER),
			linestyle=:solid,color=doublepoolingcolor)

		# Array with holes
		scatter!(arrayholessens[(n=n,m=(n==96 ? 20 : 40),q=2)],opacity=0.6,color=arrayholescolor)
		plot!(polyfit(log10.(prevalence[PLTWINDOW]),arrayholessens[(n=n,m=(n==96 ? 20 : 40),q=2)][PLTWINDOW],POLYORDER),
			linestyle=:solid,color=arrayholescolor)

		# HYPER
		scatter!(hypersens[hyperdesign],opacity=0.6,color=hypercolor[hyperdesign.m])
		plot!(polyfit(log10.(prevalence[PLTWINDOW]),hypersens[hyperdesign][PLTWINDOW],POLYORDER),
			linestyle=:dash,color=hypercolor[hyperdesign.m])

		# Individual
		scatter!(indivsens,opacity=0.6,color=indivcolor)
		plot!(polyfit(log10.(prevalence[PLTWINDOW]),indivsens[PLTWINDOW],POLYORDER),
			linestyle=:dot,color=indivcolor)
	end

	# Figure (titles, axes, annotations)
	title!.(effplts,string.("n=",first.(compfigs)," individuals"))

	plot!.(sensplts,xlims=(20,110),xticks=(20:10:110,["";["$day\n($(round(100*prevalence[day],digits=2))%)" for day in 30:10:100];""]),xlabel="day (prevalence)")
	plot!.(effplts, xlims=(20,110),xticks=(20:10:110,""))

	plot!(effplts[1],ylims=(0, 6*1.14),yticks=[0,1,96/20,96/16])
	plot!(effplts[2],ylims=(0,12*1.14),yticks=[0,1,384/40,384/48,384/32])
	plot!.(sensplts,ylims=(0.6,0.9),yticks=(0.65:0.05:0.85,string.(65:5:85,"%")))
	plot!(effplts[1], ylabel="individuals/test")
	plot!(sensplts[1],ylabel="sensitivity")

	randassignname((n,m,q)) = "Random assignment (m=$m, q=$q)"
	doublepoolingname((n,m,q)) = q == 2 ? "Double-pooling (m=$m)" : "q-pooling (m=$m, q=$q)"
	annotate!(effplts[1],[
		(21.75, 6.1, text(randassignname(last(compfigs[1])),    8, :left, :bottom, randassigncolor)),
		(21.75, 5.9, text(doublepoolingname(last(compfigs[1])), 8, :left, :top, doublepoolingcolor)),
		(21.75, 4.7, text("10x10 array with 4 holes",           8, :left, :top, arrayholescolor)),
		(84.50, 5.3, text(hypername(last(compfigs[1])),         8, :left, :bottom, hypercolor[last(compfigs[1]).m])),
		(21.75, 1.1, text("Individual testing",                 8, :left, :bottom, indivcolor))
	])

	annotate!(effplts[2],[
		(21.75, 12.2, text(randassignname(last(compfigs[2])),    8, :left, :bottom, randassigncolor)),
		(21.75, 11.6, text(doublepoolingname(last(compfigs[2])), 8, :left, :top, doublepoolingcolor)),
		(21.75,  9.4, text("20x20 array with 16 holes",          8, :left, :top, arrayholescolor)),
		(74.50, 10.2, text(hypername(last(compfigs[2])),         8, :left, :bottom, hypercolor[last(compfigs[2]).m])),
		(21.75,  1.2, text("Individual testing",                 8, :left, :bottom, indivcolor))
	])

	plot!.(effplts, bottom_margin=0mm)
	plot!.(sensplts,bottom_margin=5mm,top_margin=-2.5mm)
	plot!.([effplts; sensplts],right_margin=-1mm)

	plot(effplts...,sensplts...,layout=(2,2),link=:x)
	savefig("fig.png"); plot!()
end

# ╔═╡ Cell order:
# ╟─849adbca-08e5-11eb-1956-cd02f2a6be7a
# ╠═95eaf04a-08e5-11eb-043f-27646f1bf207
# ╠═34f9cdda-0901-11eb-1117-099323e644ef
# ╟─289fbed4-21ff-11eb-261a-25eace4cbf4e
# ╟─6539e014-543f-11eb-17c3-991de18e923c
# ╟─eb7bec9e-5453-11eb-0876-c93037e57988
# ╟─210d331a-09c3-11eb-06ed-e1f4ddc023fc
# ╟─fd643c2c-08e5-11eb-0d25-014798a505d1
# ╟─fe7cf842-08e5-11eb-0d07-41ed2902921a
# ╟─3f1b753e-0901-11eb-0060-1f1269588bef
# ╟─3f1c5288-0901-11eb-391f-dd7dd02d3d9b
# ╟─05c6923e-08e6-11eb-1d9b-21afc86a6b04
# ╟─bfbc2e08-0c3b-11eb-10c0-7bd1d3a723b9
# ╟─c6c3f3ca-0c3b-11eb-2593-439620bd9e78
# ╟─cf1b2b1a-0c3b-11eb-2159-6ffeb682d595
# ╟─122151e0-08e6-11eb-18f5-3b277145bc9d
# ╟─47d9acce-2fa6-11eb-0260-7767c0252b8c
# ╟─f8274a32-2fa6-11eb-3836-9bbe0aaf63cf
# ╟─a92db3c4-5436-11eb-181b-3f1a71204b12
# ╟─936c39ce-368e-11eb-0cb2-938376a6c18f
# ╟─72f1ae0a-544c-11eb-0a55-f18f192d64c2
# ╟─21073ca8-3684-11eb-0a25-8d16a9799cf1
# ╟─0e4bfb8a-a21d-11eb-07d9-51acf99015b3
# ╟─14d29e14-a21d-11eb-0bb8-4122abd7ed95
# ╟─18aaa1d0-a21d-11eb-225c-ed79f5246ccf
# ╟─1b3b908a-a21d-11eb-1ed1-318a1a5fd572
# ╟─203f9d60-a21d-11eb-3d17-47daf91ce8db
# ╟─24ae7948-a21d-11eb-0ed9-39ce514e5609
# ╟─13866a4a-a21d-11eb-18f0-19efe13c34d1
# ╟─2d5da91a-a21d-11eb-23c7-1157e133f235
# ╟─314dfac0-a21d-11eb-01ae-07c9c529d48a
# ╟─0fd98c50-2203-11eb-0f58-d946dcf484b4
# ╟─222d4b16-7b8c-11eb-0622-797f8f7bba13
# ╟─19411dcc-220c-11eb-1b45-4d393e53f5bf
# ╟─a6c3504c-543f-11eb-2f20-791273bc06b8
# ╟─e3dcaa08-5036-11eb-2f56-6f1d48bb9188
# ╟─6c66ebe0-543b-11eb-1619-294095a5ae7a


text:
"""
    SystemPath

A union of `PosixPath` and `WindowsPath` which is used for writing
methods that wrap base functionality.
"""
const SystemPath = Union{PosixPath, WindowsPath}

Path() = @static Sys.isunix() ? PosixPath() : WindowsPath()
Path(pieces::Tuple{Vararg{String}}) =
    @static Sys.isunix() ? PosixPath(pieces) : WindowsPath(pieces)

"""
    @__PATH__ -> SystemPath

@__PATH__ expands to a path with the directory part of the absolute path
of the file containing the macro. Returns an empty Path if run from a REPL or
if evaluated by julia -e <expr>.
"""
macro __PATH__()
    p = Path(dirname(string(__source__.file)))
    return p === nothing ? :(Path()) : :($p)
end

"""
    @__FILEPATH__ -> SystemPath

@__FILEPATH__ expands to a path with the absolute file path of the file
containing the macro. Returns an empty Path if run from a REPL or if
evaluated by julia -e <expr>.
"""
macro __FILEPATH__()
    p = Path(string(__source__.file))
    return p === nothing ? :(Path()) : :($p)
end

"""
    @LOCAL(filespec)

Construct an absolute path to `filespec` relative to the source file
containing the macro call.
"""
macro LOCAL(filespec)
    p = join(Path(dirname(string(__source__.file))), Path(filespec))
    return :($p)
end

exists(fp::SystemPath) = ispath(string(fp))
Base.real(fp::P) where {P <: SystemPath} = P(realpath(string(fp)))

#=
The following a descriptive methods for paths
built around stat
=#
Base.stat(fp::SystemPath) = Status(stat(string(fp)))
Base.lstat(fp::SystemPath) = Status(lstat(string(fp)))

"""
    mode(fp::SystemPath) -> Mode

Returns the `Mode` for the specified path.

# Example
```
julia> mode(p"src/FilePathsBase.jl")
-rw-r--r--
```
"""
mode(fp::SystemPath) = stat(fp).mode
Base.size(fp::SystemPath) = stat(fp).size

"""
    modified(fp::SystemPath) -> DateTime

Returns the last modified date for the `path`.

# Example
```
julia> modified(p"src/FilePathsBase.jl")
2017-06-20T04:01:09
```
"""
modified(fp::SystemPath) = stat(fp).mtime

"""
    created(fp::SystemPath) -> DateTime

Returns the creation date for the `path`.

# Example
```
julia> created(p"src/FilePathsBase.jl")
2017-06-20T04:01:09
```
"""
created(fp::SystemPath) = stat(fp).ctime
Base.isdir(fp::SystemPath) = isdir(mode(fp))
Base.isfile(fp::SystemPath) = isfile(mode(fp))
Base.islink(fp::SystemPath) = islink(lstat(fp).mode)
Base.issocket(fp::SystemPath) = issocket(mode(fp))
Base.isfifo(fp::SystemPath) = issocket(mode(fp))
Base.ischardev(fp::SystemPath) = ischardev(mode(fp))
Base.isblockdev(fp::SystemPath) = isblockdev(mode(fp))

"""
    isexecutable(fp::SystemPath) -> Bool

Returns whether the `path` is executable for the current user.
"""
function isexecutable(fp::SystemPath)
    s = stat(fp)
    usr = User()

    return (
        isexecutable(s.mode, :ALL) ||
        isexecutable(s.mode, :OTHER) ||
        ( usr.uid == s.user.uid && isexecutable(s.mode, :USER) ) ||
        ( usr.gid == s.group.gid && isexecutable(s.mode, :GROUP) )
    )
end

"""
    iswritable(fp::AbstractPath) -> Bool

Returns whether the `path` is writable for the current user.
"""
function Base.iswritable(fp::SystemPath)
    s = stat(fp)
    usr = User()

    return (
        iswritable(s.mode, :ALL) ||
        iswritable(s.mode, :OTHER) ||
        ( usr.uid == s.user.uid && iswritable(s.mode, :USER) ) ||
        ( usr.gid == s.group.gid && iswritable(s.mode, :GROUP) )
    )
end

"""
    isreadable(fp::SystemPath) -> Bool

Returns whether the `path` is readable for the current user.
"""
function Base.isreadable(fp::SystemPath)
    s = stat(fp)
    usr = User()

    return (
        isreadable(s.mode, :ALL) ||
        isreadable(s.mode, :OTHER) ||
        ( usr.uid == s.user.uid && isreadable(s.mode, :USER) ) ||
        ( usr.gid == s.group.gid && isreadable(s.mode, :GROUP) )
    )
end

function Base.ismount(fp::SystemPath)
    isdir(fp) || return false
    s1 = lstat(fp)
    # Symbolic links cannot be mount points
    islink(fp) && return false
    s2 = lstat(parent(fp))
    # If a directory and its parent are on different devices,  then the
    # directory must be a mount point
    (s1.device != s2.device) && return true
    (s1.inode == s2.inode) && return true
    false
end

#=
Path Operations
===============================================
The following are methods for actually manipulating the
filesystem.

NOTE: Currently, we are just wrapping base julia functions,
but in the future we'll likely be handling platform specific
code in the implementation instances.

TODO: Document these once we're comfortable with them.
=#

Base.cd(fp::SystemPath) = cd(string(fp))
function Base.cd(fn::Function, dir::SystemPath)
    old = cwd()
    try
        cd(dir)
        fn()
   finally
        cd(old)
    end
end

function Base.mkdir(fp::SystemPath; mode=0o777, recursive=false, exist_ok=false)
    if exists(fp)
        !exist_ok && error("$fp already exists.")
    else
        if hasparent(fp) && !exists(parent(fp))
			if recursive
				mkdir(parent(fp); mode=mode, recursive=recursive, exist_ok=exist_ok)
			else
				error(
					"The parent of $fp does not exist. " *
					"Pass recursive=true to create it."
				)
			end
        end

		mkdir(string(fp), mode=mode)
    end
end

function Base.symlink(src::SystemPath, dest::SystemPath; exist_ok=false, overwrite=false)
    if exists(src)
        if exists(dest) && exist_ok && overwrite
            rm(dest, recursive=true)
        end

        if !exists(dest)
            symlink(string(src), string(dest))
        elseif !exist_ok
            error("$dest already exists.")
        end
    else
        error("$src is not a valid path")
    end
end

function Base.rm(fp::SystemPath; kwargs...)
    rm(string(fp); kwargs...)
end
Base.touch(fp::SystemPath) = touch(string(fp))
function Base.mktemp(parent::SystemPath)
    fp, io = mktemp(string(parent))
    return Path(fp), io
end

Base.mktempdir(parent::SystemPath) = Path(mktempdir(string(parent)))

"""
    chown(fp::SystemPath, user::AbstractString, group::AbstractString; recursive=false)

Change the `user` and `group` of the `fp`.
"""
function Base.chown(fp::SystemPath, user::AbstractString, group::AbstractString; recursive=false)
    @static if Sys.isunix()
        chown_cmd = String["chown"]
        if recursive
            push!(chown_cmd, "-R")
        end
        append!(chown_cmd, String["$(user):$(group)", string(fp)])

        run(Cmd(chown_cmd))
    else
        error("chown is currently not supported on windows.")
    end
end

"""
    chmod(fp::SystemPath, mode::Mode; recursive=false)
    chmod(fp::SystemPath, mode::Integer; recursive=false)
    chmod(fp::SystemPath, user::UIn8=0o0, group::UInt8=0o0, other::UInt8=0o0; recursive=false)
    chmod(fp::SystemPath, symbolic_mode::AbstractString; recursive=false)

Provides various methods for changing the `mode` of a `fp`.

# Examples
```
julia> touch(p"newfile")
Base.Filesystem.File(false, RawFD(-1))

julia> mode(p"newfile")
-rw-r--r--

julia> chmod(p"newfile", 0o755)

julia> mode(p"newfile")
-rwxr-xr-x

julia> chmod(p"newfile", "-x")

julia> mode(p"newfile")
-rw-r--r--

julia> chmod(p"newfile", user=(READ+WRITE+EXEC), group=(READ+EXEC), other=READ)

julia> mode(p"newfile")
-rwxr-xr--

julia> chmod(p"newfile", mode(p"src/FilePathsBase.jl"))

julia> mode(p"newfile")
-rw-r--r--
```
"""
function Base.chmod(fp::SystemPath, mode::Mode; recursive=false)
    chmod_path = string(fp)
    chmod_mode = raw(mode)

    if isdir(fp) && recursive
        for p in readdir(fp)
            chmod(chmod_path, chmod_mode; recursive=recursive)
        end
    end

    chmod(chmod_path, chmod_mode)
end

function Base.chmod(fp::SystemPath, mode::Integer; recursive=false)
    chmod(fp, Mode(mode); recursive=recursive)
end

function Base.chmod(fp::SystemPath; user::UInt8=0o0, group::UInt8=0o0, other::UInt8=0o0, recursive=false)
    chmod(fp, Mode(user=user, group=group, other=other); recursive=recursive)
end

function Base.chmod(fp::SystemPath, symbolic_mode::AbstractString; recursive=false)
    who_char = ['u', 'g', 'o']
    who_actual = [:USER, :GROUP, :OTHER]
    act_char = ['+', '-', '=']
    perm_char = ['r', 'w', 'x']
    perm_actual = [READ, WRITE, EXEC]
    unsupported_perm_char = ['s', 't', 'X', 'u', 'g', 'o']

    tokenized = split(symbolic_mode, act_char)
    if length(tokenized) != 2
        error("Invalid symbolic string expected format <who><action><perm>.")
    end

    who_raw = tokenized[1]
    perm_raw = tokenized[2]

    who = [:ALL]
    perm = 0o0

    for i in 1:3
        if who_char[i] in who_raw
            push!(who, who_actual[i])
        end
    end

    for i in 1:3
        if perm_char[i] in perm_raw
            perm += perm_actual[i]
        end
    end

    for x in unsupported_perm_char
        if x in perm_raw
            error("$x is currently an unsupported permission char for symbolic modes.")
        end
    end

    m = mode(fp)
    new_m = Mode(perm, who...)

    if '+' in symbolic_mode
        chmod(fp, m + new_m; recursive=recursive)
    elseif '-' in symbolic_mode
        chmod(fp, m - new_m; recursive=recursive)
    elseif '=' in symbolic_mode
        chmod(fp, new_m; recursive=recursive)
    else
        error("No valid action found in symbolic mode string.")
    end
end

Base.open(fp::SystemPath, args...) = open(string(fp), args...)
function Base.open(f::Function, fp::SystemPath, args...; kwargs...)
    open(f, string(fp), args...; kwargs...)
end

Base.read(fp::SystemPath) = read(string(fp))
function Base.write(fp::SystemPath, x::Union{String, Vector{UInt8}}, mode="w")
    open(fp, mode) do f
        write(f, x)
    end
end

Base.readdir(fp::SystemPath) = readdir(string(fp))
Base.download(url::AbstractString, dest::SystemPath) = download(url, string(dest))
Base.readlink(fp::SystemPath) = Path(readlink(string(fp)))


text:
# Unit tests for integrations/dftk_integration.jl

@testset "integrations/dftk_integration.jl" begin
    # TODO
end


text:
module Substance_db_access_test
"""
Module to test the methods of module `Functions` ment to acess the project substance database and retireve data from it
"""

include("..\\src\\Analysis\\Utils\\FunctionsDB.jl")

# DAF
#   line 292
restex = Functions.texture_extract("sand", ["tot_por", "c_water_avg", "ief", "por_eff", "grain"])
#   line 307
ressub = Functions.substance_extract(4, ["c_henry", "koc_kd"])

# Runoffs
#   line 519
rescnl = Functions.cn_list_extract()

# Plumes
#   line 80
resair = Functions.air_extract("a", "c", ["sigmay1", "sigmay2", "sigmayexp", "sigmaz1", "sigmaz2", "sigmazexp"])
#   line 209
ressu2 = Functions.substance_extract(4, [ "sf_ing", "sf_inal", "iur", "rfd_ing", "rfd_inal", "rfc"])

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
# Buffer
# ======

# Data Layout
# -----------
#
# Buffered data are stored in `data` and three position fields are used to keep
# track of marked data, buffered data and margin.
#
#             marked      buffer      margin
#          |<-------->||<-------->||<-------->|
#     |....xxxxxxxxxxxxXXXXXXXXXXXX...........|
#     ^    ^           ^           ^          ^
#     1    markpos     bufferpos   marginpos  lastindex(data)
#
# `markpos` is positive iff there are marked data; otherwise it is set to zero.
# `markpos` ≤ `bufferpos` ≤ `marginpos` must hold whenever possible.

mutable struct Buffer
    # data and positions (see above)
    data::Vector{UInt8}
    markpos::Int
    bufferpos::Int
    marginpos::Int

    # the number of total bytes passed through this buffer
    total::Int64

    Buffer(size::Integer) = new(Vector{UInt8}(undef, size), 0, 1, 1, 0)
    Buffer(data::Vector{UInt8}) = new(data, 0, 1, length(data)+1, 0)
end

Buffer(data::Base.CodeUnits{UInt8}) = Buffer(Vector{UInt8}(data))
Base.length(buf::Buffer) = length(buf.data)
bufferptr(buf::Buffer) = pointer(buf.data, buf.bufferpos)
buffersize(buf::Buffer) = buf.marginpos - buf.bufferpos
buffermem(buf::Buffer) = Memory(bufferptr(buf), buffersize(buf))
marginptr(buf::Buffer) = pointer(buf.data, buf.marginpos)
marginsize(buf::Buffer) = lastindex(buf.data) - buf.marginpos + 1
marginmem(buf::Buffer) = Memory(marginptr(buf), marginsize(buf))
ismarked(buf::Buffer) = buf.markpos != 0
mark!(buf::Buffer) = buf.markpos = buf.bufferpos
unmark!(buf::Buffer) = if buf.markpos==0; false else; buf.markpos=0; true end

function reset!(buf::Buffer)
    @assert buf.markpos > 0
    buf.bufferpos = buf.markpos
    buf.markpos = 0
    return buf.bufferpos
end

consumed!(buf::Buffer, n::Integer) = (buf.bufferpos += n; buf) # Notify that `n` bytes are consumed from `buf`.
supplied!(buf::Buffer, n::Integer) = (buf.marginpos += n; buf) # Notify that `n` bytes are supplied to `buf`.
consumed2!(buf::Buffer, n::Integer) = (buf.bufferpos += n; buf.total += n; buf)
supplied2!(buf::Buffer, n::Integer) = (buf.marginpos += n; buf.total += n; buf)
initbuffer!(buf::Buffer) = (buf.markpos = 0; buf.bufferpos = buf.marginpos = 1; buf.total = 0; buf) # Discard buffered data and initialize positions.
emptybuffer!(buf::Buffer) = (buf.marginpos = buf.bufferpos; buf) # Remove all buffered data.

# Make margin with ≥`minsize` and return the size of it.
function makemargin!(buf::Buffer, minsize::Integer)
    @assert minsize ≥ 0
    if buffersize(buf) == 0 && buf.markpos == 0
        buf.bufferpos = buf.marginpos = 1
    end
    if marginsize(buf) < minsize
        # shift data to left
        if buf.markpos == 0
            datapos = buf.bufferpos
            datasize = buffersize(buf)
        else
            datapos = buf.markpos
            datasize = buf.marginpos - buf.markpos
        end
        copyto!(buf.data, 1, buf.data, datapos, datasize)
        shift = datapos - 1
        if buf.markpos > 0
            buf.markpos -= shift
        end
        buf.bufferpos -= shift
        buf.marginpos -= shift
    end
    if marginsize(buf) < minsize
        # expand data buffer
        resize!(buf.data, buf.marginpos + minsize - 1)
    end
    @assert marginsize(buf) ≥ minsize
    return marginsize(buf)
end

readbyte!(buf::Buffer) = (b = buf.data[buf.bufferpos]; consumed!(buf, 1); b) # Read a byte.
writebyte!(buf::Buffer, b::UInt8) = (buf.data[buf.marginpos] = b; supplied!(buf, 1); 1) # Write a byte.
skipbuffer!(buf::Buffer, n::Integer) = ((@assert n ≤ buffersize(buf)); consumed!(buf, n); buf) # Skip `n` bytes in the buffer.

# Take the ownership of the marked data.
function takemarked!(buf::Buffer)
    @assert buf.markpos > 0
    sz = buf.marginpos - buf.markpos
    copyto!(buf.data, 1, buf.data, buf.markpos, sz)
    initbuffer!(buf)
    return resize!(buf.data, sz)
end

# Copy data from `data` to `buf`.
function copydata!(buf::Buffer, data::Ptr{UInt8}, nbytes::Integer)
    makemargin!(buf, nbytes)
    unsafe_copyto!(marginptr(buf), data, nbytes)
    supplied!(buf, nbytes)
    return buf
end

# Copy data from `buf` to `data`.
function copydata!(data::Ptr{UInt8}, buf::Buffer, nbytes::Integer)
    # NOTE: It's caller's responsibility to ensure that the buffer has at least
    # nbytes.
    @assert buffersize(buf) ≥ nbytes
    unsafe_copyto!(data, bufferptr(buf), nbytes)
    consumed!(buf, nbytes)
    return data
end

# Insert data to the current buffer.
function insertdata!(buf::Buffer, data::Ptr{UInt8}, nbytes::Integer)
    makemargin!(buf, nbytes)
    copyto!(buf.data, buf.bufferpos + nbytes, buf.data, buf.bufferpos, buffersize(buf))
    unsafe_copyto!(bufferptr(buf), data, nbytes)
    supplied!(buf, nbytes)
    return buf
end

# Find the first occurrence of a specific byte.
function findbyte(buf::Buffer, byte::UInt8)
    p = ccall(
        :memchr,
        Ptr{UInt8},
        (Ptr{UInt8}, Cint, Csize_t),
        pointer(buf.data, buf.bufferpos), byte, buffersize(buf))
    if p == C_NULL
        return marginptr(buf)
    else
        return p
    end
end


text:
plot_line(args...; kwargs...) = @info "plot_line is deprecated, use Plots.plot instead."

@recipe function plot_line_backend(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AA{<:AbstractControlVolume,1}}

    # solution
    pltx = KS.pSpace.x[1:KS.pSpace.nx]
    plty = zeros(KS.pSpace.nx, 3)
    for i in eachindex(pltx)
        for j = 1:2
            plty[i, j] = ctr[i].prim[j]
        end

        plty[i, 3] = 1.0 / ctr[i].prim[end]
    end

    # attributes
    xguide --> "x"
    :linewidth --> 1.5

    @series begin
        label := "ρ"
        pltx, plty[:, 1]
    end

    if ctr[1].w isa AA
        @series begin
            label := "u"
            pltx, plty[:, 2]
        end
        @series begin
            label := "T"
            pltx, plty[:, 3]
        end
    end

    # user-defined
    c = get(plotattributes, :linewidth, :auto)

    return nothing

end

plot_contour(args...; kwargs...) = @info "plot_line is deprecated, use Plots.plot instead."

@recipe function plot_contour_backend(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AA{<:AbstractControlVolume,2}}

    pltx = KS.ps.x[1:KS.ps.nx, 1]
    plty = KS.ps.y[1, 1:KS.ps.ny]

    sol = zeros(size(ctr[1].w, 1), KS.pSpace.nx, KS.pSpace.ny)
    for i in axes(sol, 2)
        for j in axes(sol, 3)
            for k = 1:size(sol, 1)-1
                sol[k, i, j] = ctr[i, j].prim[k]
            end
            sol[end, i, j] = 1.0 / ctr[i, j].prim[end]
        end
    end

    layout := (2, 2)
    c = get(plotattributes, :inferno, :auto)
    for (i, l) in enumerate(("ρ", "u", "v", "T"))
        @series begin
            subplot := i
            xguide := "x"
            yguide := "y"
            fill := true
            seriescolor := c
            pltx, plty, sol[i, :, :]'
        end
    end

    return nothing

end


text:
# Load LanguageServer from the project next to this file
## Save old load path
old_load_path = copy(LOAD_PATH)
push!(empty!(LOAD_PATH), @__DIR__)
## Load packages
using LanguageServer, SymbolServer
## Restore old load path
append!(empty!(LOAD_PATH), old_load_path)

# Figure out the active project
## This configuration is a good default
project_path = let
    dirname(something(
        ## 1. Finds an explicitly set project (JULIA_PROJECT)
        Base.load_path_expand((
            p = get(ENV, "JULIA_PROJECT", nothing);
            p === nothing ? nothing : isempty(p) ? nothing : p
        )),
        ## 2. Look for a Project.toml file in the current working directory,
        ##    or parent directories, with $HOME as an upper boundary
        Base.current_project(),
        ## 3. First entry in the load path
        get(Base.load_path(), 1, nothing),
        ## 4. Fallback to default global environment,
        ##    this is more or less unreachable
        Base.load_path_expand("@v#.#"),
    ))
end

# Depot path for the server to index (empty uses default).
depot_path = get(ENV, "JULIA_DEPOT_PATH", "")

# Start the server
@info "Running julia language server" VERSION project_path depot_path
server = LanguageServer.LanguageServerInstance(stdin, stdout, project_path, depot_path)
server.runlinter = true
run(server)

max_stars_repo_name: bamzi/dotconfig

max_issues_repo_name: bamzi/dotconfig

max_forks_repo_name: bamzi/dotconfig


text:
function get_value(::Type{JacobiPolynomials}, arg0::jint, arg1::jint, arg2::jint, arg3::FieldGradient)
    return jcall(JacobiPolynomials, "getValue", FieldGradient, (jint, jint, jint, FieldGradient), arg0, arg1, arg2, arg3)
end

function get_value(::Type{JacobiPolynomials}, arg0::jint, arg1::jint, arg2::jint, arg3::Gradient)
    return jcall(JacobiPolynomials, "getValue", Gradient, (jint, jint, jint, Gradient), arg0, arg1, arg2, arg3)
end


text:
module IrisSCPA

end # module


text:
module TestSubDataFrame
    using Compat, Compat.Test, DataFrames

    @testset "view -- DataFrame" begin
        df = DataFrame(x = 1:10, y = 1.0:10.0)
        @test view(df, 1) == head(df, 1)
        @test view(df, UInt(1)) == head(df, 1)
        @test view(df, BigInt(1)) == head(df, 1)
        @test view(df, 1:2) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8))) == head(df, 2)
        @test view(df, [1, 2]) == head(df, 2)
        @test view(df, 1, :x) == head(df[[:x]], 1)
        @test view(df, 1:2, :x) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), :x) == head(df[[:x]], 2)
        @test view(df, [1, 2], :x) == head(df[[:x]], 2)
        @test view(df, 1, 1) == head(df[[:x]], 1)
        @test view(df, 1:2, 1) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), 1) == head(df[[:x]], 2)
        @test view(df, [1, 2], 1) == head(df[[:x]], 2)
        @test view(df, 1, [:x, :y]) == head(df, 1)
        @test view(df, 1:2, [:x, :y]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [:x, :y]) == head(df, 2)
        @test view(df, [1, 2], [:x, :y]) == head(df, 2)
        @test view(df, 1, [1, 2]) == head(df, 1)
        @test view(df, 1:2, [1, 2]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [1, 2]) == head(df, 2)
        @test view(df, [1, 2], [1, 2]) == head(df, 2)
        @test view(df, 1, trues(2)) == head(df, 1)
        @test view(df, 1:2, trues(2)) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), trues(2)) == head(df, 2)
        @test view(df, [1, 2], trues(2)) == head(df, 2)
        @test view(df, Integer[1, 2]) == head(df, 2)
        @test view(df, UInt[1, 2]) == head(df, 2)
        @test view(df, BigInt[1, 2]) == head(df, 2)
        @test view(df, Union{Int, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{Integer, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{UInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{BigInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, :) == df
        @test view(df, :, :) == df
        @test view(df, 1, :) == head(df, 1)
        @test view(df, :, 1) == df[:, [1]]
        @test_throws MissingException view(df, [missing, 1])
    end

    @testset "view -- SubDataFrame" begin
        df = view(DataFrame(x = 1:10, y = 1.0:10.0), 1:10)
        @test view(df, 1) == head(df, 1)
        @test view(df, UInt(1)) == head(df, 1)
        @test view(df, BigInt(1)) == head(df, 1)
        @test view(df, 1:2) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8))) == head(df, 2)
        @test view(df, [1, 2]) == head(df, 2)
        @test view(df, 1, :x) == head(df[[:x]], 1)
        @test view(df, 1:2, :x) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), :x) == head(df[[:x]], 2)
        @test view(df, [1, 2], :x) == head(df[[:x]], 2)
        @test view(df, 1, 1) == head(df[[:x]], 1)
        @test view(df, 1:2, 1) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), 1) == head(df[[:x]], 2)
        @test view(df, [1, 2], 1) == head(df[[:x]], 2)
        @test view(df, 1, [:x, :y]) == head(df, 1)
        @test view(df, 1:2, [:x, :y]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [:x, :y]) == head(df, 2)
        @test view(df, [1, 2], [:x, :y]) == head(df, 2)
        @test view(df, 1, [1, 2]) == head(df, 1)
        @test view(df, 1:2, [1, 2]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [1, 2]) == head(df, 2)
        @test view(df, [1, 2], [1, 2]) == head(df, 2)
        @test view(df, 1, trues(2)) == head(df, 1)
        @test view(df, 1:2, trues(2)) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), trues(2)) == head(df, 2)
        @test view(df, [1, 2], trues(2)) == head(df, 2)
        @test view(df, Integer[1, 2]) == head(df, 2)
        @test view(df, UInt[1, 2]) == head(df, 2)
        @test view(df, BigInt[1, 2]) == head(df, 2)
        @test view(df, Union{Int, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{Integer, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{UInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{BigInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, :) == df
        @test view(df, :, :) == df
        @test view(df, 1, :) == head(df, 1)
        @test view(df, :, 1) == df[:, [1]]
        @test_throws MissingException view(df, [missing, 1])
    end
end


text:
"""
    multinomial_hyper(α::AbstractArray{Float32,1})

[Dirichlet Distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution)
"""
struct multinomial_hyper <: distribution_hyper_params
    α::AbstractArray{Float32,1}
end

mutable struct multinomial_sufficient_statistics <: sufficient_statistics
    N::Float32
    points_sum::AbstractArray{Float32,1}
end

function calc_posterior(
    prior::multinomial_hyper,
    suff_statistics::multinomial_sufficient_statistics,
)
    if suff_statistics.N == 0
        return prior
    end
    return multinomial_hyper(prior.α + suff_statistics.points_sum)
end

function sample_distribution(hyperparams::multinomial_hyper)
    return multinomial_dist(log.(rand(Dirichlet(Float64.(hyperparams.α)))))
end

function create_sufficient_statistics(
    hyper::multinomial_hyper,
    posterior::multinomial_hyper,
    points::AbstractArray{Float32,2},
    pts_to_group = 0,
)
    pts = copy(points)
    points_sum = sum(pts, dims = 2)[:]
    S = pts * pts'
    return multinomial_sufficient_statistics(size(points, 2), points_sum)
end

function log_marginal_likelihood(
    hyper::multinomial_hyper,
    posterior_hyper::multinomial_hyper,
    suff_stats::multinomial_sufficient_statistics,
)
    D = length(suff_stats.points_sum)
    logpi = log(pi)
    val =
        logabsgamma(sum(hyper.α))[1] - logabsgamma(sum(posterior_hyper.α))[1] + sum(
            (x -> logabsgamma(x)[1]).(posterior_hyper.α) -
            (x -> logabsgamma(x)[1]).(hyper.α),
        )
    return val
end

function aggregate_suff_stats(
    suff_l::multinomial_sufficient_statistics,
    suff_r::multinomial_sufficient_statistics,
)
    return multinomial_sufficient_statistics(
        suff_l.N + suff_r.N,
        suff_l.points_sum + suff_r.points_sum,
    )
end


text:
""" :( """
GABLS_u(::Type{FT}) where {FT} = z -> FT(8.0)
""" :( """
GABLS_v(::Type{FT}) where {FT} = z -> FT(0.0)

""" :( """
GABLS_θ_liq_ice(::Type{FT}) where {FT} =
    z -> if z <= 100.0
        FT(265)
    else
        FT(265) + (z - 100) * FT(0.01)
    end

""" :( """
GABLS_q_tot(::Type{FT}) where {FT} = z -> FT(0)

""" :( """
GABLS_tke(::Type{FT}) where {FT} =
    z -> if z <= 250.0
        FT(0.4) * (1 - z / 250) * (1 - z / 250) * (1 - z / 250)
    else
        FT(0)
    end

""" :( """
GABLS_geostrophic_ug(::Type{FT}) where {FT} = z -> FT(8)
""" :( """
GABLS_geostrophic_vg(::Type{FT}) where {FT} = z -> FT(0)


text:
README(name) = """
# $name

This project was made with the [Lynx.jl](https://github.com/jorge-brito/Lynx.jl) package.

## Preview

![Preview](preview.gif)

## License

This project is part of the [Lynx Gallery](https://github.com/jorge-brito/LynxGallery) repository
which is licensed under the [MIT license](https://github.com/jorge-brito/LynxGallery.jl/blob/master/LICENSE).
"""

ENTRY_POINT(name) = """
module $name

using Lynx
using Luxor

@info "Running '$name'"

Lynx.init("$name", 400, 400)

function update(dt)
    # draw here
    background("#111")
end

run!(update, await=true)

end # $name
"""

VSCODE_TASKS(name) = """
{
    "version": "2.0.0",
    "tasks": [
        {
            "label": "Run project",
            "type": "shell",
            "command": "julia --project=.  -e 'push!(LOAD_PATH, joinpath(pwd(), \\"..\\")); include(\\"src/$name.jl\\")'",
            "problemMatcher": []
        }
    ],
}
"""

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
"""
    isosurface(θ,T,Γ)

Depth of isosurface θ=T
"""
function isosurface(θ,T,Γ)
    d=NaN*similar(θ[:,1])
    nr=size(θ,2)
    for j=1:size(d,1)
        for k=1:nr-1
            i=findall(isnan.(d[j]).&(θ[j,k].>T).&(θ[j,k+1].<=T))
            a=(θ[j,k][i] .- T)./(θ[j,k][i] .- θ[j,k+1][i])
            d[j][i]=(1 .- a).*Γ.RC[k] + a.*Γ.RC[k+1]
            i=findall(isnan.(d[j]).&(θ[j,k].<=T).&(θ[j,k+1].>T))
            a=(θ[j,k+1][i] .- T)./(θ[j,k+1][i] .- θ[j,k][i])
            d[j][i]=(1 .- a).*Γ.RC[k+1] + a.*Γ.RC[k]
        end
    end
    return d
end