text:
import Random: rand
import Base: isapprox, copy, deepcopy

"""
    Gaussian{Tm, TP}

A (multivariate) Gaussian with mean vector `m` and variance matrix `P`.
This doesn't currently conform to Distributions.jl standards. Work to make this happen
would be welcomed.

It was necessary to implement this a year or so ago for AD-related reasons. It's quite
possible that in the intervening period of time things have improved and this type is no
longer necessary in addition to the `MvNormal` type in `Distributions`. I've not had the
time to remove it though.
"""
struct Gaussian{Tm, TP}
    m::Tm
    P::TP
end

dim(x::Gaussian) = length(x.m)

AbstractGPs.mean(x::Gaussian) = Zygote.literal_getfield(x, Val(:m))

AbstractGPs.cov(x::Gaussian) = Zygote.literal_getfield(x, Val(:P))

AbstractGPs.var(x::Gaussian{<:AbstractVector}) = diag(cov(x))

AbstractGPs.var(x::Gaussian{<:Real}) = cov(x)

get_fields(x::Gaussian) = mean(x), cov(x)

Random.rand(rng::AbstractRNG, x::Gaussian) = vec(rand(rng, x, 1))

function Random.rand(rng::AbstractRNG, x::Gaussian{Tm}) where {Tm<:SVector}
    ε = randn(rng, Tm)
    return mean(x) + cholesky(Symmetric(cov(x) + UniformScaling(1e-12))).U' * ε
end

function Random.rand(rng::AbstractRNG, x::Gaussian, S::Int)
    P = cov(x) + UniformScaling(1e-12)
    return mean(x) .+ cholesky(Symmetric(P)).U' * randn(rng, length(mean(x)), S)
end

function AbstractGPs.logpdf(x::Gaussian, y::AbstractVector{<:Real})
    return first(logpdf(x, reshape(y, :, 1)))
end

function AbstractGPs.logpdf(x::Gaussian, Y::AbstractMatrix{<:Real})
    μ, C = mean(x), cholesky(Symmetric(cov(x)))
    T = promote_type(eltype(μ), eltype(C), eltype(Y))
    return -((size(Y, 1) * T(log(2π)) + logdet(C)) .+ diag_Xt_invA_X(C, Y .- μ)) ./ 2
end

Base.:(==)(x::Gaussian, y::Gaussian) = mean(x) == mean(y) && cov(x) == cov(y)

function Base.isapprox(x::Gaussian, y::Gaussian; kwargs...)
    return isapprox(mean(x), mean(y); kwargs...) && isapprox(cov(x), cov(y); kwargs...)
end

function AbstractGPs.marginals(x::Gaussian{T, T}) where {T<:Real}
    return AbstractGPs.Normal{T}(mean(x), sqrt(cov(x)))
end

function AbstractGPs.marginals(x::Gaussian{<:AbstractVector, <:AbstractMatrix})
    return AbstractGPs.Normal.(mean(x), sqrt.(diag(cov(x))))
end

storage_type(x::Gaussian{<:SVector{D, T}}) where {D, T<:Real} = SArrayStorage(T)

storage_type(gmm::Gaussian{<:Vector{T}}) where {T<:Real} = ArrayStorage(T)

storage_type(x::Gaussian{T}) where {T<:Real} = ScalarStorage(T)

function Zygote._pullback(::AContext, ::Type{<:Gaussian}, m, P)
    Gaussian_pullback(Δ::Nothing) = (nothing, nothing, nothing)
    Gaussian_pullback(Δ) = (nothing, Δ.m, Δ.P)
    return Gaussian(m, P), Gaussian_pullback
end

Base.length(x::Gaussian) = 0

# Zero-adjoint initialisation for the benefit of `scan`.
_get_zero_adjoint(x::Gaussian) = (m=_get_zero_adjoint(mean(x)), P=_get_zero_adjoint(cov(x)))


text:
@testset "reduce" begin
    m = GBMatrix([[1,2,3] [4,5,6] [7,8,9]])
    reduce(max, m, dims=2) == reduce(Monoids.MAX_MONOID, m) #this only works for dense
    reduce(Monoids.MAX_MONOID, m, dims=(1,2)) == 9
    @test_throws ArgumentError reduce(BinaryOps.TIMES, m)
end

ext: jl


text:
#!/usr/bin/env julia

"""
	Consensus checking for variant positions across samples
"""
file = ARGS[1] # .vcf/.vcf.gz file
include("../src/read.jl")

vcf = read(file)
header = []
names = split(vcf[1])[5:end]  # assumes the first line as header
push!(header, "Position")
push!(header, "Consensus")
append!(header, names)

for col in header
	print(col, '\t')
end
println()

for line in vcf[2:end]
	arr =[]
	cols = split(line)
	print(cols[2], '\t')
	alts = []
	for sample in cols[5:end]
		altRead = parse(Int64, split(sample, ':')[1])
		if altRead >= 3
			push!(alts, 1)
		else
			push!(alts, 0)
		end
	end

	@printf("%6.3f", consensus)
	for i in alts
		print('\t', i)
	end
	println()
end


text:
#!/usr/bin/env julia
""" Kullback-Leibler divergence functions and klUCB utilities, in naive Julia code.

- How to use it? From Julia, it's easy:

```
julia> using KullbackLeibler
julia> KullbackLeibler.klBern(0.5, 0.5)
0.0
julia> KullbackLeibler.klBern(0.1, 0.9)
1.757779...
julia> KullbackLeibler.klBern(0.9, 0.1)  # And this KL is symmetric
1.757779...
julia> KullbackLeibler.klBern(0.4, 0.5)
0.020135...
julia> KullbackLeibler.klBern(0.01, 0.99)
4.503217...
```

- You can also use it with the [`Distributions`](https://juliastats.github.io/Distributions.jl/) module:

```julia
julia> using Distributions
julia> Bern1 = Distributions.Bernoulli(0.33)
Distributions.Bernoulli{Float64}(p=0.33)
julia> Bern2 = Distributions.Bernoulli(0.42)
Distributions.Bernoulli{Float64}(p=0.42)
julia> ex_kl_1 = KL( Bern1, Bern2 )  # Calc KL divergence for Bernoulli R.V
0.017063...
julia> klBern(0.33, 0.42)  # same!
0.017063...
```

- Note: there is also a Python version, if you need. (Cf. https://github.com/Naereen/Kullback-Leibler-divergences-and-kl-UCB-indexes/)

- MIT Licensed, (C) 2018 Lilian Besson (Naereen), https://GitHub.com/Naereen/KullbackLeibler.jl

- Cf. https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence
- Reference: [Filippi, Cappé & Garivier - Allerton, 2011](https://arxiv.org/pdf/1004.5229.pdf), [Garivier & Cappé, 2011](https://arxiv.org/pdf/1102.2490.pdf), and [Kullback & Leibler, 1951](http://www.jstor.org/stable/2236703).
"""

__author__ = "Lilian Besson"
__version__ = "0.1"

# Name of module
module KullbackLeibler

# List of exported functions
export klBern, klBin, klPoisson, klExp, klGamma, klNegBin, klGauss, KL, klucb, klucbBern, klucbGauss, klucbPoisson, klucbExp, klucbGamma

eps = 1e-15  #: Threshold value: everything in [0, 1] is truncated to [eps, 1 - eps]

# --- Simple Kullback-Leibler divergence for known distributions

doc"""
    function klBern(x, y)

Kullback-Leibler divergence for Bernoulli distributions. https://en.wikipedia.org/wiki/Bernoulli_distribution#Kullback.E2.80.93Leibler_divergence

$$\mathrm{KL}(\mathcal{B}(x), \mathcal{B}(y)) = x \log(\frac{x}{y}) + (1-x) \log(\frac{1-x}{1-y}).$$

```julia
julia> klBern(0.5, 0.5)
0.0
julia> klBern(0.1, 0.9)
1.757779...
julia> klBern(0.9, 0.1)  # And this KL is symmetric
1.757779...
julia> klBern(0.4, 0.5)
0.020135...
julia> klBern(0.01, 0.99)
4.503217...
```

- Special values:

```julia
julia> klBern(0, 1)  # Should be +Inf, but 0 --> eps, 1 --> 1 - eps
34.539575...
```
"""
function klBern(x, y)
    x = min(max(x, eps), 1 - eps)
    y = min(max(y, eps), 1 - eps)
    return x * log(x / y) + (1 - x) * log((1 - x) / (1 - y))
end

doc"""
    function klBin(x, y, n)

Kullback-Leibler divergence for Binomial distributions. https://math.stackexchange.com/questions/320399/kullback-leibner-divergence-of-binomial-distributions

- It is simply the n times `klBern` on x and y.

$$\mathrm{KL}(\mathrm{Bin}(x, n), \mathrm{Bin}(y, n)) = n \times \left(x \log(\frac{x}{y}) + (1-x) \log(\frac{1-x}{1-y}) \right).$$

- **Warning**: The two distributions must have the same parameter n, and x, y are p, q in (0, 1).

```julia
julia> klBin(0.5, 0.5, 10)
0.0
julia> klBin(0.1, 0.9, 10)
17.57779...
julia> klBin(0.9, 0.1, 10)  # And this KL is symmetric
17.57779...
julia> klBin(0.4, 0.5, 10)
0.20135...
julia> klBin(0.01, 0.99, 10)
45.03217...
```

- Special values:

```julia
julia> klBin(0, 1, 10)  # Should be +Inf, but 0 --> eps, 1 --> 1 - eps
345.39575...
```
"""
function klBin(x, y, n)
    x = min(max(x, eps), 1 - eps)
    y = min(max(y, eps), 1 - eps)
    return n * (x * log(x / y) + (1 - x) * log((1 - x) / (1 - y)))
end

doc"""
    function klPoisson(x, y)

Kullback-Leibler divergence for Poison distributions. https://en.wikipedia.org/wiki/Poisson_distribution#Kullback.E2.80.93Leibler_divergence

$$\mathrm{KL}(\mathrm{Poisson}(x), \mathrm{Poisson}(y)) = y - x + x \times \log(\frac{x}{y}).$$

```julia
julia> klPoisson(3, 3)
0.0
julia> klPoisson(2, 1)
0.386294...
julia> klPoisson(1, 2)  # And this KL is non-symmetric
0.306852...
julia> klPoisson(3, 6)
0.920558...
julia> klPoisson(6, 8)
0.273907...
```

- Special values:

```julia
julia> klPoisson(1, 0)  # Should be +Inf, but 0 --> eps, 1 --> 1 - eps
33.538776...
julia> klPoisson(0, 0)
0.0
```
"""
function klPoisson(x, y)
    x = max(x, eps)
    y = max(y, eps)
    return y - x + x * log(x / y)
end

doc"""
    function klExp(x, y)

Kullback-Leibler divergence for exponential distributions. https://en.wikipedia.org/wiki/Exponential_distribution#Kullback.E2.80.93Leibler_divergence

$$\mathrm{KL}(\mathrm{Exp}(x), \mathrm{Exp}(y)) = \begin{cases}
\frac{x}{y} - 1 - \log(\frac{x}{y}) & \text{if} x > 0, y > 0\\
+\infty & \text{otherwise}
\end{cases}$$

```julia
julia> klExp(3, 3)
0.0
julia> klExp(3, 6)
0.193147...
julia> klExp(1, 2)  # Only the proportion between x and y is used
0.193147...
julia> klExp(2, 1)  # And this KL is non-symmetric
0.306852...
julia> klExp(4, 2)  # Only the proportion between x and y is used
0.306852...
julia> klExp(6, 8)
0.037682...
```

- x, y have to be positive:

```julia
julia> klExp(-3, 2)
Inf
julia> klExp(3, -2)
Inf
julia> klExp(-3, -2)
Inf
```
"""
function klExp(x, y)
    if (x <= 0) || (y <= 0)
        return Inf
    else
        x = max(x, eps)
        y = max(y, eps)
        return x / y - 1 - log(x / y)
    end
end

doc"""
    function klGamma(x, y, a=1)

Kullback-Leibler divergence for gamma distributions. https://en.wikipedia.org/wiki/Gamma_distribution#Kullback.E2.80.93Leibler_divergence

- It is simply the a times `klExp` on x and y.

.. math::

    \mathrm{KL}(\Gamma(x, a), \Gamma(y, a)) = \begin{cases}
    a \times \left( \frac{x}{y} - 1 - \log(\frac{x}{y}) \right) & \text{if} x > 0, y > 0\\
    +\infty & \text{otherwise}
    \end{cases}

- **Warning**: The two distributions must have the same parameter a.

```julia
julia> klGamma(3, 3)
0.0
julia> klGamma(3, 6)
0.193147...
julia> klGamma(1, 2)  # Only the proportion between x and y is used
0.193147...
julia> klGamma(2, 1)  # And this KL is non-symmetric
0.306852...
julia> klGamma(4, 2)  # Only the proportion between x and y is used
0.306852...
julia> klGamma(6, 8)
0.037682...
```

- x, y have to be positive:

```julia
julia> klGamma(-3, 2)
Inf
julia> klGamma(3, -2)
Inf
julia> klGamma(-3, -2)
Inf
```
"""
function klGamma(x, y, a=1)
    if (x <= 0) || (y <= 0)
        return Inf
    else
        x = max(x, eps)
        y = max(y, eps)
        return a * (x / y - 1 - log(x / y))
    end
end

doc"""
    function klNegBin(x, y, r=1)

Kullback-Leibler divergence for negative binomial distributions. https://en.wikipedia.org/wiki/Negative_binomial_distribution

$$\mathrm{KL}(\mathrm{NegBin}(x, r), \mathrm{NegBin}(y, r)) = r \times \log((r + x) / (r + y)) - x \times \log(y \times (r + x) / (x \times (r + y))).$$

- **Warning**: The two distributions must have the same parameter r.

```julia
julia> klNegBin(0.5, 0.5)
0.0
julia> klNegBin(0.1, 0.9)
-0.711611...
julia> klNegBin(0.9, 0.1)  # And this KL is non-symmetric
2.0321564...
julia> klNegBin(0.4, 0.5)
-0.130653...
julia> klNegBin(0.01, 0.99)
-0.717353...
```

- Special values:

```julia
julia> klBern(0, 1)  # Should be +Inf, but 0 --> eps, 1 --> 1 - eps
    34.539575...
```

- With other values for `r`:

```julia
julia> klNegBin(0.5, 0.5, r=2)
0.0
julia> klNegBin(0.1, 0.9, r=2)
-0.832991...
julia> klNegBin(0.1, 0.9, r=4)
-0.914890...
julia> klNegBin(0.9, 0.1, r=2)  # And this KL is non-symmetric
2.3325528...
julia> klNegBin(0.4, 0.5, r=2)
-0.154572...
julia> klNegBin(0.01, 0.99, r=2)
-0.836257...
```
"""
function klNegBin(x, y, r=1)
    x = max(x, eps)
    y = max(y, eps)
    return r * log((r + x) / (r + y)) - x * log(y * (r + x) / (x * (r + y)))
end

doc"""
    function klGauss(x, y, sig2x=0.25, sig2y=0.25)

Kullback-Leibler divergence for Gaussian distributions of means ``x`` and ``y`` and variances ``sig2x`` and ``sig2y``, $\nu_1 = \mathcal{N}(x, \sigma_x^2)$ and $\nu_2 = \mathcal{N}(y, \sigma_x^2)$:

$$\mathrm{KL}(\nu_1, \nu_2) = \frac{(x - y)^2}{2 \sigma_y^2} + \frac{1}{2}\left( \frac{\sigma_x^2}{\sigma_y^2} - 1 \log\left(\frac{\sigma_x^2}{\sigma_y^2}\right) \right).$$

See https://en.wikipedia.org/wiki/Normal_distribution#Other_properties

- By default, sig2y is assumed to be sig2x (same variance).

```julia
julia> klGauss(3, 3)
0.0
julia> klGauss(3, 6)
18.0
julia> klGauss(1, 2)
2.0
julia> klGauss(2, 1)  # And this KL is symmetric
2.0
julia> klGauss(4, 2)
8.0
julia> klGauss(6, 8)
8.0
```

- x, y can be negative:

```julia
julia> klGauss(-3, 2)
50.0
julia> klGauss(3, -2)
50.0
julia> klGauss(-3, -2)
2.0
julia> klGauss(3, 2)
2.0
```

- With other values for `sig2x`:

```julia
julia> klGauss(3, 3, sig2x=10)
0.0
julia> klGauss(3, 6, sig2x=10)
0.45
julia> klGauss(1, 2, sig2x=10)
0.05
julia> klGauss(2, 1, sig2x=10)  # And this KL is symmetric
0.05
julia> klGauss(4, 2, sig2x=10)
0.2
julia> klGauss(6, 8, sig2x=10)
0.2
```

- With different values for `sig2x` and `sig2y`:

```julia
julia> klGauss(0, 0, sig2x=0.25, sig2y=0.5)
-0.0284...
julia> klGauss(0, 0, sig2x=0.25, sig2y=1.0)
0.2243...
julia> klGauss(0, 0, sig2x=0.5, sig2y=0.25)  # not symmetric here!
1.1534...

julia> klGauss(0, 1, sig2x=0.25, sig2y=0.5)
0.9715...
julia> klGauss(0, 1, sig2x=0.25, sig2y=1.0)
0.7243...
julia> klGauss(0, 1, sig2x=0.5, sig2y=0.25)  # not symmetric here!
3.1534...

julia> klGauss(1, 0, sig2x=0.25, sig2y=0.5)
0.9715...
julia> klGauss(1, 0, sig2x=0.25, sig2y=1.0)
0.7243...
julia> klGauss(1, 0, sig2x=0.5, sig2y=0.25)  # not symmetric here!
3.1534...
```

- **Warning:** Using [`Policies.klUCB`](https://smpybandits.github.io/docs/Policies.klUCB.html) (and variants) with `klGauss` is equivalent to use [`Policies.UCB`](https://smpybandits.github.io/docs/Policies.UCB.html), so prefer the simpler version. (**this is only for the Python version**)
"""
function klGauss(x, y, sig2x=0.25, sig2y=0.25)
    if - eps < (sig2y - sig2x) < eps
        return (x - y) ^ 2 / (2.0 * sig2x)
    else
        return (x - y) ^ 2 / (2.0 * sig2y) + 0.5 * ((sig2x/sig2y)^2 - 1 - log(sig2x/sig2y))
    end
end

# --- generic KL divergence functions for Distributions
# See https://github.com/Naereen/KullbackLeibler.jl/issues/2

using Distributions

doc"""
    function KL( D1, D2 )

Common function to compute the Kullback-Leibler for some continuous and discrete 1D parametric distributions (of module Distributions).

See [doc for continuous distributions](https://juliastats.github.io/Distributions.jl/latest/univariate.html#Continuous-Distributions-1) and [doc for discrete distributions](https://juliastats.github.io/Distributions.jl/latest/univariate.html#Discrete-Distributions-1).

Currently supported: `Bernoulli`, `Binomial`, `Poisson`, `NegativeBinomial` (discrete) and `Exponential`, `Gamma`, `Gaussian` (continuous).

Examples:

- Bernoulli:
```julia
julia> Bern1 = Distributions.Bernoulli(0.33)
Distributions.Bernoulli{Float64}(p=0.33)
julia> Bern2 = Distributions.Bernoulli(0.42)
Distributions.Bernoulli{Float64}(p=0.42)
julia> ex_kl_1 = KL( Bern1, Bern2 )  # Calc KL divergence for Bernoulli R.V
0.017063...
julia> klBern(0.33, 0.42)  # same!
0.017063...
```

- Binomial:
```julia
julia> Bin1 = Distributions.Binomial(13, 0.33)
Distributions.Binomial{Float64}(n=13, p=0.33)
julia> Bin2 = Distributions.Binomial(13, 0.42)   # must have same parameter n
Distributions.Binomial{Float64}(n=13, p=0.42)
julia> ex_kl_2 = KL( Bin1, Bin2 )  # Calc KL divergence for Binomial R.V
0.221828...
```

- Poisson:
```julia
julia> Pois1 = Distributions.Poisson(0.33)
Distributions.Poisson{Float64}(λ=0.33)
julia> Pois2 = Distributions.Poisson(0.92)
Distributions.Poisson{Float64}(λ=0.92)
julia> ex_kl_3 = KL( Pois1, Pois2 )  # Calc KL divergence for Poisson R.V
0.251657
```

- Exponential:
```julia
julia> Exp1 = Distributions.Exponential( 0.33 )
Distributions.Exponential{Float64}(θ=0.33)
julia> Exp2 = Distributions.Exponential( 0.42 )
Distributions.Exponential{Float64}(θ=0.42)
julia> ex_kl_4 = KL( Exp1, Exp2 )  # Calc KL div for Exponential R.V
0.0268763
```

- Gamma:
```julia
julia> Gam1 = Distributions.Gamma(0.5, 0.33)
Distributions.Gamma{Float64}(α=0.5, θ=0.33)
julia> Gam2 = Distributions.Gamma(0.5, 0.42)     # must have same parameter α
Distributions.Gamma{Float64}(α=0.5, θ=0.42)
julia> ex_kl_5 = KL( Gam1, Gam2 )  # Calc KL divergence for Gamma R.V
0.0134381
```

- NegativeBinomial:
```julia
julia> NegBin1 = Distributions.NegativeBinomial(40, 0.33)
Distributions.NegativeBinomial{Float64}(r=40, p=0.33)
julia> NegBin2 = Distributions.NegativeBinomial(40, 0.99)  # must have same parameter r
Distributions.NegativeBinomial{Float64}(r=40, p=0.99)
julia> ex_kl_6 = KL( NegBin1, NegBin2 )       # Calc KL divergence for NegativeBinomial R.V
17.277824
```

- Gaussian:
```julia
julia> Gauss1 = Distributions.Gaussian(0.33, 1)
Distributions.Normal{Float64}(μ=0.33, σ=1.0)
julia> Gauss2 = Distributions.Gaussian(0.42, 1)      # can have same σ
Distributions.Normal{Float64}(μ=0.42, σ=1.0)
julia> ex_kl_7 = KL( Gauss1, Gauss2 )  # Calc KL divergence for Gaussian R.V
0.0040499...

julia> Gauss3 = Distributions.Gaussian(0.42, 10)     # but can also have a different σ
Distributions.Normal{Float64}(μ=0.33, σ=10.0)
julia> ex_kl_7 = KL( Gauss1, Gauss3 )  # Calc KL divergence for Gaussian R.V
0.656697...
```
"""
function KL(D1, D2)
    throw(ErrorException("Unsupported distributions D1 or D2. Only Bernoulli, Binomial, Poisson, Exponential, Gamma, NegativeBinomial, Gaussian are supported."))
end

KL( D1::Distributions.Bernoulli, D2::Distributions.Bernoulli ) = klBern( D1.p, D2.p )

function KL( D1::Distributions.Binomial, D2::Distributions.Binomial )
    assert(D1.n == D2.n)
    return klBin( D1.p, D2.p, D1.n )
end

KL( D1::Distributions.Poisson, D2::Distributions.Poisson ) = klPoisson( D1.λ, D2.λ )

KL( D1::Distributions.Exponential, D2::Distributions.Exponential ) = klExp( D1.θ, D2.θ )

function KL( D1::Distributions.Gamma, D2::Distributions.Gamma )
    assert(D1.α == D2.α)
    return klGamma( D1.θ, D2.θ, D1.α )
end

function KL( D1::Distributions.NegativeBinomial, D2::Distributions.NegativeBinomial )
    assert(D1.r == D2.r)
    return klGamma( D1.p, D2.p, D1.r )
end

KL( D1::Distributions.Gaussian, D2::Distributions.Gaussian ) = klGauss( D1.μ, D2.μ, D1.σ, D2.σ )

# --- klucb functions, for the kl-UCB policy

doc"""
    function klucb(x, d, kl, upperbound, lowerbound=-Inf, precision=1e-6, max_iterations=50)

The generic kl-UCB index computation.

- x: value of the cum reward,
- d: upper bound on the divergence,
- kl: the KL divergence to be used (`klBern`, `klGauss`, etc),
- upperbound, lowerbound=-Inf: the known bound of the values x,
- precision=1e-6: the threshold from where to stop the research,
- max_iterations: max number of iterations of the loop (safer to bound it to reduce time complexity).

- **Note**: It uses a **bisection search**, and one call to ``kl`` for each step of the bisection search.

For example, for `klucbBern`, the two steps are to first compute an upperbound (as precise as possible) and the compute the kl-UCB index:

```julia
julia> x, d = 0.9, 0.2   # mean x, exploration term d
julia> upperbound = min(1.0, klucbGauss(x, d, sig2x=0.25))  # variance 1/4 for [0,1] bounded distributions
julia> upperbound
1.0
julia> klucb(x, d, klBern, upperbound, lowerbound=0, precision=1e-3, max_iterations=10)
0.9941...
julia> klucb(x, d, klBern, upperbound, lowerbound=0, precision=1e-6, max_iterations=10)
0.994482...
julia> klucb(x, d, klBern, upperbound, lowerbound=0, precision=1e-3, max_iterations=50)
0.9941...
julia> klucb(x, d, klBern, upperbound, lowerbound=0, precision=1e-6, max_iterations=100)  # more and more precise!
0.994489...
```

- **Note**: See below for more examples for different KL divergence functions.
"""
function klucb(x, d, kl, upperbound, lowerbound=-Inf, precision=1e-6, max_iterations=50)
    value = max(x, lowerbound)
    u = upperbound
    _count_iteration = 0
    while (_count_iteration < max_iterations) && (u - value > precision)
        _count_iteration += 1
        m = (value + u) / 2.0
        if kl(x, m) > d
            u = m
        else
            value = m
        end
    end
    return (value + u) / 2.0
end

doc"""
    function klucbBern(x, d, precision=1e-6)

kl-UCB index computation for Bernoulli distributions, using `klucb`.

- Influence of x:

```julia
julia> klucbBern(0.1, 0.2)
0.378391...
julia> klucbBern(0.5, 0.2)
0.787088...
julia> klucbBern(0.9, 0.2)
0.994489...
```

- Influence of d:

```julia
julia> klucbBern(0.1, 0.4)
0.519475...
julia> klucbBern(0.1, 0.9)
0.734714...

julia> klucbBern(0.5, 0.4)
0.871035...
julia> klucbBern(0.5, 0.9)
0.956809...

julia> klucbBern(0.9, 0.4)
0.999285...
julia> klucbBern(0.9, 0.9)
0.999995...
```
"""
function klucbBern(x, d, precision=1e-6)
    upperbound = min(1.0, klucbGauss(x, d, sig2x=0.25))  # variance 1/4 for [0,1] bounded distributions
    # upperbound = min(1.0, klucbPoisson(x, d))  # also safe, and better ?
    return klucb(x, d, klBern, upperbound, precision)
end

doc"""
    function klucbGauss(x, d, sig2x=0.25, precision=0.0)

kl-UCB index computation for Gaussian distributions.

- **Note**: it does not require any search.

- **Warning**: it works only if the good variance constant is given.

- Influence of x:

```julia
julia> klucbGauss(0.1, 0.2)
0.416227...
julia> klucbGauss(0.5, 0.2)
0.816227...
julia> klucbGauss(0.9, 0.2)
1.216227...

- Influence of d:

```julia
julia> klucbGauss(0.1, 0.4)
0.547213...
julia> klucbGauss(0.1, 0.9)
0.770820...

julia> klucbGauss(0.5, 0.4)
0.947213...
julia> klucbGauss(0.5, 0.9)
1.170820...

julia> klucbGauss(0.9, 0.4)
1.347213...
julia> klucbGauss(0.9, 0.9)
1.570820...
```

- **Warning**: Using [`Policies.klUCB`](https://smpybandits.github.io/docs/Policies.klUCB.html) (and variants) with `klucbGauss` is equivalent to use [`Policies.UCB`](https://smpybandits.github.io/docs/Policies.UCB.html), so prefer the simpler version (**this is only for the Python version**).
"""
function klucbGauss(x, d, sig2x=0.25, precision=0.0)
    return x + sqrt(2 * sig2x * d)
end

doc"""
    function klucbPoisson(x, d, precision=1e-6)

kl-UCB index computation for Poisson distributions, using `klucb`.

- Influence of x:

```julia
julia> klucbPoisson(0.1, 0.2)
0.450523...
julia> klucbPoisson(0.5, 0.2)
1.089376...
julia> klucbPoisson(0.9, 0.2)
1.640112...
```

- Influence of d:

```julia
julia> klucbPoisson(0.1, 0.4)
0.693684...
julia> klucbPoisson(0.1, 0.9)
1.252796...

julia> klucbPoisson(0.5, 0.4)
1.422933...
julia> klucbPoisson(0.5, 0.9)
2.122985...

julia> klucbPoisson(0.9, 0.4)
2.033691...
julia> klucbPoisson(0.9, 0.9)
2.831573...
```
"""
function klucbPoisson(x, d, precision=1e-6)
    upperbound = x + d + sqrt(d * d + 2 * x * d)  # looks safe, to check: left (Gaussian) tail of Poisson dev
    return klucb(x, d, klPoisson, upperbound, precision)
end

doc"""
    function klucbExp(x, d, precision=1e-6)

kl-UCB index computation for exponential distributions, using `klucb`.

- Influence of x:

```julia
julia> klucbExp(0.1, 0.2)
0.202741...
julia> klucbExp(0.5, 0.2)
1.013706...
julia> klucbExp(0.9, 0.2)
1.824671...
```

- Influence of d:

```julia
julia> klucbExp(0.1, 0.4)
0.285792...
julia> klucbExp(0.1, 0.9)
0.559088...

julia> klucbExp(0.5, 0.4)
1.428962...
julia> klucbExp(0.5, 0.9)
2.795442...

julia> klucbExp(0.9, 0.4)
2.572132...
julia> klucbExp(0.9, 0.9)
5.031795...
```
"""
function klucbExp(x, d, precision=1e-6)
    if d < 0.77
        # XXX where does this value come from?
        upperbound = x / (1 + 2.0 / 3 * d - sqrt(4.0 / 9 * d * d + 2 * d))
        # safe, klexp(x,y) >= e^2/(2*(1-2e/3)) if x=y(1-e)
    else
        upperbound = x * exp(d + 1)
    end
    if d > 1.61
        # XXX where does this value come from?
        lowerbound = x * exp(d)
    else
        lowerbound = x / (1 + d - sqrt(d * d + 2 * d))
    end
    return klucb(x, d, klGamma, upperbound, lowerbound, precision)
end

doc"""
    function klucbGamma(x, d, precision=1e-6)

kl-UCB index computation for Gamma distributions, using `klucb`.

- Influence of x:

```julia
julia> klucbGamma(0.1, 0.2)
0.202...
julia> klucbGamma(0.5, 0.2)
1.013...
julia> klucbGamma(0.9, 0.2)
1.824...
```

- Influence of d:

```julia
julia> klucbGamma(0.1, 0.4)
0.285...
julia> klucbGamma(0.1, 0.9)
0.559...

julia> klucbGamma(0.5, 0.4)
1.428...
julia> klucbGamma(0.5, 0.9)
2.795...

julia> klucbGamma(0.9, 0.4)
2.572...
julia> klucbGamma(0.9, 0.9)
5.031...
```
"""# FIXME this one is wrong!
function klucbGamma(x, d, precision=1e-6)
    if d < 0.77
        # XXX where does this value come from?
        upperbound = x / (1 + 2.0 / 3 * d - sqrt(4.0 / 9 * d * d + 2 * d))
        # safe, klexp(x,y) >= e^2/(2*(1-2e/3)) if x=y(1-e)
    else
        upperbound = x * exp(d + 1)
    end
    if d > 1.61
        # XXX where does this value come from
        lowerbound = x * exp(d)
    else
        lowerbound = x / (1 + d - sqrt(d * d + 2 * d))
    end
    # FIXME specify the value for a !
    return klucb(x, d, klGamma, max(upperbound, 1e2), min(-1e2, lowerbound), precision)
end

end

# --- Debugging

if "test" in ARGS
    # Code for debugging purposes.
    println("FIXME write these tests!")
    println("\nDone for tests of 'kullback_leibler.jl' ...")
end


text:
function _lobatto_legendre_fast(s, T)
    c, b = FastGaussQuadrature.gausslobatto(s)
    shift!(b,c)
    QuadratureRule(2s-2, c, b, T)
end

"""
Lobatto-Legendre quadrature.
"""
function LobattoLegendreQuadrature(::Type{T}, s::Integer; IT=BigFloat, fast=false) where {T}
    if s == 1
        throw(ErrorException("Lobatto quadrature is not defined for one stage."))
    end

    if fast
        return _lobatto_legendre_fast(s, T)
    end

    D(k) = Polynomials.derivative(Polynomial(IT[0, 1, -1])^(k-1), k-2)
    P(k,x) = Polynomials.derivative(Polynomial(IT[-1, 0, 1])^k, k)(x) / factorial(k) / 2^k
    c = sort(IT.(Polynomials.roots(D(s)))); c[begin] = 0; c[end] = 1;
    b = [ 1 / ( s*(s-1) * P(s-1, 2c[i] - 1)^2 ) for i in 1:s ]
    return QuadratureRule(2s-2, c, b, T)
end

LobattoLegendreQuadrature(s; kwargs...) = LobattoLegendreQuadrature(Float64, s; kwargs...)


text:
"""
    function vader_koala

Another canonical koala?

# Example

```jldoctest
julia> cowsay("Luke, you are my joey!", cow=Cowsay.vader_koala)
 ________________________
< Luke, you are my joey! >
 ------------------------
   \\
    \\        .
     .---.  //
    Y|o o|Y//
   /_(i=i)K/
   ~()~*~()~
    (_)-(_)

     Darth
     Vader
     koala
```
"""
function vader_koala(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
   $thoughts
    $thoughts        .
     .---.  //
    Y|o o|Y//
   /_(i=i)K/
   ~()~*~()~
    (_)-(_)

     Darth
     Vader
     koala
"""
return the_cow
end


text:
# This file investigates mean square error (MSE) with respect to free variables in hidden 1D interpolation

using FractalTools
using Makie

# Construct interpolation data
f(x) = sin(2π * x)
g(x) = cos(2π * x)
dx   = 0.1
xi   = 0.
xf   = 1.
x    = collect(xi : dx : xf)
y    = f.(x)
z    = g.(x)
pts  = collect.(zip(x, y, z))

# Construct test data
xt   = collect(xi : 0.1dx : xf)
npts = length(xt)

# Compute errors
fvals = f.(xt)
freevars = 0.001 : 0.01 : 0.999
mse = map(freevars) do freevar
    interp = interpolate(pts, HInterp1D(fill(freevar, 2, 2)))
    ivals = getindex.(interp.(xt), 1)
    sum((fvals - ivals).^2) / npts
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Free Variable", ylabel="MSE", title="1D Hidden Interpolation MSE")
stem!(ax, freevars, mse, color=:black)
save(joinpath(@__DIR__, "hinterp1d_mse.png"), fig)
display(fig)


text:
const named_array_plot_types = [BarPlot, Heatmap, Volume]

struct NamedSparseArray{S<:Tuple, V<:AbstractVector}
    keys::S
    values::V
    function NamedSparseArray(args::AbstractVector...)
        keys = Base.front(args)
        values = last(args)
        return new{typeof(keys), typeof(values)}(keys, values)
    end
end

# this logic may be moved to AbstractPlotting
function dense(kv::NamedSparseArray)
    keys, values = kv.keys, kv.values
    labels = map(collect∘uniquesorted, keys)
    indices = map(Base.OneTo∘length, labels)
    converter = map((k, v) -> Dict(zip(k, v)), labels, indices)
    d = zeros(eltype(values), indices)
    for (k, v) in zip(StructArray(keys), values)
        I = map(getindex, converter, k)
        d[I...] = v
    end
    return labels, d
end

function convert_arguments(P::PlotFunc, m::NamedSparseArray)
    labels, values = dense(m)
    ptype = plottype(P, named_array_plot_types[length(labels)])
    to_plotspec(ptype, convert_arguments(ptype, labels..., values))
end

# TODO optimize with refarray and refvalue for pooled data
function _frequency(args...)
    s = StructArray(map(vec, args))
    gp = GroupPerm(s)
    sp = sortperm(gp)
    itr = (s[sp[first(range)]] => length(range) for range in gp)
    keys, values = fieldarrays(StructArray(itr, unwrap = t -> t <: Tuple))
    return NamedSparseArray(fieldarrays(keys)..., values)
end

const frequency = Analysis(_frequency)


text:
# generate examples
import Literate

EXAMPLEDIR = joinpath(@__DIR__, "src", "literate")
GENERATEDDIR = joinpath(@__DIR__, "src", "examples")
mkpath(GENERATEDDIR)

# Copy supplementary files first
suplementary_fileextensions = [".inp", ".svg", ".png", ".jpg", ".gif"]
for example in readdir(EXAMPLEDIR)
    if any(endswith.(example, suplementary_fileextensions))
        cp(joinpath(EXAMPLEDIR, example), joinpath(GENERATEDDIR, example); force=true)
    end
end

for example in readdir(EXAMPLEDIR)
    if endswith(example, ".jl")
        input = abspath(joinpath(EXAMPLEDIR, example))
        script = Literate.script(input, GENERATEDDIR)
        code = strip(read(script, String))

        # remove "hidden" lines which are not shown in the markdown
        line_ending_symbol = occursin(code, "\r\n") ? "\r\n" : "\n"
        code_clean = join(filter(x->!endswith(x,"#hide"),split(code, r"\n|\r\n")), line_ending_symbol)

        mdpost(str) = replace(str, "@__CODE__" => code_clean)
        Literate.markdown(input, GENERATEDDIR, postprocess = mdpost)
        Literate.notebook(input, GENERATEDDIR, execute = is_ci) # Don't execute locally
    else
        if !any(endswith.(example, suplementary_fileextensions))
            @warn "ignoring $example"
        end
    end
end

# remove any .vtu files in the generated dir (should not be deployed)
cd(GENERATEDDIR) do
    foreach(file -> endswith(file, ".vtu") && rm(file), readdir())
    foreach(file -> endswith(file, ".pvd") && rm(file), readdir())
end


text:
module Storage

include("module_common.jl")

const _libcloud_storage_types = PyCall.PyNULL()
const _libcloud_storage_providers = PyCall.PyNULL()
const _libcloud_storage_base = PyCall.PyNULL()

immutable StorageDriver
    o::PyObject

    function StorageDriver(provider::Compat.String, args...; kwargs...)
        cls = _libcloud_storage_providers[:get_driver](provider)
        StorageDriver(cls(args...; kwargs...))
    end
    StorageDriver(driver::PyObject) = new(driver)
end
show(io::IO, c::StorageDriver) = print(io, c.o[:__str__]())

immutable Container
    o::PyObject

    name::Compat.String
    extra::Dict
    driver::StorageDriver

    function Container(o::PyObject)
        new(o, o[:name], o[:extra], StorageDriver(o[:driver]))
    end
end
PyObject(c::Container) = c.o
convert(::Type{Container}, o::PyObject) =  Container(o)
show(io::IO, o::Container) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_storage_base, "Container") Container

immutable Object
    o::PyObject

    name::Compat.String
    size::Int
    hash::Nullable{Compat.String}
    container::Container
    extra::Dict
    meta_data::Dict
    driver::StorageDriver

    function Object(o::PyObject)
        sz = o[:size]
        szint = isa(sz, Integer) ? Int(sz) : parse(Int, sz)
        new(o, o[:name], szint, o[:hash], Container(o[:container]), o[:extra], o[:meta_data] , StorageDriver(o[:driver]))
    end
end
PyObject(o::Object) = o.o
convert(::Type{Object}, o::PyObject) =  Object(o)
show(io::IO, o::Object) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_storage_base, "Object") Object

const _storagedriver_fns = [
    :create_container, :delete_container, :delete_object, :download_object, :download_object_as_stream,
    :enable_container_cdn, :enable_object_cdn, :get_container, :get_container_cdn_url, :get_object,
    :get_object_cdn_url, :iterate_container_objects, :iterate_containers, :list_container_objects,
    :list_containers, :upload_object, :upload_object_via_stream
]

const _container_fns = [
    :iterate_objects, :list_objects, :get_cdn_url, :enable_cdn, :get_object, :upload_object,
    :upload_object_via_stream, :download_object, :download_object_as_stream, :delete_object, :delete
]

const _object_fns = [
    :get_cdn_url, :enable_cdn, :download, :as_stream, :delete
]

for f in union(Set(_storagedriver_fns), Set(_base_driver_fns))
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_storage_base, "StorageDriver", $sf) $(f)(storage::StorageDriver, args...; kwargs...) = storage.o[$(sf)](args..., kwargs...)
end

for f in _container_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_storage_base, "Container", $sf) $(f)(cont::Container, args...; kwargs...) = cont.o[$(sf)](args..., kwargs...)
end

for f in _object_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_storage_base, "Object", $sf) $(f)(obj::Object, args...; kwargs...) = obj.o[$(sf)](args..., kwargs...)
end

# Initialize cloud storage
function __init__()
    copy!(_libcloud_storage_types, pyimport("libcloud.storage.types"))
    copy!(_libcloud_storage_providers, pyimport("libcloud.storage.providers"))
    copy!(_libcloud_storage_base, pyimport("libcloud.storage.base"))

    global const StorageProvider = pywrap(_libcloud_storage_types[:Provider])

    _map_types(_libcloud_storage_base, (Container, Object))
end

# types
export StorageProvider, StorageDriver, Container, Object
# storage driver functions
export create_container, delete_container, delete_object, download_object, download_object_as_stream,
       enable_container_cdn, enable_object_cdn, get_container, get_container_cdn_url, get_object,
       get_object_cdn_url, iterate_container_objects, iterate_containers, list_container_objects,
       list_containers, upload_object, upload_object_via_stream
# container functions
export iterate_objects, list_objects, get_cdn_url, enable_cdn, get_object, upload_object,
       upload_object_via_stream, download_object, download_object_as_stream, delete_object, delete
# object functions
export get_cdn_url, enable_cdn, download, as_stream, delete

end # module Storage


text:
@enum TerritoryStatus begin
    black_stone = 2
    black_territory = 1
    dame_territory = 0
    white_territory = -1
    white_stone = -2
end

struct Territory
    num_black_territory::Int
    num_white_territory::Int
    num_black_stones::Int
    num_white_stones::Int
    num_dame_territory::Int
    dame_points::Set{Point}

    function Territory(territorymap::Dict{Point, TerritoryStatus})
        num_black_territory = 0
        num_white_territory = 0
        num_black_stones = 0
        num_white_stones = 0
        num_dame_territory = 0
        dame_points = Set()

        for (point, status) in territorymap
            if status === black_stone
                num_black_stones += 1
            elseif status === white_stone
                num_white_stones += 1
            elseif status === black_territory
                num_black_territory += 1
            elseif status == white_territory
                num_white_territory += 1
            elseif status == dame_territory
                num_dame_territory += 1
                push!(dame_points, point)
            end
        end

        new(num_black_territory, num_white_territory, num_black_stones,
            num_white_stones, num_dame_territory, dame_points)
    end
end

struct GameResult
    b::Int
    w::Int
    komi::Float64
end

score(result::GameResult) = result.b  - result.w - result.komi
winnercolor(result::GameResult) = score(result) > 0 ? black : white

function Base.show(io::IO, result::GameResult)
    s = score(result)
    if s > 0
        print(io, "B+$s")
    else
        print(io, "W+$(-s)")
    end
end

function evaluate_territory(board::Board)
    territorymap = Dict{Point, TerritoryStatus}()
    for r in 1:board.num_rows
        for c in 1:board.num_cols
            point = (row=r, col=c)

            if point in keys(territorymap) # <1>
                continue
            end

            string = board[point]
            if !isvoid(string) # <2>
                territorymap[point] = (string.color === black) ? black_stone : white_stone
            else
                group, nbrs = collect_region(point, board)
                if length(nbrs) == 1 # <3>
                    nbr_string = pop!(nbrs)
                    status = (nbr_string.color === black) ? black_territory : white_territory
                else
                    status = dame_territory # <4>
                end
                for q in group
                    territorymap[q] = status
                end
            end
        end
    end
    Territory(territorymap)
end

# <1> Skip the point, if you already visited this as part of a different group.
# <2> If the point is a stone, add it as status.
# <3> If a point is completely surrounded by black or white stones, count it as territory.
# <4> Otherwise the point has to be a neutral point, so we add it to dame.
# end::scoring_evaluate_territory[]

function collect_region(start_pos::Point, board::Board, visited=Set{Point}())
    if start_pos in visited
        return Point[], Set{Point}()
    end
    all_points = [start_pos]
    all_borders = Set{GoString}()
    push!(visited, start_pos)
    here = board[start_pos]
    for next_p in nbrs(start_pos)
        !is_on_grid(board, next_p) && continue
        nbr_string = board[next_p]
        if nbr_string === here
            points, borders = collect_region(next_p, board, visited)
            append!(all_points, points)
            union!(all_borders, borders)
        else
            push!(all_borders, nbr_string)
        end
    end
    return all_points, all_borders
end

function compute_game_result(game_state::GameState, komi::Float64 = 7.5)
    territory = evaluate_territory(game_state.board)
    return GameResult(
        territory.num_black_territory + territory.num_black_stones,
        territory.num_white_territory + territory.num_white_stones,
        komi)
end


text:
module SIAN

using Nemo
using StructuralIdentifiability: PreprocessODE, eval_at_nemo
using LinearAlgebra
using Singular
using GroebnerBasis
using MacroTools
using OrderedCollections
using ModelingToolkit
using Logging # * for the @warn macro
include("util.jl")
include("ODE.jl")

# exports are required for "using"
export identifiability_ode, PreprocessODE
export @ODEmodel
export ODE
export Nemo, OrderedDict, Generic, macroexpand, macrohelper_extract_vars, macrohelper_clean,  fmpq_mpoly, get_parameters
# ------------------------------------------------------------------------------
"""
    func make_derivative(var_name, der_order)

Given a variable name `var_name` add a derivative order `der_order`.
"""
function make_derivative(var_name, der_order)
    return(string(var_name, "_", der_order))
end

# ------------------------------------------------------------------------------
"""
    func add_to_var(vr, ring, r)

Convert a variable `vr` to a derivative of order `r` and convert the result to symbol.
"""
function add_to_var(vr, ring, r)
    return str_to_var(make_derivative(vr, r), ring)
end

# ------------------------------------------------------------------------------
"""
    func create_jet_ring(var_list, param_list, max_ord)

Given a list of variables `var_list` and a list of parameters `param_list`, create a jet ring of derivatives up to order `max_ord`.
"""
function create_jet_ring(var_list, param_list, max_ord)
    varnames = vcat(vec(["$(s)_$i" for s in var_list, i in 0:max_ord]), "z_aux", ["$(s)_0" for s in param_list])
    return Nemo.PolynomialRing(Nemo.QQ, varnames)[1]
end

# ------------------------------------------------------------------------------
"""
    func differentiate_all(diff_poly, var_list, shft, max_ord)

Differentiate a polynomial `diff_poly` with respect to `var_list` up to `max_ord` order.
"""
function differentiate_all(diff_poly, var_list, shft, max_ord)
    result = 0
    for i in 1:(shft * (max_ord + 1))
        result = result + derivative(diff_poly, var_list[i]) * var_list[i + shft]
    end
    return(result)
end

# ------------------------------------------------------------------------------
"""
    func sample_point(bound, x_vars, y_vars, u_variables, all_params, X_eq, Y_eq, Q)

Sample random values for parameters of the polynomial system.
"""
function sample_point(bound, x_vars, y_vars, u_variables, all_params, X_eq, Y_eq, Q)
    local u_hat, theta_hat, all_subs

    s = length(all_params)
    y_hat_vars = Array{fmpq_mpoly}(undef, 0)
    y_hat_vals = Array{fmpq}(undef, 0)

    while true
        theta_hat = [fmpq(rnd) for rnd in rand(0:bound, s)]
        u_hat =  [fmpq(rnd) for rnd in rand(0:bound, length(u_variables))]
        all_subs = [vcat(all_params, u_variables), vcat(theta_hat, u_hat)]
        if evaluate(Q, all_subs[1], all_subs[2]) == 0
            next
        end
        vls = insert_zeros_to_vals(all_subs[1], all_subs[2])
        for i in 0:(s + 1)
            for j in 1:length(y_vars)
                eq = Y_eq[(j - 1) * (s + 2) + i + 1][2]
                vl = evaluate(unpack_fraction(eq)[1], vls) // evaluate(unpack_fraction(eq)[2], vls)
                y_hat_vars = vcat(y_hat_vars, Y_eq[(j - 1) * (s + 2) + i + 1][1])
                y_hat_vals = vcat(y_hat_vals, vl)
                vls[var_index(Y_eq[(j - 1) * (s + 2) + i + 1][1])] = vl
                all_subs = [vcat(all_subs[1], Y_eq[(j - 1) * (s + 2) + i + 1][1]), vcat(all_subs[2], vl)]
            end
            for j in 1:length(x_vars)
                eq = X_eq[(j - 1) * (s + 2) + i + 1][2]
                vl = evaluate(unpack_fraction(eq)[1], vls) // evaluate(unpack_fraction(eq)[2], vls)
                all_subs = [vcat(all_subs[1], X_eq[(j - 1) * (s + 2) + i + 1][1]),vcat(all_subs[2], vl)]
                vls[var_index(X_eq[(j - 1) * (s + 2) + i + 1][1])] = vl
            end
        end
        break
    end
    return([[y_hat_vars, y_hat_vals], [u_variables, u_hat], [all_params, theta_hat], all_subs])
end

# ------------------------------------------------------------------------------
"""
    func jacobi_matrix(pol_arr, vrs, vals)

Generate a Jacobi matrix from a given array of polynomial `pol_arr`,
with respect to variables `vars`.
The matrix is evaluated at `vals` from a symbolic to numeric representation.
"""
function jacobi_matrix(pol_arr, vrs, vals)
    m = Nemo.MatrixSpace(Nemo.QQ, length(pol_arr), length(vrs))()
    for i in 1:length(pol_arr)
        for j in 1:length(vrs)
            m[i, j] = evaluate(derivative(pol_arr[i], vrs[j]), vals)
        end
    end
    return(m)
end

# ------------------------------------------------------------------------------
"""
    func get_order_var2(diff_var, non_jet_vars, shft, s)
"""
function get_order_var2(diff_var, non_jet_vars, shft, s)
    idx = var_index(diff_var)
    if idx <= shft * (s + 3)
        return([non_jet_vars[rem(idx - 1, shft) + 1], div(idx - 1, shft)])
    else
        return([non_jet_vars[idx - shft * (s + 2) - 1], 0])
    end
end

# ------------------------------------------------------------------------------
"""
    func get_order_var(diff_var, non_jet_ring)

A helper function to obtain derivative order from string.
"""
function get_order_var(diff_var, non_jet_ring)
    rex = match(r"^(.*_)([0-9]+)$", string(diff_var))
    if rex === nothing
        return(["", ""])
    else
        return([str_to_var(first(rex[1], length(rex[1]) - 1), non_jet_ring), parse(Int, rex[2])])
    end
end

# ------------------------------------------------------------------------------
"""
    func get_vars(diff_poly, var_list, non_jet_vars, shft, s)

Get variables from `diff_poly` based on the intersection with `var_list`.
"""
function get_vars(diff_poly, var_list, non_jet_vars, shft, s)
    return [v for v in vars(diff_poly) if get_order_var2(v, non_jet_vars, shft, s)[1] in var_list]
end

# ------------------------------------------------------------------------------
"""
    func compare_diff_var(dvl, dvr, non_jet_vars, shft, s)

Comparison method of variables based on order.
"""
function compare_diff_var(dvl, dvr, non_jet_vars, shft, s)
    vl, hl = get_order_var2(dvl, non_jet_vars, shft, s)
    vr, hr = get_order_var2(dvr, non_jet_vars, shft, s)
    if hl != hr
        return (hl > hr)
    end
    if length(string(vl)) != length(string(vr))
        return (length(string(vl)) > length(string(vr)))
    end
    return (vr >= vl)
end

# ------------------------------------------------------------------------------
"""
    func get_parameters(ode; initial_conditions=true)

Retrieve parameters from the `ode` system. Retrieve initial conditions if `initial_conditions` is set `true`.

## Input
    - `ode::ODE` - an ODE system
    - `initial_conditions::Bool` - whether to extract initial conditions. Default `true`.

## Output
        - Array of parameters (and initial conditions).
"""
function get_parameters(ode; initial_conditions=true)
    if initial_conditions
        return vcat(ode.parameters, ode.x_vars)
    else
        return ode.parameters
    end
end

# ------------------------------------------------------------------------------
# Main Code
# ------------------------------------------------------------------------------
#

function identifiability_ode(ode::ModelingToolkit.ODESystem, params_to_assess=[]; p=0.99, p_mod=0, nthrds=1)
    ode_prep, input_syms, gens_ = PreprocessODE(ode)
    t = ModelingToolkit.arguments(ModelingToolkit.states(ode)[1])[1]
    if length(params_to_assess)==0
        params_to_assess_ = SIAN.get_parameters(ode_prep)
        nemo2mtk = Dict(gens_.=>input_syms)
    else
        params_to_assess_ = [eval_at_nemo(each, Dict(syms.=>gens_)) for each in params_to_assess]
        nemo2mtk = Dict(params_to_assess_.=>params_to_assess)
    end

    res = identifiability_ode(ode_prep, params_to_assess_; p=p, p_mod=p_mod, nthrds=1)

    @info "Post-Processing: Converting Nemo output to ModelingToolkit types"
    out = Dict()
    for (id_type, pars) in pairs(res)
        out[id_type] = [ModelingToolkit.Num(substitute(nemo2mtk[each], t=>0)) for each in pars]
    end
    return out
end

"""
    func identifiability_ode(ode, params_to_assess; p=0.99, p_mod=0, nthrds=1)

Perform identifiability check for a given `ode` system with respect to parameters in `params_to_assess` list.

## Input

  - `ode` - an ODE system returned by the `@ODEmodel` macro.
  - `params_to_assess` - an array of parameters returned by `get_parameters` function.
  - `p` - probability of correctness, default `0.99`.
  - `p_mod` - a prime characteristic, default `0`.
  - `nthrds` - number of threads for concurrency, default `1`.
  - `infolevel` - an integer, controls the verbosity of Groebner Basis computation, default `0` (no output).
                 See GroebnerBasis.jl documentation for details.
"""
function identifiability_ode(ode, params_to_assess; p=0.99, p_mod=0, nthrds=1, infolevel=0)

    println("Solving the problem")

    if p_mod != 0
        @warn "Using `p_mod!=0` does not guarantee the same probability of correctness but allows to run the program faster. This warning was raised by `p_mod = $p_mod`."
    end
    # 1.Construct the maximal system

    # (a) ---------------

    println("Constructing the maximal system")

    non_jet_ring = ode.poly_ring
    all_indets = gens(non_jet_ring)
    x_vars = ode.x_vars
    y_vars = ode.y_vars
    u_vars = ode.u_vars
    mu = ode.parameters

    n = length(x_vars)
    m = length(y_vars)
    u = length(u_vars)
    s = length(mu) + n

    Rjet = create_jet_ring(vcat(x_vars, y_vars, u_vars), mu, s + 2)
    gens_Rjet = gens(Rjet)
    z_aux = gens_Rjet[end - length(mu)]
    x_eqs = collect(values(ode.x_equations))
    y_eqs = collect(values(ode.y_equations))

    x_eqs = [[add_to_var(x_vars[i], Rjet, 1), add_zero_to_vars(unpack_fraction(x_eqs[i])[1], Rjet) // add_zero_to_vars(unpack_fraction(x_eqs[i])[2], Rjet)] for i in 1:n]
    y_eqs = [[add_to_var(y_vars[i], Rjet, 0), add_zero_to_vars(unpack_fraction(y_eqs[i])[1], Rjet) // add_zero_to_vars(unpack_fraction(y_eqs[i])[2], Rjet)] for i in 1:m]

    eqs = vcat(x_eqs, y_eqs)
    Q = foldl(lcm, [unpack_fraction(ex[2])[2] for ex in eqs])

    not_int_cond_params = gens_Rjet[(end - length(ode.parameters) + 1):end]
    all_params = vcat(not_int_cond_params, gens_Rjet[1:n])
    x_variables = gens_Rjet[1:n]
    for i in 1:(s + 1)
        x_variables = vcat(x_variables, gens_Rjet[(i * (n + m + u) + 1):(i * (n + m + u) + n)])
    end
    u_variables = gens_Rjet[(n + m + 1):(n + m + u)]
    for i in 1:(s + 1)
        u_variables = vcat(u_variables, gens_Rjet[((n + m + u) * i + n + m + 1):((n + m + u) * (i + 1))])
    end

    # (b,c) -------------
    X = Array{fmpq_poly}(undef, 0)
    X_eq = Array{fmpq_poly}(undef, 0)
    for i in 1:n
        X = vcat(X, [Array{fmpq_poly}(undef, 0)])
        poly_d = unpack_fraction(x_eqs[i][1] - x_eqs[i][2])[1]
        for j in 0:s + 1
            if j > 0
                poly_d = differentiate_all(poly_d, gens_Rjet, n + m + u, j)
            end
            leader = gens_Rjet[i + (n + m + u) * (j + 1)]
            separant = derivative(poly_d, leader)
            X[i] = vcat(X[i], poly_d)
            X_eq = vcat(X_eq, [[leader,-(poly_d - separant * leader) // separant]])
        end
    end

    # (d,e) -----------
    Y = Array{fmpq_poly}(undef, 0)
    Y_eq = Array{fmpq_poly}(undef, 0)
    for i in 1:m
        Y = vcat(Y, [Array{fmpq_poly}(undef, 0)])
        poly_d = unpack_fraction(y_eqs[i][1] - y_eqs[i][2])[1]
        for j in 0:s + 1
            if j > 0
                poly_d = differentiate_all(poly_d, gens_Rjet, n + m + u, j - 1)
            end
            leader = gens_Rjet[i + n + (n + m + u) * j]
            separant = derivative(poly_d, leader)
            Y[i] = vcat(Y[i], poly_d)
            Y_eq = vcat(Y_eq, [[leader,-(poly_d - separant * leader) // separant]])
        end
    end

    # 2.Truncate
    println("Truncating")

    # (a) -----------------------
    d0 = BigInt(maximum(vcat([total_degree(unpack_fraction(Q * eq[2])[1]) for eq in eqs], total_degree(Q))))

    # (b) -----------------------
    D1 = floor(BigInt, (length(params_to_assess) + 1) * 2 * d0 * s * (n + 1) * (1 + 2 * d0 * s) / (1 - p))

    # (c, d) ---------------
    sample = sample_point(D1, x_vars, y_vars, u_variables, all_params, X_eq, Y_eq, Q)
    all_subs = sample[4]
    u_hat = sample[2]
    y_hat = sample[1]

    # (e) ------------------
    alpha = [1 for i in 1:n]
    beta = [0 for i in 1:m]
    Et = Array{fmpq_poly}(undef, 0)
    x_theta_vars = all_params
    prolongation_possible = [1 for i in 1:m]

    # (f) ------------------
    all_x_theta_vars_subs = insert_zeros_to_vals(all_subs[1], all_subs[2])
    eqs_i_old = Array{fmpq_mpoly}(undef, 0)
    evl_old = Array{fmpq_mpoly}(undef, 0)
    while sum(prolongation_possible) > 0
        for i in 1:m
            if prolongation_possible[i] == 1
                eqs_i = vcat(Et, Y[i][beta[i] + 1])
                evl     = [evaluate(eq, vcat(u_hat[1], y_hat[1]), vcat(u_hat[2], y_hat[2])) for eq in eqs_i if !(eq in eqs_i_old)]
                evl_old = vcat(evl_old, evl)
                JacX    = jacobi_matrix(evl_old, x_theta_vars, all_x_theta_vars_subs)
                eqs_i_old = eqs_i
                if LinearAlgebra.rank(JacX) == length(eqs_i)
                    Et = vcat(Et, Y[i][beta[i] + 1])
                    beta[i] = beta[i] + 1
                    # adding necessary X-equations
                    polys_to_process = vcat(Et, [Y[k][beta[k] + 1] for k in 1:m])
                    while length(polys_to_process) != 0
                        new_to_process = Array{fmpq_mpoly}(undef, 0)
                        vrs = Set{fmpq_mpoly}()
                        for poly in polys_to_process
                            vrs = union(vrs, [v for v in vars(poly) if v in x_variables])
                        end
                        vars_to_add = Set{fmpq_mpoly}(v for v in vrs if !(v in x_theta_vars))
                        for v in vars_to_add
                            x_theta_vars = vcat(x_theta_vars, v)
                            ord_var = get_order_var2(v, all_indets, n + m + u, s)
                            var_idx = var_index(ord_var[1])
                            poly = X[ var_idx ][ ord_var[2] ]
                            Et = vcat(Et, poly)
                            new_to_process = vcat(new_to_process, poly)
                            alpha[ var_idx ] = max(alpha[ var_idx ], ord_var[2] + 1)
                        end
                        polys_to_process = new_to_process
                    end
                else
                    prolongation_possible[i] = 0
                end
            end
        end
    end

    println("Assessing local identifiability")

    max_rank = length(Et)
    for i in 1:m
        for j in (beta[i] + 1):length(Y[i])
            to_add = true
            for v in get_vars(Y[i][j], x_vars, all_indets, n + m + u, s)
                if !(v in x_theta_vars)
                    to_add = false
                end
            end
            if to_add
                beta[i] = beta[i] + 1
                Et = vcat(Et, Y[i][j])
            end
        end
    end

    theta_l = Array{fmpq_mpoly}(undef, 0)
    params_to_assess = [add_to_var(param, Rjet, 0) for param in params_to_assess]
    Et_eval_base = [evaluate(e, vcat(u_hat[1], y_hat[1]), vcat(u_hat[2], y_hat[2])) for e in Et]
    for param_0 in params_to_assess
        other_params = [v for v in x_theta_vars if v != param_0]
        Et_subs = [evaluate(e, [param_0], [evaluate(param_0, all_x_theta_vars_subs)]) for e in Et_eval_base]
        JacX = jacobi_matrix(Et_subs, other_params, all_x_theta_vars_subs)
        if LinearAlgebra.rank(JacX) != max_rank
            theta_l = vcat(theta_l, param_0)
        end
    end

    if length(theta_l) == 0
        println("\n=== Summary ===")
        println("Globally identifiable parameters:                 []")
        println("Locally but not globally identifiable parameters: []")
        println("Not identifiable parameters:                      [", join(params_to_assess, ", "), "]")
    else
        println("Locally identifiable parameters: [", join([get_order_var(th, non_jet_ring)[1] for th in theta_l], ", "), "]")
        println("Not identifiable parameters:     [", join([get_order_var(th, non_jet_ring)[1] for th in setdiff(params_to_assess, theta_l)], ", "), "]")

        # 3. Randomize.
        println("Randomizing")
        # (a) ------------
        deg_variety =  foldl(*, [BigInt(total_degree(e)) for e in Et])
        D2 = floor(BigInt, 6 * length(theta_l) * deg_variety * (1 + 2 * d0 * maximum(beta)) / (1 - p))
        # (b, c) ---------
        sample = sample_point(D2, x_vars, y_vars, u_variables, all_params, X_eq, Y_eq, Q)
        y_hat = sample[1]
        u_hat = sample[2]
        theta_hat = sample[3]

        # (d) ------------
        Et_hat = [evaluate(e, vcat(y_hat[1], u_hat[1]), vcat(y_hat[2], u_hat[2])) for e in Et]
        Et_x_vars = Set{fmpq_mpoly}()
        for poly in Et_hat
            Et_x_vars = union(Et_x_vars, get_vars(poly, x_vars, all_indets, n + m + u, s))
        end
        Q_hat = evaluate(Q, u_hat[1], u_hat[2])
        vrs_sorted = vcat(sort([e for e in Et_x_vars], lt=(x, y) -> compare_diff_var(x, y, all_indets, n + m + u, s)), z_aux, sort(not_int_cond_params, rev=true))

        # 4. Determine.
        println("GB computation")

        if p_mod > 0
            Et_hat = [_reduce_poly_mod_p(e, p_mod) for e in Et_hat]
            z_aux = _reduce_poly_mod_p(z_aux, p_mod)
            Q_hat = _reduce_poly_mod_p(Q_hat, p_mod)
            Rjet_new, vrs_sorted = Singular.PolynomialRing(Singular.Fp(p_mod), [string(v) for v in vrs_sorted], ordering=:degrevlex)
        else
            Rjet_new, vrs_sorted = Singular.PolynomialRing(Singular.QQ, [string(v) for v in vrs_sorted], ordering=:degrevlex)
        end

        theta_g = Array{spoly}(undef, 0)
        Et_hat = [parent_ring_change(e, Rjet_new) for e in Et_hat]
        gb = GroebnerBasis.f4(Ideal(Rjet_new, vcat(Et_hat, parent_ring_change(z_aux * Q_hat, Rjet_new) - 1)), nthrds=nthrds, infolevel=infolevel)
        println("Remainder computation")

        if p_mod > 0
            theta_l_new = [parent_ring_change(_reduce_poly_mod_p(th, p_mod), Rjet_new) for th in theta_l]

            for i in 1:length(theta_l)
                if Singular.reduce(theta_l_new[i], gb) == parent_ring_change(_reduce_poly_mod_p(Rjet(theta_hat[2][findfirst(isequal(theta_l[i]), theta_hat[1])]), p_mod), Rjet_new)
                    theta_g = vcat(theta_g, theta_l_new[i])
                end
            end
        else
            theta_l_new = [parent_ring_change(th, Rjet_new) for th in theta_l]

            for i in 1:length(theta_l)
                if Singular.reduce(theta_l_new[i], gb) == parent_ring_change(Rjet(theta_hat[2][findfirst(isequal(theta_l[i]), theta_hat[1])]), Rjet_new)
                    theta_g = vcat(theta_g, theta_l_new[i])
                end
            end

        end
        result = Dict(
            "globally" => Set([get_order_var(th, non_jet_ring)[1] for th in theta_g]),
            "locally_not_globally" => Set([get_order_var(th, non_jet_ring)[1] for th in setdiff(theta_l_new, theta_g)]),
            "nonidentifiable" => Set([get_order_var(th, non_jet_ring)[1] for th in setdiff(params_to_assess, theta_l)])
        )
        println("\n=== Summary ===")
        println("Globally identifiable parameters:                 [", join(result["globally"], ", "), "]")
        println("Locally but not globally identifiable parameters: [", join(result["locally_not_globally"], ", "), "]")
        println("Not identifiable parameters:                      [", join(result["nonidentifiable"], ", "), "]")
        println("===============")
        return result
    end
end
end

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2021 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

# includes -
include("Include.jl")

"""
    solve_static_problem()
"""
function solve_static_problem()::VLStaticSolution

    # create a default data dictionary -
    data_dictionary = generate_default_data_dictionary()

    # optional: do we need to change the problem setup to maximize the objective function?
    # is_minimum_flag = update_is_minimum_flag(data_dictionary)
    # data_dictionary["is_minimum_flag"] = is_minimum_flag

    # update the objective coefficient array -
    obj_coeff_array = update_objective_coefficient_array(data_dictionary)
    data_dictionary["objective_coefficient_array"] = obj_coeff_array

    # update the initial state arrary -
    species_concentration_array = update_species_concentration_array(data_dictionary)
    data_dictionary["species_concentration_array"] = species_concentration_array

    # update the species bounds array -
    species_bounds_array = update_species_bounds_array(data_dictionary)
    data_dictionary["species_bounds_array"] = species_bounds_array

    # update the flux bounds array -
    flux_bounds_array = update_flux_bounds_array(data_dictionary)
    data_dictionary["flux_bounds_array"] = flux_bounds_array

    # create a problem object (use the factory method in utilities ...)
    problem_object = build_static_fba_problem_object(data_dictionary)

    # solve the problem -
    return solve_simulation_problem(problem_object)
end

# execute -
@info("Start: solving the static flux balance analysis problem ...")
static_soln_object = solve_static_problem();
@info("Stop: calculation completed. The solution is stored in the static_soln_object workspace variable (which is of type VLStaticFBAProblem defined in Types.jl)")

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
struct AuxRotation{T}
    buff::Vector{T}  # each Thread must have an AuxRotation struct
end

export AuxRotation
export rotate!

"""
    rotate!(::AuxRotation, A::AbstractVector, tailpos::Integer)

Rotates `A` on `tailpos` using at most `n / 2` extra memory, 3n/2 copy-write ops.
"""
function rotate!(aux::AuxRotation, A::AbstractVector, tailpos::Integer)
    m = tailpos - 1
    n = length(A) - m
    if m <= n
        resize!(aux.buff, m)
        _copy!(aux.buff, 1, A, 1, m)
        _left_shift!(A, tailpos)
        _copy!(A, n+1, aux.buff, 1, m)
    else
        resize!(aux.buff, n)
        _copy!(aux.buff, 1, A, tailpos, n)
        _right_shift!(A, m)
        _copy!(A, 1, aux.buff, 1, n)
    end

    A
end


text:
__precompile__(true)
module Mongo

@static if VERSION < v"0.7-"
    const Nothing = Void
    const Cvoid   = Void
end

const deps_script = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(deps_script)
    error("Mongo.jl is not installed properly, run Pkg.build(\"Mongo\") and restart Julia.")
end
include(deps_script)
check_deps()

using LibBSON

ccall(
    (:mongoc_init, libmongoc),
    Cvoid, ()
    )

atexit() do
    ccall((:mongoc_cleanup, libmongoc), Cvoid, ())
end

const NakedDict = Union{Pair,Tuple}

include("MongoClient.jl")
include("MongoCollection.jl")
include("MongoCursor.jl")
include("query.jl")

end


text:
module MMS_Example

using EventSimulation
using Random
using Test
using Future: randjump

# Objectives of the example:
# * show how SimQueue and SimResource objects can be used
# * show that SimQueue gives a fine grained control than SimResource
#   but is more expensive (actually SimResource use is a hack)
# * show how MersenneTwister with randjump can be used to control RNG streams

# exact formula for M/M/s queue
function run_mms_exact(ar, sr, s)
    ρ = ar/sr
    P₀= 1 / (ρ^s/(factorial(s)*(1-ρ/s)) +
             sum(ρ^n/factorial(n) for n in 0:(s-1)))
    (ρ+(ρ^(s+1)/(factorial(s-1)*(s-ρ)^2))*P₀)/ar
end

# hardcoded M/M/1 queue
function run_mm1_fast(until, ar, sr, seed)
    start_time = time_ns()
    ma = MersenneTwister(seed)
    ms = randjump(ma, big(10)^20)
    queue = Vector{Float64}()
    nextArrival = randexp(ma) / ar
    nextDeparture = Inf
    totalWait = 0.0
    totalCount = 0
    now = 0.0
    msg = String[]
    while now < until
        if nextArrival < nextDeparture
            now = nextArrival
            if now <= until
                push!(msg, "A $now")
            end
            if isempty(queue)
                nextDeparture = nextArrival + randexp(ms) / sr
            end
            push!(queue, nextArrival)
            nextArrival += randexp(ma) / ar
        else
            now = nextDeparture
            if now <= until
                push!(msg, "D $now")
                totalWait += nextDeparture - popfirst!(queue)
                totalCount += 1
            end
            if isempty(queue)
                nextDeparture = Inf
            else
                nextDeparture += randexp(ms) / sr
            end
        end
    end
    println("MM1 fast time: ", (time_ns()-start_time)/10^9)
    return totalWait/totalCount, totalCount, msg
end

# Implementation of M/M/s queue using SimQueue
mutable struct StateQ <: AbstractState
    ar::Float64
    sr::Float64
    q::SimQueue{Float64}
    tcust::Int
    tin::Float64
    ma::MersenneTwister
    ms::MersenneTwister
    msg::Vector{String}
    function StateQ(ar, sr)
        ma = MersenneTwister(1)
        ms = randjump(ma, big(10)^20)
        new(Float64(ar), Float64(sr), SimQueue{Float64}(),
            0, 0.0, ma, ms, String[])
    end
end

function arrivalQ!(s::Scheduler)
    push!(s.state.msg, "A $(s.now)")
    provide!(s, s.state.q, s.now)
end

function startserviceQ!(s::Scheduler, a::Float64)
    register!(s, x -> endserviceQ!(x, a) , randexp(s.state.ms) / s.state.sr)
end

function endserviceQ!(s::Scheduler, a::Float64)
    s.state.tcust += 1
    s.state.tin += s.now - a
    push!(s.state.msg, "D $(s.now)")
    request!(s, s.state.q, startserviceQ!)
end

function run_mms_queue(until, ar, sr, count)
    start_time = time_ns()
    s = Scheduler(StateQ(ar, sr))
    repeat_register!(s, arrivalQ!, x -> randexp(x.state.ma) / x.state.ar)
    for i in 1:count
        request!(s, s.state.q, startserviceQ!)
    end
    go!(s, until)
    println("MM$count queue time: ", (time_ns() - start_time)/10^9)
    return s.state.tin/s.state.tcust, s.state.tcust, s.state.msg
end

# Implementation of M/M/s queue using SimResource
# Notice that if s>1 we have no control of arrival vs departure order
mutable struct StateR <: AbstractState
    ar::Float64
    sr::Float64
    r::SimResource{Int}
    ars::Vector{Float64}
    tcust::Int
    tin::Float64
    ma::MersenneTwister
    ms::MersenneTwister
    msg::Vector{String}
    function StateR(ar, sr)
        ma = MersenneTwister(1)
        ms = randjump(ma, big(10)^20)
        new(Float64(ar), Float64(sr), SimResource{Int}(),
            Vector{Float64}(), 0, 0.0, ma, ms, String[])
    end
end

function arrivalR!(s::Scheduler)
    register!(s, arrivalR!, randexp(s.state.ma) / s.state.ar)
    push!(s.state.ars, s.now)
    push!(s.state.msg, "A $(s.now)")
    provide!(s, s.state.r, 1)
end

function startserviceR!(s::Scheduler)
    register!(s, endserviceR!, randexp(s.state.ms) / s.state.sr)
end

function endserviceR!(s::Scheduler)
    s.state.tcust += 1
    # use s.state.tin only for count==1, otherwise it is inexact
    # as customer ordering migh have changed during the service
    s.state.tin += s.now - popfirst!(s.state.ars)
    push!(s.state.msg, "D $(s.now)")
    request!(s, s.state.r, 1, startserviceR!)
end

function run_mms_resource(until, ar, sr, count)
    start_time = time_ns()
    s = Scheduler(StateR(ar, sr))
    register!(s, arrivalR!, randexp(s.state.ma) / s.state.ar)
    for i in 1:count
        request!(s, s.state.r, 1, startserviceR!)
    end
    go!(s, until)
    println("MM$count resource time: ", (time_ns()-start_time)/10^9)
    # have to handle shortcomming of resource
    return (count == 1 ? s.state.tin/s.state.tcust : NaN,
            s.state.tcust, s.state.msg)
end

@testset "Test of M/M/1 queue" begin
    t_f = run_mm1_fast(100_000, 0.8, 1.4, 1)
    t_q = run_mms_queue(100_000, 0.8, 1.4, 1)
    t_r = run_mms_resource(100_000, 0.8, 1.4, 1)
    t_x = run_mms_exact(0.8, 1.4, 1)

    @test t_f == t_q == t_r
    @test isapprox(t_x, t_r[1], atol=0.1)
end

@testset "Test of M/M/2 queue" begin
    t2_q = run_mms_queue(100_000, 0.8, 0.7, 2)
    t2_r = run_mms_resource(100_000, 0.8, 0.7, 2)
    t2_x = run_mms_exact(0.8, 0.7, 2)

    @test t2_q[2:3] == t2_r[2:3]
    @test isapprox(t2_x, t2_q[1], atol=0.1)
end

@testset "Test of M/M/5 queue" begin
    t5_q = run_mms_queue(100_000, 0.8, 0.2, 5)
    t5_r = run_mms_resource(100_000, 0.8, 0.2, 5)
    t5_x = run_mms_exact(0.8, 0.2, 5)

    @test t5_q[2:3] == t5_r[2:3]
    @test isapprox(t5_x, t5_q[1], atol=0.2)
end

end # module


text:
#=
Direct estimation of the probability of for a diffusion with X(0) = x₀ satisfying
X(T) ∈ B for the Muller potential.
=#

using Statistics
using HypothesisTests
using Printf

include("muller_setup.jl");
opts = MDOptions(n_iters=nΔt)

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


text:
using SARA
using GaussianDistributions

# set up temperature profile
dT = 25 # uncertainty in peak temperature in C
dx = 50/1000 # uncertainty in position in mm
P = TemperatureProfile(dT, dx) # TODO: this needs to be updated before the run!

# define the range of temperatures that are taken from each stripe
t_low, t_high = 200, 10 # offset from T_peak
nout = 16
relevant_T = get_relevant_T(t_low, t_high, nout)

# define kernel used to calculate gradients per stripe (χAI)
l = .2 # length scale of inner GP
char_kernel = iSARA(l) # characterization kernel with length scale l

# TODO: interpolate xrd patterns

# convert stripe xrd data to global gradient data
# stripe_to_global(x, y, σ, k, P, T_max, log10_τ, relevant_T, input_noise)
# x is position
# Y is matrix of integrated xrd spectrograms as a function of position
# σ is measurement noise (can be learned using the estimate_noise function)
σ = 0.05
# T_max, log10_τ are stripe conditions
# relevant_T is function defined above
input_noise = Val(true) # switches on input noise propagation
condition = T_max, log10_τ
xrd_to_global(x, Y, σ, char_kernel, P, condition, relevant_T, input_noise)

# σAI
# stripe uncertainty acquisition function
policy = stripe_uncertainty_sampling(relevant_T)

# set up outer GP
T_length = 30.
τ_length = .3
c_length = .1 # disregard for systems without composition dimension
synth_kernel = oSARA(T_length, τ_length, c_length)
G = Gaussian(synth_kernel)

# standardize global gradient data before passing it to GP
y_normalization = std(y)
y /= y_normalization
σ² /= y_normalization^2 # noise variance

# G is GP
# x is vector / array of condition tuples
# y are corresponding gradients (vector of floats)
# variance is uncertainty in gradients (2nd return value of xrd_to_global)
# stripe_conditions is vector of potential conditions (tuples)
# policy is active learning aquisition function
next, G = next_stripe(G, x, y, variance, stripe_conditions, policy)

# evaluates GP
# function evaluate(G::Gaussian, n::Int = 256,
#                   T = range(200, 1200, length = n + 1),
#                   τ = range(2.5, 4, length = n),
#                   c = range(0, 1, length = 32))
#     return mean(G).(tuple.(T, τ', reshape(c, 1, 1, :)))
# end

# using HDF5
#
# # test workflow with LaMnOx data
# datapath = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/LaMnOx/"
# file = "LaMnOx_outer_xrd_data.h5"
# f = h5open(datapath * file)
# positions = read(f, "positions")
# temperatures = read(f, "temperatures")
# dwelltimes = read(f, "dwelltimes")
# compositions = read(f, "compositions")
# spectrograms = read(f, "spectrograms")
# close(f)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export container

"""
    container(;kwargs...)
    container(children::Any;kwargs...)
    container(children_maker::Function;kwargs...)

A Container component.
Center content horizontally with predefined max-width. For more information, see: https://mantine.dev/core/container/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Content that should be centered vertically and horizontally
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `fluid` (Bool; optional): If fluid is set to true, size prop is ignored and Container always take 100% of width
- `padding` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Horizontal padding defined in theme.spacing, or number value for padding in px
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Predefined container max-width or number for max-width in px
- `style` (Dict; optional): Inline style override
"""
function container(; kwargs...)
        available_props = Symbol[:children, :id, :className, :fluid, :padding, :size, :style]
        wild_props = Symbol[]
        return Component("container", "Container", "dash_mantine_components", available_props, wild_props; kwargs...)
end

container(children::Any; kwargs...) = container(;kwargs..., children = children)
container(children_maker::Function; kwargs...) = container(children_maker(); kwargs...)


text:
#=
Bit counting operations
=======================

Copyright (c) 2017 Ben J. Ward & Luis Yanes

This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
=#

"""
    count_nonzero_bitpairs(x::Unsigned)

Count the number of bit pairs in an unsigned integer `x` which have at least one bit set.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 19.
"""
@inline function count_nonzero_bitpairs(x::Unsigned)
    return count_ones((x | (x >> 1)) & repeatpattern(typeof(x), 0x55))
end

"""
    count_00_bitpairs(x::Unsigned)

Count the number of bit pairs in an unsigned integer `x` which have no bits set.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 13.
"""
@inline function count_00_bitpairs(x::Unsigned)
    return count_11_bitpairs(~x)
end

function count_zero_bitpairs(x::Unsigned)
    Base.depwarn("count_zero_bitpairs is deprecated, use count_00_bitpairs instead", :count_zero_bitpairs)
    return count_00_bitpairs(x)
end

"""
    count_11_bitpairs(x::Unsigned)

Count the number of bit pairs in an unsigned integer `x` which have both bits set.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 8.
"""
@inline function count_11_bitpairs(x::Unsigned)
    return count_ones((x & (x >> 1)) & repeatpattern(typeof(x), 0x55))
end

function count_one_bitpairs(x::Unsigned)
    Base.depwarn("count_one_bitpairs is deprecated, use count_11_bitpairs instead", :count_one_bitpairs)
    return count_11_bitpairs(x)
end

"""
    count_01_bitpairs(x::Unsigned)

Count the number of bit pairs in an unsigned integer `x` which have just their LSB set.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 11.
"""
@inline function count_01_bitpairs(x::Unsigned)
    return count_ones((((~x) >> 1) & x) & repeatpattern(typeof(x), 0x55))
end

"""
    count_10_bitpairs(x::Unsigned)

Count the number of bit pairs in an unsigned integer `x` which have just their MSB set.

E.g. An input of:

0x8F11F111F11A11F1

Would give the answer: 3.
"""
@inline function count_10_bitpairs(x::Unsigned)
    return count_ones(((x >> 1) & (~x)) & repeatpattern(typeof(x), 0x55))
end

"""
    count_nonzero_nibbles(x::Unsigned)

Count the number of nibbles (aligned 4 bit segments) in an unsigned integer `x`
which have at least one bit set.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 15.
"""
@inline function count_nonzero_nibbles(x::Unsigned)
    return count_ones((x | (x >> 1) | (x >> 2) | (x >> 3)) & repeatpattern(typeof(x), 0x11))
end

"""
    count_0000_nibbles(x::Unsigned)

Counts the number of nibbles (aligned 4 bit segments) in an unsigned integer `x`
that have all their bits unset i.e.
nibbles of 0000.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 1.
"""
@inline function count_0000_nibbles(x::Unsigned)
    return count_1111_nibbles(~x)
end

@inline function count_zero_nibbles(x::Unsigned)
    Base.depwarn("count_zero_nibbles is deprecated, use count_0000_nibbles instead", :count_zero_nibbles)
    return count_0000_nibbles(x)
end

"""
    count_1111_nibbles(x::Unsigned)

Counts the number of nibbles (aligned 4 bit segments) in an unsigned integer `x`
that have all their bits set i.e. counts all nibbles of 1111 in an integer.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 4.
"""
@inline function count_1111_nibbles(x::Unsigned)
    return count_ones(x & (x >> 1) & (x >> 2) & (x >> 3) & repeatpattern(typeof(x), 0x11))
end

@inline function count_one_nibbles(x::Unsigned)
    Base.depwarn("count_one_nibbles is deprecated, use count_1111_nibbles instead", :count_one_nibbles)
    return count_1111_nibbles(x)
end


text:
global tophap_ccf_mask_scale_factor=1.6

global max_spectra_to_use = 200
global fits_target_str
if fits_target_str == "Solar"
   global linelist_for_ccf_filename = "G2.espresso.mas"
   hostname = gethostname()
   if occursin("aci.ics.psu.edu",hostname)
      global ancilary_solar_data_path = "/gpfs/group/ebf11/default/ebf11/neid_solar"
   elseif occursin("nuc8",hostname)  # Eric's home machine :)
      global ancilary_solar_data_path = "/home/eford/Data/SolarSpectra/NEID_solar/"
   end
   global ccf_mid_velocity = 0
   global bjd_first_good = 2458745.1296134139
   global bjd_last_good = 2458745.283
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @filter(bjd_first_good <= _.bjd < bjd_last_good) |>
      @take(max_spectra_to_use) |>
      DataFrame
elseif fits_target_str == "101501"
   global linelist_for_ccf_filename = "G8.espresso.mas"
   global ccf_mid_velocity = -5e3
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end


text:
using Test
include("components.jl")
include("registry.jl")

max_stars_repo_name: waralex/DashBase

max_issues_repo_name: waralex/DashBase

max_forks_repo_name: waralex/DashBase


text:
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Stats: I, mean, var, std, cor
@testset "Stats" begin
    @testset "Identity Matrix" begin
        x = I(3)

        @test x isa Array{GreyNumber{Float64}, 2}
        @test x[1, 1] isa GreyNumber

        @test x[1, 1] == GreyNumber{Float64}(1.0, 1.0)
        @test x[2, 2] == GreyNumber{Float64}(1.0, 1.0)
        @test x[3, 3] == GreyNumber{Float64}(1.0, 1.0)
        @test x[1, 2] == GreyNumber{Float64}(0.0, 0.0)
        @test x[1, 3] == GreyNumber{Float64}(0.0, 0.0)
        @test x[2, 1] == GreyNumber{Float64}(0.0, 0.0)
        @test x[2, 3] == GreyNumber{Float64}(0.0, 0.0)
        @test x[3, 1] == GreyNumber{Float64}(0.0, 0.0)
        @test x[3, 2] == GreyNumber{Float64}(0.0, 0.0)
    end

    @testset "Testing mean" begin
        x = [
            GreyNumber{Float64}(1.0, 1.0),
            GreyNumber{Float64}(2.0, 2.0),
            GreyNumber{Float64}(3.0, 3.0),
            GreyNumber{Float64}(4.0, 4.0),
            GreyNumber{Float64}(5.0, 5.0),
        ]
        mx = mean(x)
        @test mx == GreyNumber{Float64}(3.0, 3.0)
    end

    @testset "Testing var" begin
        x = [
            GreyNumber{Float64}(1.0, 1.0),
            GreyNumber{Float64}(2.0, 2.0),
            GreyNumber{Float64}(3.0, 3.0),
            GreyNumber{Float64}(4.0, 4.0),
            GreyNumber{Float64}(5.0, 5.0),
        ]
        varx = var(x)
        @test varx == GreyNumber{Float64}(2.5, 2.5)
    end

    @testset "Testing std" begin
        tol = 0.001
        x = [
            GreyNumber{Float64}(1.0, 1.0),
            GreyNumber{Float64}(2.0, 2.0),
            GreyNumber{Float64}(3.0, 3.0),
            GreyNumber{Float64}(4.0, 4.0),
            GreyNumber{Float64}(5.0, 5.0),
        ]
        stdx = std(x)
        @test isapprox(GreyNumber{Float64}(1.58113, 1.58113), stdx, atol = tol)
    end

    @testset "Testing cor" begin
        tol = 0.001
        x = [
            GreyNumber{Float64}(1.0, 1.0),
            GreyNumber{Float64}(2.0, 2.0),
            GreyNumber{Float64}(3.0, 3.0),
            GreyNumber{Float64}(4.0, 4.0),
            GreyNumber{Float64}(5.0, 5.0),
        ]
        y = reverse(x)

        cor1 = cor(x, x)
        @test isapprox(GreyNumber{Float64}(1.0, 1.0), cor1, atol = tol)

        cor2 = cor(x, y)
        @test isapprox(GreyNumber{Float64}(-1.0, -1.0), cor2, atol = tol)
    end
end


text:
module Schistoxpkg
include("parameters.jl")
include("functions.jl")

export Environment
export Human
export Parameters
export create_population
export make_age_contact_rate_array
export make_death_rate_array
export birth_of_human
export cercariae_uptake!
export worm_maturity!
export calculate_worm_pairs
export calculate_total_worms
export egg_production!
export egg_production_increasing!
export miracidia_production!
export death_of_human
export birth_of_human
export count_eggs
export human_cercariae_maturity
export update_env
export find_death_rate
export cercariae_death!
export miracidia_death!
export mda
export administer_drug
export administer_vaccine
export update_contact_rate
export update_death_rate
export update_vaccine
export vaccinate
export vac_decay!
export update_mda
export load_population_from_file
export kato_katz
export get_prevalences!
export save_population_to_file
export run_simulation_from_loaded_population
export find_mean_worms_by_age
export run_simulation
export generate_age_distribution
export specified_age_distribution
export create_population_specified_ages
export update_env_to_equilibrium
export create_mda
export vaccine_information
export mda_information
export update_env_no_births_deaths
export update_env_keep_population_same
export collect_prevs
export run_repeated_sims_no_population_change
# export run_repeated_sims_random_births_deaths
export run_repeated_sims_no_population_change_increasing
# export run_repeated_sims_random_births_deaths_increasing
export create_contact_settings
export get_death_age
export generate_ages_and_deaths
export out
export update_env_keep_population_same_save_predisposition
export death_of_human_save_predis
export birth_of_human_specified_predis
export update_env_with_mda_no_births_deaths
export run_repeated_sims_no_births_deaths_increasing
export run_repeated_sims_no_births_deaths
export update_env_constant_population
export update_env_to_equilibrium_increasing
export update_env_constant_population_increasing
export update_env_no_births_deaths_increasing
export cercariae_uptake_with_human_larvae!
export human_larvae_maturity
export update_env_to_equilibrium_human_larvae
export update_env_constant_population_human_larvae
export update_env_no_births_deaths_human_larvae
export run_repeated_sims_no_population_change_human_larvae
export run_repeated_sims_no_births_deaths_human_larvae
export enact_maturity_function_false
export enact_maturity_function_true
# using Distributions
# using Random
# using PoissonRandom

end # module


text:
# sub(aln, :, i) is ~ 0.015 seconds faster than aln[:,i] and better with the memmory (PF00085)
# aln.msa' takes ~ 0.015 seconds in PF00085
# Using aln.msa instead of aln in estimate_on_column_pairs is ~ 0.53 seconds faster for PF00085

"""
`estimateincolumns(aln, [count,] use, [α, β,] measure, [pseudocount,] [weight,] [usediagonal, diagonalvalue])`

This function `estimate` a `AbstractMeasure` in columns or pair of columns of a MSA.

- `aln` : This argument is mandatory and it should be a `Matrix{Residue}`. Use the function `getresidues` (from the MSA module) over a MSA object to get the needed matrix.
- `count` : This argument is optional. It should be defined when `use` is a `ResidueProbability` object. It indicates the element type of the counting table.
- `use` : This argument is mandatory and indicates the sub-type of `ResidueContingencyTables` used by `estimate` inside the function.
If the table has one dimension (`N`=`1`), the occurrences/probabilities are counted for each sequence/column.
If the table has two dimension (`N`=`2`), pairs of sequences/columns are used.
The dimension `N` and the `UseGap` parameter of `Residueount{T, N, UseGap}` or `ResidueProbability{T, N, UseGap}` determines the output and behaviour of this functions.
If `UseGap` is true, gaps are used in the estimations.
- `α` : This argument is optional, and indicates the weight of real frequencies to apply BLOSUM62 based pseudo frequencies.
- `β` : This argument is optional, and indicates the weight of BLOSUM62 based pseudo frequencies.
- `measure` : This argument is mandatory and indicates the measure to be used by `estimate` inside the function.
- `pseudocount` : This argument is optional. It should be an `AdditiveSmoothing` instance (default to zero).
- `weight` : This argument is optional. It should be an instance of `ClusteringResult` or `AbstractVector` (vector of weights).
Each sequence has weight 1 (`NoClustering()`) by default.
- `usediagonal` : This functions return a `Vector` in the one dimensional case, or a `PairwiseListMatrix` in the bidimensional case.
This argument only have sense in the bidimensional case and indicates if the list on the `PairwiseListMatrix` should include the diagonal (default to `true`).
- `diagonalvalue` : This argument is optional (default to zero). Indicates the value of output diagonal elements.
"""
function estimateincolumns{T, TP, UseGap}(aln::AbstractMatrix{Residue}, use::Type{ResidueCount{T, 1, UseGap}}, measure::AbstractMeasure{TP},
                                          pseudocount::Pseudocount{T}=zero(AdditiveSmoothing{T}), weight::SequenceWeights=NoClustering())
  N = ResidueCount{T, 1, UseGap}()
  ncol = ncolumns(aln)
  scores = Array(TP, ncol)
  @inbounds for i in 1:ncol
    fill!(N, pseudocount) # instead of fill! with 0
    count!(N, weight, sub(aln,:,i))
    scores[i] = estimate(measure, N)
  end
  scores
end

function estimateincolumns{T <: Real, TP, UseGap}(aln::AbstractMatrix{Residue}, count::Type{T}, use::Type{ResidueProbability{TP, 1, UseGap}}, measure::AbstractMeasure{TP},
                                                  pseudocount::Pseudocount{T}=zero(AdditiveSmoothing{T}), weight::SequenceWeights=NoClustering())
  N = ResidueCount{T, 1, UseGap}()
  P = ResidueProbability{TP, 1, UseGap}()
  ncol = ncolumns(aln)
  scores = Array(TP, ncol)
  @inbounds for i in 1:ncol
    fill!(N, pseudocount) # instead of fill! with 0
    count!(N, weight, sub(aln,:,i))
    fill!(P, N) # count! calls update!
    scores[i] = estimate(measure, P)
  end
  scores
end

function estimateincolumns{T, TP, UseGap}(aln::AbstractMatrix{Residue}, use::Type{ResidueCount{T, 2, UseGap}}, measure::SymmetricMeasure{TP},
                                          pseudocount::Pseudocount{T}=zero(AdditiveSmoothing{T}), weight::SequenceWeights=NoClustering(), usediagonal::Bool=true, diagonalvalue::TP=zero(TP))
  Nab = ResidueCount{T, 2, UseGap}()
  ncol = ncolumns(aln)
  scores = columnpairsmatrix(aln, TP, usediagonal, diagonalvalue) # Array(TP, ncol, ncol)
  @inbounds for i in 1:ncol
    a = sub(aln,:,i)
    for j in i:ncol
      if !usediagonal && i == j
        continue
      end
      b = sub(aln,:,j)
      fill!(Nab, pseudocount) # instead of fill! with 0
      count!(Nab, weight, a, b)
      scores[i,j] = estimate(measure, Nab)
    end
  end
  scores
end

function estimateincolumns{T <: Real, TP, UseGap}(aln::AbstractMatrix{Residue}, count::Type{T}, use::Type{ResidueProbability{TP, 2, UseGap}},
                                                  measure::SymmetricMeasure{TP}, pseudocount::Pseudocount{T}=zero(AdditiveSmoothing{T}),
                                                  weight::SequenceWeights=NoClustering(), usediagonal::Bool=true, diagonalvalue::TP=zero(TP))
  Nab = ResidueCount{T, 2, UseGap}()
  Pab = ResidueProbability{TP, 2, UseGap}()
  ncol = ncolumns(aln)
  scores = columnpairsmatrix(aln, TP, usediagonal, diagonalvalue) # zeros(TP, ncol, ncol)
  @inbounds for i in 1:ncol
    a = sub(aln,:,i)
    for j in i:ncol
      if !usediagonal && i == j
        continue
      end
      b = sub(aln,:,j)
      fill!(Nab, pseudocount) # instead of fill! with 0
      count!(Nab, weight, a, b)
      fill!(Pab, Nab)
      scores[i,j] = estimate(measure, Pab)
    end
  end
  scores
end

function estimateincolumns{T <: Real, TP}(aln::AbstractMatrix{Residue}, count::Type{T}, use::Type{ResidueProbability{TP, 2, false}}, α, β,
                                          measure::SymmetricMeasure{TP}, pseudocount::Pseudocount{T}=zero(AdditiveSmoothing{T}),
                                          weight::SequenceWeights=NoClustering(), usediagonal::Bool=true, diagonalvalue::TP=zero(TP))
  Nab = ResidueCount{T, 2, false}()
  Pab = ResidueProbability{TP, 2, false}()
  Gab = ResidueProbability{TP, 2, false}()
  ncol = ncolumns(aln)
  scores = columnpairsmatrix(aln, TP, usediagonal, diagonalvalue) # zeros(TP, ncol, ncol)
  @inbounds for i in 1:ncol
    a = sub(aln,:,i)
    for j in i:ncol
      if !usediagonal && i == j
        continue
      end
      b = sub(aln,:,j)
      fill!(Nab, pseudocount) # instead of fill! with 0
      count!(Nab, weight, a, b)
      fill!(Pab, Nab)
      blosum_pseudofrequencies!(Gab, Pab)
      apply_pseudofrequencies!(Pab, Gab, α, β)
      scores[i,j] = estimate(measure, Pab)
    end
  end
  scores
end

"""
This function `estimate` a `measure` over sequences or sequence pairs.
It has the same arguments than `estimateincolumns`, look the documentation of the last.
"""
estimateinsequences(aln::AbstractMatrix{Residue}, args...) = estimateincolumns(transpose(aln), args...)

# cMI
# ===

"""
`cumulative` allows to calculate cumulative scores (i.e. cMI) as defined in Buslje et. al. 2010

*"We calculated a cumulative mutual information score (cMI) for each residue
as the sum of MI values above a certain threshold for every amino acid pair where the particular residue appears.
This value defines to what degree a given amino acid takes part in a mutual information network."*
Buslje, Cristina Marino, Elin Teppa, Tomas Di Doménico, José María Delfino, and Morten Nielsen.
*Networks of high mutual information define the structural proximity of catalytic sites: implications for catalytic residue identification.*
PLoS Comput Biol 6, no. 11 (2010): e1000978.
"""
function cumulative{T, UseDiag}(mat::PairwiseListMatrix{T, UseDiag}, threshold::T)
    nrow, ncol = size(mat)
    cum = zeros(T, ncol)
    @inbounds for col in 1:ncol
        for row in 1:nrow
            if col != row
                elem = mat[row, col]
                if !isnan(elem) && elem >= threshold
                    cum[col] += elem
                end
            end
        end
    end
    cum
end