text:
using NonlinearEigenproblems.RKHelper
using LinearAlgebra
using Random

export nleigs, NleigsSolutionDetails

"""
    nleigs(nep::NEP, Σ::AbstractVector{Complex{T}})

Find a few eigenvalues and eigenvectors of a nonlinear eigenvalue problem, using the `nleigs` algorithm.

# Arguments
- `nep`: An instance of a nonlinear eigenvalue problem.
- `Σ`: A vector containing the points of a polygonal target set in the complex plane.
- `Ξ`: A vector containing a discretization of the singularity set.
- `logger`: Level of display (0, 1, 2).
- `maxdgr`: Max degree of approximation.
- `minit`: Min number of iterations after linearization is converged.
- `maxit`: Max number of total iterations.
- `tol`: Tolerance for residual.
- `tollin`: Tolerance for convergence of linearization.
- `v`: Starting vector.
- `errmeasure`: Function for error measure (residual norm). Called with arguments (λ,v).
- `isfunm` : Whether to use matrix functions.
- `static`: Whether to use static version of NLEIGS.
- `leja`: Use of Leja-Bagby points (0 = no, 1 = only in expansion phase, 2 = always).
- `nodes`: Prefixed interpolation nodes (only when leja is 0 or 1).
- `reusefact`: Reuse of matrix factorizations (0 = no, 1 = only after converged linearization, 2 = always).
- `blksize`: Block size for pre-allocation.
- `return_details`: Whether to return solution details (see NleigsSolutionDetails).
- `check_error_every`: Check for convergence / termination every this number of iterations.

See [`augnewton`](@ref) for other parameters.

# Return values
- `λ`: Vector of eigenvalues of the nonlinear eigenvalue problem NLEP inside the target set Σ.
- `X`: Corresponding matrix of eigenvectors.
- `res`: Corresponding residuals.
- `details`: Solution details, if requested (see NleigsSolutionDetails).

# Example
```julia-repl
julia> nep=nep_gallery("dep0");
julia> unit_square = float([1+1im, 1-1im, -1-1im,-1+1im])
julia> (λ,v)=nleigs(nep,unit_square);
julia> norm(compute_Mlincomb(nep,λ[1],v[:,1]))
5.028313023882492e-14
julia> norm(compute_Mlincomb(nep,λ[2],v[:,2]))
1.1937025845487509e-13
```
# References
- S. Guettel, R. Van Beeumen, K. Meerbergen, and W. Michiels. NLEIGS: A class
  of fully rational Krylov methods for nonlinear eigenvalue problems. SIAM J.
  Sci. Comput., 36(6), A2842-A2864, 2014.
- [NLEIGS Matlab toolbox](http://twr.cs.kuleuven.be/research/software/nleps/nleigs.php) (GPL License)
- [NLEIGS Matlab toolbox](https://bitbucket.org/roelvb/nleigs) (MIT License)
"""
nleigs(nep, Σ; params...) = nleigs(Float64, nep, Σ; params...)
function nleigs(
        ::Type{T},
        nep::NEP,
        Σ::AbstractVector{CT}=Vector{CT}([-1.0-1im,-1+1im,+1+1im,1-1im]);
        Ξ::Vector{T} = [T(Inf)],
        logger = 0,
        maxdgr::Int = 100,
        minit::Int = 20,
        maxit::Int = 200,
        linsolvercreator=DefaultLinSolverCreator(),
        tol::T = 1e-10,
        tollin::T = max(tol/10, 100*eps(T)),
        v::Vector{CT} = CT.(randn(T, size(nep, 1))),
        errmeasure::ErrmeasureType = ResidualErrmeasure(nep),
        isfunm::Bool = true,
        static::Bool = false,
        leja::Int = 1,
        nodes::Vector{CT} = Vector{CT}(),
        reusefact::Int = 1,
        blksize::Int = 20,
        return_details::Bool = false,
        check_error_every::Int = 5) where {T<:Real, CT<:Complex{T}}

    @parse_logger_param!(logger)

    # The following variables are used when creating the return values, so put them in scope
    D = Vector{Matrix{CT}}()
    conv = BitVector()
    lam = Vector{CT}()
    X = Matrix{CT}(undef, 0, 0)
    res = Vector{T}()

    P = get_rk_nep(T, nep)
    n = size(nep, 1)
    n == 1 && (maxdgr = maxit + 1)
    computeD = (n <= 400) # for small problems, explicitly use generalized divided differences
    b = blksize

    lin_solver_cache = LinSolverCache(CT, nep, linsolvercreator)

    # Initialization
    if static
        V = zeros(CT, n, 1)
    elseif !P.spmf
        V = zeros(CT, (b+1)*n, b+1)
    else
        if !P.is_low_rank || b < P.p
            V = zeros(CT, (b+1)*n, b+1)
        else
            V = zeros(CT, P.p*n+(b-P.p+1)*P.r, b+1)
        end
    end
    H = zeros(CT, b+1, b)
    K = zeros(CT, b+1, b)
    nrmD = Array{T}(undef, 1)
    if return_details
        Lam = zeros(CT, b, b)
        Res = zeros(T, b, b)
    end

    # Discretization of Σ --> Gamma & Leja-Bagby points
    if leja == 0 # use no leja nodes
        if isempty(nodes)
            error("Interpolation nodes must be provided via 'nodes' when no Leja-Bagby points ('leja' == 0) are used.")
        end
        gamma,_ = discretizepolygon(Σ)
        max_count = static ? maxit+maxdgr+2 : max(maxit,maxdgr)+2
        σ = repeat(reshape(nodes, :, 1), ceil(Int, max_count/length(nodes)), 1)
        _,ξ,β = lejabagby(σ[1:maxdgr+2], Ξ, gamma, maxdgr+2, true, P.p)
    elseif leja == 1 # use leja nodes in expansion phase
        if isempty(nodes)
            gamma,nodes = discretizepolygon(Σ, true)
        else
            gamma,_ = discretizepolygon(Σ)
        end
        nodes = repeat(reshape(nodes, :, 1), ceil(Int, (maxit+1)/length(nodes)), 1)
        σ,ξ,β = lejabagby(gamma, Ξ, gamma, maxdgr+2, false, P.p)
    else # use leja nodes in both phases
        gamma,_ = discretizepolygon(Σ)
        max_count = static ? maxit+maxdgr+2 : max(maxit,maxdgr)+2
        σ,ξ,β = lejabagby(gamma, Ξ, gamma, max_count, false, P.p)
    end
    ξ[maxdgr+2] = NaN # not used
    if (!P.spmf || !isfunm) && length(σ) != length(unique(σ))
        error("All interpolation nodes must be distinct when no matrix " *
            "functions are used for computing the generalized divided differences.")
    end

    # Rational Newton coefficients
    range = 1:maxdgr+2
    if !P.spmf
        D = ratnewtoncoeffs(λ -> compute_Mder(nep, λ[1]), σ[range], ξ[range], β[range])
        nrmD[1] = norm(D[1]) # Frobenius norm
        sgdd = Matrix{CT}(undef, 0, 0)
    else
        # Compute scalar generalized divided differences
        sgdd = scgendivdiffs(σ[range], ξ[range], β[range], maxdgr, isfunm, get_fv(nep))
        # Construct first generalized divided difference
        computeD && push!(D, constructD(0, P, sgdd))
        # Norm of first generalized divided difference
        nrmD[1] = maximum(abs.(sgdd[:,1]))
    end
    if !isfinite(nrmD[1]) # check for NaN
        error("The generalized divided differences must be finite.");
    end

    # Rational Krylov
    v = solve(lin_solver_cache, σ[1], v/norm(v), reusefact == 2)
    V[1:n,1] .= v ./ norm(v)
    expand = true
    kconv = trunc(Int, typemax(Int)/2)
    kn = n   # length of vectors in V
    l = 0    # number of vectors in V
    N = 0    # degree of approximations
    nbconv = 0 # number of converged lambdas inside Σ
    nblamin = 0 # number of lambdas inside Σ, converged or not
    kmax = static ? maxit + maxdgr : maxit
    k = 1
    while k <= kmax
        # resize matrices if we're starting a new block
        if l > 0 && (b == 1 || mod(l+1, b) == 1)
            nb = round(Int, 1 + l/b)
            Vrows = size(V, 1)
            if !P.spmf
                Vrows = kn+b*n
            else
                if expand
                    if !P.is_low_rank || l + b < P.p
                        Vrows = kn+b*n
                    elseif l < P.p-1
                        Vrows = P.p*n+(nb*b-P.p+1)*P.r
                    else # l >= P.p-1
                        Vrows = kn+b*P.r
                    end
                end
            end
            V = resize_matrix(V, Vrows, nb*b+1)
            H = resize_matrix(H, size(H, 1) + b, size(H, 2) + b)
            K = resize_matrix(K, size(K, 1) + b, size(K, 2) + b)
            if return_details
                Lam = resize_matrix(Lam, size(Lam, 1) + b, size(Lam, 2) + b)
                Res = resize_matrix(Res, size(Res, 1) + b, size(Res, 2) + b)
            end
        end

        if expand
            # set length of vectors
            if !P.is_low_rank || k < P.p
                kn += n
            else
                kn += P.r
            end

            # rational divided differences
            if P.spmf && computeD
                push!(D, constructD(k, P, sgdd))
            end
            N += 1

            # monitoring norms of divided difference matrices
            if !P.spmf
                push!(nrmD, norm(D[k+1])) # Frobenius norm
            else
                # The below can cause out of bounds in sgdd if there's
                # no convergence (also happens in MATLAB implementation)
                push!(nrmD, maximum(abs.(sgdd[:,k+1])))
            end
            if !isfinite(nrmD[k+1]) # check for NaN
                error("The generalized divided differences must be finite.");
            end
            if n > 1 && k >= 5 && k < kconv
                if sum(nrmD[k-3:k+1]) < 5*tollin
                    kconv = k - 1
                    if static
                        kmax = maxit + kconv
                    end
                    expand = false
                    if leja == 1
                        # TODO: can we pre-allocate σ?
                        if length(σ) < kmax+1
                            resize!(σ, kmax+1)
                        end
                        σ[k+1:kmax+1] = nodes[1:kmax-k+1]
                    end
                    if !P.spmf || computeD
                        D = D[1:k]
                    end
                    ξ = ξ[1:k]
                    β = β[1:k]
                    nrmD = nrmD[1:k]
                    if static
                        if !P.is_low_rank || k < P.p
                            kn -= n
                        else
                            kn -= P.r
                        end
                        V = resize_matrix(V, kn, b+1)
                    end
                    N -= 1
                    push_info!(logger,
                               "Linearization converged after $kconv iterations")
                    push_info!(logger,
                               "--> freeze linearization")
                elseif k == maxdgr+1
                    kconv = k
                    expand = false
                    if leja == 1
                        if length(σ) < kmax+1
                            resize!(σ, kmax+1)
                        end
                        σ[k+1:kmax+1] = nodes[1:kmax-k+1]
                    end
                    if static
                        V = resize_matrix(V, kn, b+1)
                    end
                    N -= 1
                    @warn "NLEIGS: Linearization not converged after $maxdgr iterations"
                    push_info!(logger,
                               "--> freeze linearization")
                end
            end
        end

        l = static ? k - N : k

        if !static || (static && !expand)
            # shift-and-invert
            t = [zeros(l-1); 1]    # continuation combination
            wc = V[1:kn, l]        # continuation vector
            w = backslash(wc, P, lin_solver_cache, reusefact, computeD, σ, k, D, β, N, ξ, expand, kconv, sgdd)

            # orthogonalization
            Vview = view(V, 1:kn, 1:l)
            H[l+1,l] = orthogonalize_and_normalize!(Vview, w, view(H, 1:l,l), DGKS)
            K[1:l,l] .= view(H, 1:l, l) .* σ[k+1] .+ t
            K[l+1,l] = H[l+1,l] * σ[k+1]
            V[1:kn,l+1] = w
#            @printf("new vector V: size = %s, sum = %s\n", size(V[1:kn,l+1]), sum(sum(V[1:kn,l+1])))
        end

        function check_convergence(all)
            lambda, S = eigen(K[1:l,1:l], H[1:l,1:l])

            # select eigenvalues
            if !all
                lamin = in_Σ(lambda, Σ, tol)
                ilam = [1:l;][lamin]
                lam = lambda[ilam]
            else
                ilam = [1:l;][isfinite.(lambda)]
                lam = lambda[ilam]
                lamin = in_Σ(lam, Σ, tol)
            end

            nblamin = sum(lamin)
            for i = ilam
                S[:,i] /= norm(H[1:l+1,1:l] * S[:,i])
            end
            X = V[1:n,1:l+1] * (H[1:l+1,1:l] * S[:,ilam])
            for i = 1:size(X,2)
                X[:,i] /= norm(X[:,i])
            end

            # compute residuals & check for convergence
            res = map(i -> estimate_error(errmeasure,lam[i], X[:,i]), 1:length(lam))
            conv = abs.(res) .< tol
            if all
                resall = fill(T(NaN), l, 1)
                resall[ilam] = res
                # sort complex numbers by magnitude, then angle
                si = sortperm(lambda, lt = (a,b) -> abs(a) < abs(b) || (abs(a) == abs(b) && angle(a) < angle(b)))
                Res[1:l,l] = resall[si]
                Lam[1:l,l] = lambda[si]
                conv .&= lamin
            end

            nbconv = isempty(conv) ? 0 : sum(conv)

            iteration = static ? k - N : k
            push_info!(logger,
                       "  iteration $iteration: $nbconv of $nblamin < $tol")
        end

        # Ritz pairs
        if !return_details && (
            (!expand && k >= N + minit && mod(k-(N+minit), check_error_every) == 0) ||
            (k >= kconv + minit && mod(k-(kconv+minit), check_error_every) == 0) || k == kmax)
            check_convergence(false)
        elseif return_details && (!static || (static && !expand))
            check_convergence(true)
        end

        # stopping
        if ((!expand && k >= N + minit) || k >= kconv + minit) && nblamin == nbconv
            break
        end

        # increment k
        k += 1
    end

    details = NleigsSolutionDetails{T,CT}()

    if return_details
        Lam = Lam[1:l,1:l]
        Res = Res[1:l,1:l]
        σ = σ[1:k]
        if expand
            ξ = ξ[1:k]
            β = β[1:k]
            nrmD = nrmD[1:k]
            @warn "NLEIGS: Linearization not converged after $maxdgr iterations"
        end
        details = NleigsSolutionDetails(Lam, Res, σ, ξ, β, nrmD, kconv)
    end

    return lam[conv], X[:,conv], res[conv], details
end

"Construct generalized divided difference for number `nb`."
function constructD(nb, P, sgdd::AbstractMatrix{CT}) where CT<:Complex{<:Real}
    n = size(P.nep, 1)
    BC = get_Av(P.nep)
    if !P.is_low_rank || nb <= P.p
        D = spzeros(CT, n, n)
        for ii = 1:(P.p+1+P.q)
            D += sgdd[ii,nb+1] * BC[ii]
        end
    else
        D = []
        for ii = 1:P.q
            d = sgdd[P.p+1+ii,nb+1] * P.L[ii]
            D = ii == 1 ? d : hcat(D, d)
        end
    end
    return D
end

"Backslash or left matrix divide for continuation vector `wc`."
function backslash(wc, P, lin_solver_cache, reusefact, computeD, σ, k, D, β, N, ξ, expand, kconv, sgdd)
    n = size(P.nep, 1)::Int
    shift = σ[k+1]

    # construction of B*wc
    Bw = zeros(eltype(wc), size(wc))
    # first block (if low rank)
    if P.is_low_rank
        i0b = (P.p-1)*n + 1
        i0e = P.p*n
        if !P.spmf || computeD
            Bw[1:n] = -D[P.p+1] * wc[i0b:i0e] / β[P.p+1]
        else
            Bw[1:n] = -sum(reshape(P.BBCC * wc[i0b:i0e], n, :) .* transpose(sgdd[:,P.p+1]), dims = 2) / β[P.p+1];
        end
    end
    # other blocks
    i0b = 1
    i0e = n
    for ii = 1:N
        # range of block i+1
        i1b = i0e + 1
        if !P.is_low_rank || ii < P.p
            i1e = i0e + n
        else
            i1e = i0e + P.r
        end
        # compute block i+1
        if !P.is_low_rank || ii != P.p
            Bw[i1b:i1e] = wc[i0b:i0e] + β[ii+1]/ξ[ii]*wc[i1b:i1e]
        else
            Bw[i1b:i1e] = P.UU' * wc[i0b:i0e] + β[ii+1]/ξ[ii]*wc[i1b:i1e]
        end
        # range of next block i
        i0b = i1b
        i0e = i1e
    end

    # construction of z0
    z = copy(Bw)
    i1b = n + 1
    if !P.is_low_rank || P.p > 1
        i1e = 2*n
    else
        i1e = n + P.r
    end
    nu = β[2]*(1 - shift/ξ[1])
    z[i1b:i1e] = 1/nu * z[i1b:i1e]
    for ii = 1:N
        # range of block i+2
        i2b = i1e + 1
        if !P.is_low_rank || ii < P.p-1
            i2e = i1e + n
        else
            i2e = i1e + P.r
        end
        # add extra term to z0
        if !P.spmf || computeD
            if !P.is_low_rank || ii != P.p
                z[1:n] -= D[ii+1] * z[i1b:i1e]
            end
        else
            if !P.is_low_rank || ii < P.p
                z[1:n] -= sum(reshape(P.BBCC * z[i1b:i1e], n, :) .* transpose(sgdd[:,ii+1]), dims = 2)
            elseif ii > P.p
                dd = sgdd[P.p+2:end,ii+1]
                @inbounds for i = 1:length(P.iLr)
                    for j = 1:length(P.LL[i].nzval)
                        ind = P.LL[i].nzind[j]
                        z[P.iLr[i]] -= P.LL[i].nzval[j] * z[i1b+ind-1] * dd[P.iL[ind]]
                    end
                end
            end
        end
        # update block i+2
        if ii < N
            mu = shift - σ[ii+1]
            nu = β[ii+2] * (1 - shift/ξ[ii+1])
            if !P.is_low_rank || ii != P.p-1
                z[i2b:i2e] = 1/nu * z[i2b:i2e] + mu/nu * z[i1b:i1e]
            else # i == p-1
                z[i2b:i2e] = 1/nu * z[i2b:i2e] + mu/nu * P.UU'*z[i1b:i1e]
            end
        end
        # range of next block i+1
        i1b = i2b
        i1e = i2e
    end

    # solving Alam x0 = z0
    w = zeros(eltype(wc), size(wc))
    add_to_cache = ((!expand || k > kconv) && reusefact == 1) || reusefact == 2
    w[1:n] = solve(lin_solver_cache, shift, z[1:n]/β[1], add_to_cache)

    # substitutions x[i+1] = mu/nu*x[i] + 1/nu*Bw[i+1]
    i0b = 1
    i0e = n
    for ii = 1:N
        # range of block i+1
        i1b = i0e + 1
        if !P.is_low_rank || ii < P.p
            i1e = i0e + n
        else
            i1e = i0e + P.r
        end
        # compute block i+1
        mu = shift - σ[ii]
        nu = β[ii+1] * (1 - shift/ξ[ii])
        if !P.is_low_rank || ii != P.p
            w[i1b:i1e] = mu/nu * w[i0b:i0e] + 1/nu * Bw[i1b:i1e]
        else
            w[i1b:i1e] = mu/nu * P.UU'*w[i0b:i0e] + 1/nu * Bw[i1b:i1e]
        end
        # range of next block i
        i0b = i1b
        i0e = i1e
    end

    return w
end

"True for complex points `z` inside polygonal set `Σ`."
function in_Σ(z::AbstractVector{CT}, Σ::AbstractVector{CT}, tol::T) where {T<:Real, CT<:Complex{T}}
    if length(Σ) == 2 && isreal(Σ)
        realΣ = real([Σ[1]; Σ[1]; Σ[2]; Σ[2]])
        imagΣ = [-tol; tol; tol; -tol]
    else
        realΣ = real(Σ)
        imagΣ = imag(Σ)
    end
    return map(p -> inpolygon(real(p), imag(p), realΣ, imagΣ), z)
end

function resize_matrix(A, rows, cols)
    resized = zeros(eltype(A), rows, cols)
    resized[1:size(A, 1), 1:size(A, 2)] = A
    return resized
end

struct NleigsSolutionDetails{T<:Real, CT<:Complex{T}}
    "matrix of Ritz values in each iteration"
    Lam::AbstractMatrix{CT}

    "matrix of residuals in each iteraion"
    Res::AbstractMatrix{T}

    "vector of interpolation nodes"
    σ::AbstractVector{CT}

    "vector of poles"
    ξ::AbstractVector{T}

    "vector of scaling parameters"
    β::AbstractVector{T}

    "vector of norms of generalized divided differences (in function handle
    case) or maximum of absolute values of scalar divided differences in
    each iteration (in matrix function case)"
    nrmD::AbstractVector{T}

    "number of iterations until linearization converged"
    kconv::Int
end

NleigsSolutionDetails{T,CT}() where {T<:Real, CT<:Complex{T}} = NleigsSolutionDetails(
    Matrix{CT}(undef, 0, 0), Matrix{T}(undef, 0, 0), Vector{CT}(),
    Vector{T}(), Vector{T}(), Vector{T}(), 0)


text:
"""
Split an iterable evenly into N chunks, which will be returned.
"""
function split_evenly(itr, N)
    count     = length(itr) ÷ N
    remainder = length(itr) % N

    map(0:N-1) do i
        # The first remainder chunks get count + 1 tasks
        if i < remainder
            start = 1 + i * (count + 1)
            stop  = start + count
        else
            start = 1 + i * count + remainder
            stop  = start + count - 1
        end
        itr[start:stop]
    end
end

ext: jl


text:
using Flux: onecold, onehotbatch
using Flux.Losses: logitcrossentropy
using Flux.Data: DataLoader
include("../preprocessing/ohe_cat_features.jl")

function loss_and_accuracy(data_loader::DataLoader, model, device;
    set=""::String, verbose=true)
    acc = 0
    ls = 0.0f0
    num = 0
    testmode!(model, true)
    for (y, X_cont, X_cat...) in data_loader
        y = y |> device
        X_cat = X_cat|> device
        X_cont = X_cont |> device
        pred = model(X_cat, X_cont)
        ls += logitcrossentropy(pred, y, agg=sum)
        acc += sum(onecold(cpu(pred)) .== onecold(cpu(y)))
        num +=  size(y, 2)
    end
    testmode!(model, false)

    ls = ls / num
    acc = acc / num
    if verbose
        println("  $set loss = $ls, $set accuracy = $acc")
    end
    return ls, acc
end


text:
@testset "1578.minimum-deletion-cost-to-avoid-repeating-letters.jl" begin
    @test min_cost_1578("abaac", [1, 2, 3, 4, 5]) == 3
    @test min_cost_1578("abc", [1, 2, 3]) == 0
    @test min_cost_1578("aabaa", [1, 2, 3, 4, 1]) == 2
    @test min_cost_1578("", Int[]) == 0
end


text:
using Documenter
using ForwardDiff
using GalacticOptim
using JETTest
using MeasureTheory
using Optim
using PointProcesses
using Random
using Test

Random.seed!(63)

include("doctests.jl")

include("history.jl")

include("markov.jl")

include("point_processes.jl")

include("models.jl")

include("hmm.jl")

include("utils.jl")


text:
using ClusterComplex
import ComputationalHomology: ComputationalHomology, filtration
import Clustering
import LMCLUS
import TOML
import Statistics: median, mean
import Random

const SEED = 921369837
Random.seed!(SEED)
prngs = [Random.MersenneTwister(SEED)]

KMeansPart = (Y, k)->Clustering.kmeans(collect(Y), k)
LMCLUSPart = (Y, bb)->let p = LMCLUS.Parameters(size(Y, 1))
           p.sampling_heuristic=2
           p.best_bound = abs(bb)
           p.basis_alignment = true
           p.bounded_cluster = true
           p.min_cluster_size = 30
           p.sampling_factor = 0.1
           res = LMCLUS.LMCLUSResult(LMCLUS.lmclus(Y, p, prngs)...)
           if bb < 0
                dfun = (X,m)  -> LMCLUS.distance_to_manifold(X, mean(m), LMCLUS.projection(m))
                efun = (X,ms) -> LMCLUS.MDL.calculate(LMCLUS.MDL.OptimalQuant, ms, X, 32, 16)
                res = LMCLUS.refine(res, Y, dfun, efun, bounds=true)
           end
           res
end

constractions = Dict(
    "vietorisrips" => (X, r)->ComputationalHomology.vietorisrips(X, r, maxoutdim=size(X,1)-1),
    "witness" => (X, r, l, hom=size(X,1)-1)->ComputationalHomology.witness(X, l, r, maxoutdim=hom),
    "lmcdist_km" => (X, r, k, hom=size(X,1)-1)->clustercomplex(X, KMeansPart(X, k), r, maxoutdim=hom, expansion = :inductive, method=:subspacemahalonobis),
    "lmcdist_lmc" => (X, r, bb, hom=size(X,1)-1)->clustercomplex(X, LMCLUSPart(X, bb), r, maxoutdim=hom, expansion = :inductive, method=:subspacemahalonobis),
    "mahalonobis_km" => (X, r, k, hom=size(X,1)-1)->clustercomplex(X, KMeansPart(X, k), r, maxoutdim=hom, expansion = :inductive),
    "mahalonobis_lmc" => (X, r, bb, hom=size(X,1)-1)->clustercomplex(X, LMCLUSPart(X, bb), r, maxoutdim=hom, expansion = :inductive),
    "mahalonobis_km_nu1" => (X, r, k, hom=size(X,1)-1)->clustercomplex(X, KMeansPart(X, k), r, maxoutdim=hom, expansion = :inductive, ν=1),
    "mahalonobis_lmc_nu1" => (X, r, bb, hom=size(X,1)-1)->clustercomplex(X, LMCLUSPart(X, bb), r, maxoutdim=hom, expansion = :inductive,ν=1),
)

function prepareparams(params, pnames = ["r", "l", "bb", "k", "hom"])
    res = Any[]
    for k in pnames
        haskey(params, k) && push!(res, params[k])
    end
    return tuple(res...)
end

function experiment(datasets, mthds; N=100)
    constraction_params = TOML.parsefile("experiments.toml")
    results = Dict(ds => Dict(m => Any[] for m in mthds) for ds in datasets)

    for ds in datasets
        X, H, L = ClusterComplex.dataset(ds)
        println("$ds: $H")

        for m in mthds
            print("Run $m:")
            n = N
            if m == "vietorisrips"
                n = 1
                if ds == "OIP300" || ds == "OIP15"
                    n = 0
                end
            end
            rdom = zeros(n)
            clen = zeros(Int, n)
            psrt = zeros(n)
            for i in 1:n
                divs = get(constraction_params[ds][m], "divs", Inf)
                try
                    res = constractions[m](X, prepareparams(constraction_params[ds][m])...)
                    flt = filtration(res[1:2]...; divisions=divs)
                    RD, R0, _ = ClusterComplex.dominance(flt, H)
                    rdom[i] = RD
                    psrt[i] = R0
                    if !isnan(RD)
                        cplx = complex(flt, R0)
                        clen[i] = length(cplx)
                    end
                    print(isnan(RD) ? "x" : ".")
                catch err
                    print("$i")
                end
            end
            succ = n == 1 ? [1] : findall(rdom .> 0)
            push!(results[ds][m], 100*length(succ)/N)
            push!(results[ds][m], length(succ) > 0 ? median(rdom[succ]) : 0.0)
            push!(results[ds][m], length(succ) > 0 ? mean(rdom[succ]) : 0.0)
            push!(results[ds][m], length(succ) > 0 ? median(clen[succ]) : 0.0)
            push!(results[ds][m], length(succ) > 0 ? mean(clen[succ]) : 0.0)
            push!(results[ds][m], length(succ) > 0 ? mean(psrt[succ]) : 0.0)
            println("\nmedian relative dominance: $(results[ds][m][2]) [$(results[ds][m][1])%]")
        end
    end

    return results
end

datasets = ["Sphere","TwoMoons","Circles","OIP300", "OIP15"]
# datasets = ["OIP300", "OIP15"]
# datasets = ["OIP15"]

mthds = ["witness","lmcdist_lmc","mahalonobis_lmc","lmcdist_km","mahalonobis_km","vietorisrips"]
# mthds = ["lmcdist_km","mahalonobis_km"]
# mthds = ["lmcdist_lmc","mahalonobis_lmc"]
# mthds = ["mahalonobis_lmc","mahalonobis_km"]
# mthds = ["witness", "vietorisrips"]
# mthds = ["witness","mahalonobis_lmc","mahalonobis_km","vietorisrips"]
# mthds = ["witness"]
# mthds = ["vietorisrips"]

res = experiment(datasets, mthds)
open(io->TOML.print(io, res), "results.toml", "w")
mv("results.toml", "results-$(time_ns()).toml")
# success rate, median dominance, mean dominance, median complex size, mean complex size, profile start F-value


text:
using Base.Test
using Compat
using FixedPointNumbers

import Images: Image
import VideoIO

testdir = dirname(@__FILE__)
videodir = joinpath(testdir, "..", "videos")

VideoIO.TestVideos.available()
VideoIO.TestVideos.download_all()

swapext(f, new_ext) = "$(splitext(f)[1])$new_ext"

println(STDERR, "Testing file reading...")

if !isdefined(Main, :UFixed8)
    UFixed8 = Ufixed8
end

function notblank(img)
    all(Images.green(img) .== 0x00uf8) || all(Images.blue(img) .== 0x00uf8) || all(Images.red(img) .== 0x00uf8) || maximum(reinterpret(UFixed8, img)) < 0xcfuf8
end

if isdefined(Images, :load)
    imload = Images.load
else
    imload = Images.imread
end

for name in VideoIO.TestVideos.names()
    is_apple() && startswith(name, "crescent") && continue
    println(STDERR, "   Testing $name...")

    first_frame_file = joinpath(testdir, swapext(name, ".png"))
    first_frame = imload(first_frame_file) # comment line when creating png files

    f = VideoIO.testvideo(name)
    v = VideoIO.openvideo(f)

    img = read(v, Image)

    # Find the first non-trivial image
    while notblank(img)
        read!(v, img)
    end

    #imwrite(img, first_frame_file)        # uncomment line when creating png files

    @test img == first_frame               # comment line when creating png files

    while !eof(v)
        read!(v, img)
    end

    # read first frames again, and compare
    seekstart(v)

    read!(v, img)

    while notblank(img)
        read!(v, img)
    end

    @test img == first_frame

    close(v)
end

println(STDERR, "Testing IO reading...")
for name in VideoIO.TestVideos.names()
    is_apple() && startswith(name, "crescent") && continue
    # TODO: fix me?
    (startswith(name, "ladybird") || startswith(name, "NPS")) && continue

    println(STDERR, "   Testing $name...")
    first_frame_file = joinpath(testdir, swapext(name, ".png"))
    first_frame = imload(first_frame_file) # comment line when creating png files

    filename = joinpath(videodir, name)
    v = VideoIO.openvideo(open(filename))

    img = read(v, Image)

    # Find the first non-trivial image
    while notblank(img)
        read!(v, img)
    end

    #imwrite(img, first_frame_file)        # uncomment line when creating png files

    @test img == first_frame               # comment line when creating png files

    while !eof(v)
        read!(v, img)
    end
end

VideoIO.testvideo("ladybird") # coverage testing
@test_throws ErrorException VideoIO.testvideo("rickroll")
@test_throws ErrorException VideoIO.testvideo("")

#VideoIO.TestVideos.remove_all()


text:
module SimpleTensorNetworks

using Requires
using LinearAlgebra

include("tensors.jl")
include("tensorcontract.jl")
include("simplify.jl")
include("contractionorder/contractionorder.jl")

function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
    @require Viznet = "52a3aca4-6234-47fd-b74a-806bdf78ede9" include("viz.jl")
end

end


text:
mutable struct RelPerms{T <: Float64}
    swr::T
    sor::T
    krw0::T
    kro0::T
    nw::T
    no::T

    function RelPerms(;swr::T, sor::T, krw0::T, kro0::T, nw::T, no::T) where T <: Float64
        new{Float64}(swr, sor, krw0, kro0, nw, no)
    end
end

struct WaterFlooding{T <: Real}
    S̃::T
    Si::T
    Sj::T
    kr::RelPerms
    μw::T
    μo::T
end

struct PolymerFlooding{T <: Real}
    S̃::T
    Si::T
    Sj::T
    kr1::RelPerms
    kr2::RelPerms
    μw::T
    μo::T
    Sb::T
    D::T
    Vb::T
    VΔc::T
end

struct Tracer
    phase::Symbol
    sc::Float64
    v::Float64
end


text:
module TestIndexing

using Test, DataFrames

@testset "getindex DataFrame" begin
    df = DataFrame(a=1:3, b=4:6, c=7:9)

    @test df[1] == [1, 2, 3]
    @test df[1] === eachcol(df)[1]
    @test df[1:2] == DataFrame(a=1:3, b=4:6)
    @test eachcol(df, false)[1] === df[1]
    @test eachcol(view(df,1:2), false)[1] == eachcol(df, false)[1]
    @test eachcol(df[1:2], false)[1] == eachcol(df, false)[1]
    @test eachcol(df[1:2], false)[1] !== eachcol(df, false)[1]
    @test df[:] == df
    @test df[:] !== df
    @test eachcol(view(df, :), false)[1] == eachcol(df, false)[1]
    @test eachcol(df[:], false)[1] == eachcol(df, false)[1]
    @test eachcol(df[:], false)[1] !== eachcol(df, false)[1]
    @test eachcol(df)[1] === last(eachcol(df, true)[1])
    @test eachcol(df)[1] === last(eachcol(df, true)[1])

    @test df[1, 1] == 1
    @test df[1, 1:2] isa DataFrameRow
    @test copy(df[1, 1:2]) == (a=1, b=4)
    @test df[1, :] isa DataFrameRow
    @test copy(df[1, :]) == (a=1, b=4, c=7)
    @test parent(df[1, :]) === df

    @test df[1:2, 1] == [1, 2]
    @test df[1:2, 1:2] == DataFrame(a=1:2, b=4:5)
    @test df[1:2, :] == DataFrame(a=1:2, b=4:5, c=7:8)

    @test df[:, 1] == [1, 2, 3]
    @test df[:, 1] !== df[1]
    @test df[:, 1:2] == DataFrame(a=1:3, b=4:6)
    @test eachcol(df[:, 1:2])[1] !== df[1]
    @test df[:, :] == df
    @test eachcol(df[:, :])[1] !== df[1]
end

@testset "getindex df[col] and df[cols]" begin
    x = [1, 2, 3]
    df = DataFrame(x=x, copycols=false)
    @test df.x === x
    @test df[:x] === x
    @test df[[:x]].x !== x
    @test df[:].x !== x
end

@testset "view DataFrame" begin
    df = DataFrame(a=1:3, b=4:6, c=7:9)

    @test view(df, 1) == [1, 2, 3]
    @test view(df, 1) isa SubArray
    @test view(df, 1:2) isa SubDataFrame
    @test view(df, 1:2) == df[1:2]
    @test view(df, :) isa SubDataFrame
    @test view(df, :) == df
    @test parent(view(df, :)) === df

    @test view(df, 1, 1) isa SubArray
    @test view(df, 1, 1)[] == 1
    @test view(df, 1, 1:2) isa DataFrameRow
    @test copy(view(df, 1, 1:2)) == (a=1, b=4)
    @test view(df, 1, :) isa DataFrameRow
    @test copy(view(df, 1, :)) == (a=1, b=4, c=7)
    @test parent(view(df, 1, :)) === df

    @test view(df, 1:2, 1) == [1, 2]
    @test view(df, 1:2, 1) isa SubArray
    @test view(df, 1:2, 1:2) isa SubDataFrame
    @test view(df, 1:2, 1:2) == DataFrame(a=1:2, b=4:5)
    @test view(df, 1:2, :) isa SubDataFrame
    @test view(df, 1:2, :) == df[1:2, :]
    @test parent(view(df, 1:2, :)) === df

    @test view(df, :, 1) == [1, 2, 3]
    @test view(df, :, 1) isa SubArray
    @test view(df, :, 1:2) isa SubDataFrame
    @test view(df, :, 1:2) == DataFrame(a=1:3, b=4:6)
    @test view(df, :, :) isa SubDataFrame
    @test view(df, :, :) == df[:, :]
    @test parent(view(df, :, :)) === df
end

@testset "getindex SubDataFrame" begin
    df = DataFrame(x=-1:3, a=0:4, b=3:7, c=6:10, d=9:13)
    sdf = view(df, 2:4, 2:4)

    @test sdf[1] == [1, 2, 3]
    @test sdf[1] isa SubArray
    @test sdf[1:2] == DataFrame(a=1:3, b=4:6)
    @test sdf[1:2] isa SubDataFrame
    @test sdf[:] == df[2:4, 2:4]
    @test sdf[:] isa SubDataFrame
    @test parent(sdf[:]) === parent(sdf)

    @test sdf[1, 1] == 1
    @test sdf[1, 1:2] isa DataFrameRow
    @test copy(sdf[1, 1:2]) == (a=1, b=4)
    @test sdf[1, :] isa DataFrameRow
    @test copy(sdf[1, :]) == (a=1, b=4, c=7)
    @test parent(sdf[1, :]) === parent(sdf)

    @test sdf[1:2, 1] == [1, 2]
    @test sdf[1:2, 1] isa Vector
    @test sdf[1:2, 1:2] == DataFrame(a=1:2, b=4:5)
    @test sdf[1:2, 1:2] isa DataFrame
    @test sdf[1:2, :] == DataFrame(a=1:2, b=4:5, c=7:8)
    @test sdf[1:2, :] isa DataFrame

    @test sdf[:, 1] == [1, 2, 3]
    @test sdf[:, 1] isa Vector
    @test sdf[:, 1] !== df[1]
    @test sdf[:, 1:2] == DataFrame(a=1:3, b=4:6)
    @test sdf[:, 1:2] isa DataFrame
    @test sdf[:, :] == df[2:4, 2:4]
    @test sdf[:, :] isa DataFrame
end

@testset "view SubDataFrame" begin
    df = DataFrame(x=-1:3, a=0:4, b=3:7, c=6:10, d=9:13)
    sdf = view(df, 2:4, 2:4)

    @test view(sdf, 1) == [1, 2, 3]
    @test view(sdf, 1) isa SubArray
    @test view(sdf, 1:2) isa SubDataFrame
    @test view(sdf, 2:3) == df[2:4, 3:4]
    @test view(sdf, :) isa SubDataFrame
    @test view(sdf, :) == df[2:4, 2:4]
    @test parent(view(sdf, :)) == parent(sdf)

    @test view(sdf, 1, 1) isa SubArray
    @test view(sdf, 1, 1)[] == 1
    @test view(sdf, 1, 1:2) isa DataFrameRow
    @test copy(view(sdf, 1, 1:2)) == (a=1, b=4)
    @test view(sdf, 1, :) isa DataFrameRow
    @test copy(view(sdf, 1, :)) == (a=1, b=4, c=7)
    @test parent(view(sdf, 1, :)) === parent(sdf)

    @test view(sdf, 1:2, 1) == [1, 2]
    @test view(sdf, 1:2, 1) isa SubArray
    @test view(sdf, 1:2, 1:2) isa SubDataFrame
    @test view(sdf, 1:2, 1:2) == DataFrame(a=1:2, b=4:5)
    @test view(sdf, 1:2, :) isa SubDataFrame
    @test view(sdf, 1:2, :) == df[2:3, 2:4]
    @test parent(view(sdf, 1:2, :)) === parent(sdf)

    @test view(sdf, :, 1) == [1, 2, 3]
    @test view(sdf, :, 1) isa SubArray
    @test view(sdf, :, 1:2) isa SubDataFrame
    @test view(sdf, :, 1:2) == DataFrame(a=1:3, b=4:6)
    @test view(sdf, :, :) isa SubDataFrame
    @test parent(view(sdf, :, :)) === parent(sdf)
    @test view(sdf, :, :) == df[2:4, 2:4]
end

@testset "getindex DataFrameRow" begin
    df = DataFrame(a=1:4, b=4:7, c=7:10)
    dfr = df[1, :]

    @test dfr[1] == 1
    @test dfr[1:2] isa DataFrameRow
    @test copy(dfr[1:2]) == (a=1, b=4)
    @test dfr[:] isa DataFrameRow
    @test copy(dfr[:]) == (a=1, b=4, c=7)
    @test parent(dfr[:]) === df
end

@testset "view DataFrameRow" begin
    df = DataFrame(a=1:4, b=4:7, c=7:10)
    dfr = df[1, :]

    @test view(dfr, 1)[] == 1
    @test view(dfr, 1) isa SubArray
    @test view(dfr, 1:2) isa DataFrameRow
    @test copy(dfr[1:2]) == (a=1, b=4)
    @test dfr[:] isa DataFrameRow
    @test copy(dfr[:]) == (a=1, b=4, c=7)
    @test parent(dfr[:]) === df
end

end # module


text:
"""
**Bayesian SVM**

The [Bayesian SVM](https://arxiv.org/abs/1707.05532) is a Bayesian interpretation of the classical SVM.
``p(y|f) \\propto \\exp\\left(2\\max(1-yf,0)\\right)``

```julia
BayesianSVM()
```
---
For the analytic version of the likelihood, it is augmented via:
```math
p(y|f,\\omega) = \\frac{1}{\\sqrt{2\\pi\\omega}}\\exp\\left(-\\frac{1}{2}\\frac{(1+\\omega-yf)^2}{\\omega}\\right)
```
where ``\\omega\\sim 1_{[0,\\infty]}`` has an improper prior (his posterior is however has a valid distribution (Generalized Inverse Gaussian)). For reference [see this paper](http://ecmlpkdd2017.ijs.si/papers/paperID502.pdf)
"""
struct BayesianSVM{T<:Real} <: ClassificationLikelihood{T}
    ω::AbstractVector{AbstractVector{T}}
    θ::AbstractVector{AbstractVector{T}}
    function BayesianSVM{T}() where {T<:Real}
        new{T}()
    end
    function BayesianSVM{T}(ω::AbstractVector{<:AbstractVector{<:Real}},θ::AbstractVector{<:AbstractVector{<:Real}}) where {T<:Real}
        new{T}(ω,θ)
    end
end

function BayesianSVM()
    BayesianSVM{Float64}()
end

function init_likelihood(likelihood::BayesianSVM{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Integer,nFeatures::Integer) where T
    BayesianSVM{T}([abs.(rand(T,nSamplesUsed)) for _ in 1:nLatent],[zeros(T,nSamplesUsed) for _ in 1:nLatent])
end
function pdf(l::BayesianSVM,y::Real,f::Real)
    svmlikelihood(y*f)
end

function Base.show(io::IO,model::BayesianSVM{T}) where T
    print(io,"Bayesian SVM")
end

"""Return likelihood equivalent to SVM hinge loss"""
function svmlikelihood(f::Real)
    pos = svmpseudolikelihood(f)
    return pos./(pos.+svmpseudolikelihood(-f))
end

"""Return the pseudo likelihood of the SVM hinge loss"""
function svmpseudolikelihood(f::Real)
    return exp(-2.0*max.(1.0-f,0))
end

function compute_proba(l::BayesianSVM{T},μ::Vector{T},σ²::Vector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    sig_pred = zeros(T,N)
    for i in 1:N
        if σ²[i] <= 0.0
            pred[i] = svmlikelihood(μ[i])
        else
            nodes = pred_nodes.*sqrt2.*sqrt.(σ²[i]).+μ[i]
            pred[i] =  dot(pred_weights,svmlikelihood.(nodes))
            sig_pred[i] = dot(pred_weights,svmlikelihood.(nodes).^2)-pred[i]^2
        end
    end
    return pred, sig_pred
end

###############################################################################

function local_updates!(model::VGP{T,<:BayesianSVM,<:AnalyticVI}) where {T}
    model.likelihood.ω .= broadcast((μ,Σ,y)->abs2.(one(T) .- y.*μ) + Σ ,model.μ,diag.(model.Σ),model.y)
    model.likelihood.θ .= broadcast(b->one(T)./sqrt.(b),model.likelihood.ω)
end

function local_updates!(model::SVGP{T,<:BayesianSVM,<:AnalyticVI}) where {T}
    model.likelihood.ω .= broadcast((κ,μ,Σ,y,K̃)->abs2.(one(T) .- y.*(κ*μ)) + opt_diag(κ*Σ,κ) + K̃,model.κ,model.μ,model.Σ,model.inference.y,model.K̃)
    model.likelihood.θ .= broadcast(b->one(T)./sqrt.(b),model.likelihood.ω)
end

@inline ∇E_μ(model::AbstractGP{T,<:BayesianSVM, <:GibbsorVI}) where {T} = broadcast((y,θ)->y.*(θ.+one(T)),model.inference.y,model.likelihood.θ)
@inline ∇E_μ(model::AbstractGP{T,<:BayesianSVM, <:GibbsorVI},i::Int) where {T} = model.inference.y[i].*(model.likelihood.θ[i].+one(T))

@inline ∇E_Σ(model::AbstractGP{T,<:BayesianSVM,<:GibbsorVI}) where {T} = 0.5.*model.likelihood.θ
@inline ∇E_Σ(model::AbstractGP{T,<:BayesianSVM,<:GibbsorVI},i::Int) where {T} = 0.5.*model.likelihood.θ[i]

function ELBO(model::AbstractGP{T,<:BayesianSVM,<:AnalyticVI}) where {T}
    return expecLogLikelihood(model) - GaussianKL(model) - GIGEntropy(model)
end

function expecLogLikelihood(model::VGP{T,<:BayesianSVM,<:AnalyticVI}) where {T}
    tot = -model.nLatent*(0.5*model.nSample*logtwo)
    tot += sum(broadcast((μ,y,θ,Σ)->(sum(μ.*y)-0.5*dot(θ,Σ+abs2.(one(T).-y.*μ))),
                        model.μ,model.y,model.likelihood.θ,diag.(model.Σ)))
    return tot
end

function expecLogLikelihood(model::SVGP{T,<:BayesianSVM,<:AnalyticVI}) where {T}
    tot = -model.nLatent*(0.5*model.nSample*logtwo)
    tot += sum(broadcast((κμ,y,θ,κΣκ,K̃)->(sum(κμ.*y)-0.5*dot(θ,K̃+κΣκ+abs2.(one(T).-y.*κμ))),
                        model.κ.*model.μ,model.inference.y,model.likelihood.θ,opt_diag.(model.κ.*model.Σ,model.κ),model.K̃))
    return model.inference.ρ*tot
end


text:
function doFDR(pvalues,FDR=0.1)
    tri = sortperm(pvalues)
    m = length(pvalues)
    iBH = findall([0;pvalues[tri]].<= collect(0:m)/m*FDR)
    iBH = iBH[length(iBH)] # the last
    if (iBH<=m)
        H0 = tri[iBH:m]
    else H0 = Int64[]
    end
    H0
end

getChildren(pixIdx,infoDen,dataProj,dataVar,thrs,toCluster) = getChildren(CPU1(),pixIdx,infoDen,dataProj,dataVar,thrs,toCluster)

function getChildren(resource::CPU1,pixIdx,infoDen,dataProj,dataVar,thrs,toCluster)

    if (length(infoDen[pixIdx]) == 0)
        return nothing
    end
    ## the projection of the original dynamic of pix.idx
    pixOrg = dataProj[:,pixIdx]
    ## the projection of the denoised dynamic of pixIdx
    pixDen = infoDen[pixIdx]["Px"]
    ## the length of pix.idx neighborhood
    pixLen = length(infoDen[pixIdx]["Vx"])

    ## initialize the children of pix.idx to its neighbors
    children = getActive(infoDen[pixIdx]["Vx"],toCluster)

    i = 1
    while (i < length(children))
        i = i+1
        infoChild = infoDen[children[i]]
        varChild = dataVar[children[i]]
        ## test if denoised version are time coherent
        if multitestH0(infoChild["Px"]-pixOrg,varChild/length(infoChild["Vx"]) + dataVar[pixIdx],thrs)[1]
            ## Time coherent : we add non existing neighbors to children
            children = [children;setdiff(getActive(infoChild["Vx"],toCluster),children)]
        else
            ## not time coherent : remove from children
            deleteat!(children,i)
            i = i-1
        end
    end
    children::Array{Int64}
end

function getActive(indexes,toCluster)
    indexes[toCluster[indexes].==1]
end

checkClusterNew(center,centClusts,clusterDict,dataVar,quant,iter;pvalueMax=false,useFdr=true)=checkClusterNew(CPU1(),center,centClusts,clusterDict,dataVar,quant,iter,pvalueMax=pvalueMax,useFdr=useFdr)

function checkClusterNew(resource::CPU1,center,centClusts,clusterDict,dataVar,quant,iter;pvalueMax=false,useFdr=true)
    ## test new potential center against the centers from already build clusters
    ## if pvalue.max=T returns the index of the cluster giving the largest pvalue
    ## if pvalue.max=F returns the indexes of the clusters with pvalue corresponding to H0 (default)
    ## if use.fdr=TRUE (default) correct multiplicity with False Discovery Rate otherwise use Bonferroni correction
    sizeto = length(centClusts)

    if length(clusterDict)==0
        if !pvalueMax
            return([])
        else
            return(false)
        end
    end

    centVar = mean(dataVar[centClusts])/sizeto
    clusterCenters = Array{Float64}(undef,(length(center),length(clusterDict)))
    clusterVar = Array{Float64}(undef,length(clusterDict))
    clusterSizes = Array{Float64}(undef,length(clusterDict))
    for (i,x) in enumerate(clusterDict)
        clusterCenters[:,i] = x["center"]
        clusterSizes[i] = length(x["cluster"])
        clusterVar[i] = mean(x["cluster"])/length(x["cluster"])
    end

    ## test thresholds with Bonferroni correction adapted to both partition number and number of existing clusters
    if (pvalueMax || useFdr)
        thrs = nothing
    else
        thrs = [quantile(Chisq(i),1-quant/length(clusterDict)) for i=exp2(0:iter)]
    end
    tmp = multitestH0(clusterCenters.-center[:],clusterVar .+ centVar,thrs)
    ## return pvalues if thrs==NULL and booleans otherwise

    if pvalueMax
        return(argmax(tmp))
    elseif useFdr
        return(doFDR(tmp))
    else
        return(find(tmp))
    end
end

robustMean(children,dataProj,dataVar)=robustMean(CPU1(),children,dataProj,dataVar)

function robustMean(resource::CPU1,children,dataProj,dataVar)
    ## Build a robust 1-mean for a children list
    ## Initialize with the dynamic average over all children
    childProj = dataProj[:,children]
    cluster = children
    center = Inf
    i = 0
    iMax = 1000
    cont = true
    ## repeat until center does not change anymore
    while (cont && i<iMax)
        i = i+1
        ## compute the new cluster center
        centerNew = mean(dataProj[:,cluster],dims=2)
        ## test coherence between cluster center and children dynamics
        pvalues = multitestH0(childProj.-centerNew,dataVar[children].+mean(dataVar[cluster])/length(cluster),nothing)
        #println(pvalues)
        ## the new cluster
        cluster = children[doFDR(pvalues)]
        ## check if the cluster center has changed or not
        cont = (length(cluster)>0) && (mean((centerNew .- center).^2) > 0.0001)
        ## update
        center = centerNew
    end
    outliers=setdiff(children,cluster)

    if i==iMax
        println("stop RobustMean without stabilization")
    end

    Dict("cluster" => cluster,"center" => center,"outliers" => outliers,"OK" => (i==iMax))
end


text:
using Test, NetCDF, Random

"Function that returns a valid file name of a non-existing file for all OSs"
function tempname2()
    f, h = mktemp()
    close(h)
    rm(f)
    f
end

include("intermediate.jl")
include("high.jl")
include("array-like.jl")
include("chunks.jl")


text:
module Percolation

using Gtk, Gtk.ShortNames, StatsBase, Plots, Random
gr()
import Plots: heatmap, heatmap!
import Images: label_components, label_components!

include("latticetype.jl")
include("connected.jl")
include("clusterproperties.jl")
include("forestfire.jl")
include("visualize.jl")
include("ui.jl")

end # module


text:
function KS_solve_Emin_sd( pw::PWGrid, V_ionic, Focc, Nstates::Int;
                         psi0=nothing, Potentials0 = nothing, E_NN = nothing,
                         α_t = 3e-5, NiterMax=1000, verbose=false )
    Ns = pw.Ns
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx

    if psi0 == nothing
        srand(2222)
        psi = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)
        psi = ortho_gram_schmidt(psi)
    else
        psi = copy(psi0)
    end

    if Potentials0 == nothing
        Potentials = PotentialsT( V_ionic, zeros(Npoints), zeros(Npoints) )
        rho = calc_rho( pw, Focc, psi )
        Potentials.Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
        Potentials.XC = excVWN( rho ) + rho .* excpVWN( rho )
    else
        Potentials = PotentialsT( Potentials0.Ionic,
                                  Potentials0.Hartree,
                                  Potentials0.XC )
    end

    Etot = 0.0
    Etot_old = 0.0
    if( E_NN != nothing )
        Energies = EnergiesT( 0.0, 0.0, 0.0, 0.0, 0.0, E_NN )
    else
        Energies = EnergiesT( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
    end

    for iter = 1:NiterMax
        psi = psi - α_t*calc_grad( pw, Potentials, Focc, psi )

        psi = ortho_gram_schmidt(psi)

        # Update potentials
        rho = calc_rho( pw, Focc, psi)
        Potentials.Hartree = real( G_to_R( Ns, Poisson_solve( pw, rho ) ) )
        Potentials.XC = excVWN( rho ) + rho .* excpVWN( rho )

        Energies = calc_Energies( pw, Potentials, Focc, psi, Energies.NN )
        Etot = Energies.Total

        conv = abs(Etot-Etot_old)
        @printf("Emin SD: (s, step, conv) %8d %18.10f %18.10e\n", iter, Etot, conv )
        if conv < 1e-6
          print("Convergence achieved\n")
          break
        end
        Etot_old = Etot
    end
    return psi, Energies, Potentials
end


text:
function print_build_tarballs(io::IO, state::WizardState)
    urlfiles = zip(state.source_urls, state.source_files)

    sources_strings = map(urlfiles) do x
        url_string = repr(x[1])
        if endswith(x[1], ".git")
            "GitSource($(url_string), $(repr(x[2].hash)))"
        elseif any(endswith(x[1], ext) for ext in archive_extensions)
            "ArchiveSource($(url_string), $(repr(x[2].hash)))"
        else
            "FileSource($(url_string), $(repr(x[2].hash)))"
        end
    end
    if !isempty(state.patches)
        push!(sources_strings, "DirectorySource(\"./bundled\")")
    end
    sources_string = strip(join(sources_strings, ",\n    "))
    if Set(state.platforms) == Set(supported_platforms())
        platforms_string = "supported_platforms()"
    else
        platforms_string = """
        [
            $(strip(join(state.platforms,",\n    ")))
        ]
        """
    end

    stuff = collect(zip(state.files, state.file_kinds, state.file_varnames))
    products_string = join(map(stuff) do x
        file, kind, varname = x

        if kind == :executable
            # It's very common that executable products are in `bin`,
            # so we special-case these to just print the basename:
            dir_path = dirname(file)
            if isempty(dir_path) || dirname(file) == "bin"
                return "ExecutableProduct($(repr(basename(file))), $(repr(varname)))"
            else
                return "ExecutableProduct($(repr(basename(file))), $(repr(varname)), $(repr(dir_path)))"
            end
        elseif kind == :library
            # Normalize library names, since they are almost always in `lib`/`bin`
            return "LibraryProduct($(repr(normalize_name(file))), $(repr(varname)))"
        else
            return "FileProduct($(repr(file)), $(repr(varname)))"
        end
    end,",\n    ")

    if length(state.dependencies) >= 1
        psrepr(ps) = "\n    Dependency(PackageSpec(name=\"$(getname(ps))\", uuid=\"$(getpkg(ps).uuid)\"))"
        dependencies_string = "[" * join(map(psrepr, state.dependencies)) * "\n]"
    else
        dependencies_string = "Dependency[\n]"
    end

    # Keyword arguments to `build_tarballs()`
    kwargs_vec = String[]
    if state.compilers != [:c] && length(state.compilers) >= 1
        push!(kwargs_vec, "compilers = [$(join(map(x -> ":$(x)", state.compilers), ", "))]")
    end
    # Default GCC version is the oldest one
    if state.preferred_gcc_version != getversion(available_gcc_builds[1])
        push!(kwargs_vec, "preferred_gcc_version = v\"$(state.preferred_gcc_version)\"")
    end
    # Default LLVM version is the latest one
    if state.preferred_llvm_version != getversion(available_llvm_builds[end])
        push!(kwargs_vec, "preferred_llvm_version = v\"$(state.preferred_llvm_version)\"")
    end
    kwargs = ""
    if length(kwargs_vec) >= 1
        kwargs = "; " * join(kwargs_vec, ", ")
    end

    print(io, """
    # Note that this script can accept some limited command-line arguments, run
    # `julia build_tarballs.jl --help` to see a usage message.
    using BinaryBuilder, Pkg

    name = $(repr(state.name))
    version = $(repr(state.version))

    # Collection of sources required to complete build
    sources = [
        $(sources_string)
    ]

    # Bash recipe for building across all platforms
    script = raw\"\"\"
    $(strip(state.history))
    \"\"\"

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = $(platforms_string)

    # The products that we will ensure are always built
    products = [
        $(products_string)
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = $(dependencies_string)

    # Build the tarballs, and possibly a `build.jl` as well.
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies$(kwargs))
    """)
end

"""
    yggdrasil_deploy(state::WizardState)

Write out a WizardState to a `build_tarballs.jl` in an `Yggdrasil` clone, then
open a pull request against `Yggdrasil`.
"""
function yggdrasil_deploy(state::WizardState, open_pr::Bool = true)
    btb = IOBuffer()
    print_build_tarballs(btb, state)
    seek(btb, 0)
    build_tarballs_content = String(read(btb))

    return yggdrasil_deploy(
        state.name,
        state.version,
        state.patches,
        build_tarballs_content;
        open_pr=open_pr
    )
end

function yggdrasil_deploy(name, version, patches, build_tarballs_content; open_pr::Bool=false, branch_name=nothing)
    # First, fork Yggdrasil (this just does nothing if it already exists)
    gh_auth = github_auth(;allow_anonymous=false)
    fork = GitHub.create_fork("JuliaPackaging/Yggdrasil"; auth=gh_auth)

    mktempdir() do tmp
        # Clone our bare Yggdrasil out to a temporary directory
        repo = LibGit2.clone(get_yggdrasil(), tmp)

        # Check out a unique branch name
        @info("Checking temporary Yggdrasil out to $(tmp)")
        if branch_name === nothing
            recipe_hash = bytes2hex(sha256(build_tarballs_content)[end-3:end])
            branch_name = "wizard/$(name)-v$(version)_$(recipe_hash)"
        end
        LibGit2.branch!(repo, branch_name)

        # Spit out the buildscript to the appropriate file:
        rel_bt_path = yggdrasil_build_tarballs_path(name)
        @info("Generating $(rel_bt_path)")
        output_path = joinpath(tmp, rel_bt_path)
        mkpath(dirname(output_path))
        open(output_path, "w") do io
            write(io, build_tarballs_content)
        end

        # If needed, add patches
        if !isempty(patches)
            rel_patches_path = joinpath(dirname(rel_bt_path), "bundled", "patches")
            mkpath(joinpath(tmp, rel_patches_path))
            for patch in patches
                patch_path = joinpath(rel_patches_path, patch.name)
                open(f->write(f, patch.patch), joinpath(tmp, patch_path), "w")
                LibGit2.add!(repo, patch_path)
            end
        end

        # Commit it and push it up to our fork
        @info("Committing and pushing to $(fork.full_name)#$(branch_name)...")
        LibGit2.add!(repo, rel_bt_path)
        LibGit2.commit(repo, "New Recipe: $(name) v$(version)")
        creds = LibGit2.UserPasswordCredential(
            dirname(fork.full_name),
            deepcopy(gh_auth.token)
        )
        try
            LibGit2.push(
                repo,
                refspecs=["+HEAD:refs/heads/$(branch_name)"],
                remoteurl="https://github.com/$(fork.full_name).git",
                credentials=creds,
                # This doesn't work :rage: instead we use `+HEAD:` at the beginning of our
                # refspec: https://github.com/JuliaLang/julia/issues/23057
                #force=true,
            )
        finally
            Base.shred!(creds)
        end

        if open_pr
            # Open a pull request against Yggdrasil
            @info("Opening a pull request against JuliaPackaging/Yggdrasil...")
            params = Dict(
                "base" => "master",
                "head" => "$(dirname(fork.full_name)):$(branch_name)",
                "maintainer_can_modify" => true,
                "title" => "Wizard recipe: $(name)-v$(version)",
                "body" => """
                This pull request contains a new build recipe I built using the BinaryBuilder.jl wizard:

                * Package name: $(name)
                * Version: v$(version)

                @staticfloat please review and merge.
                """
            )
            pr = create_or_update_pull_request("JuliaPackaging/Yggdrasil", params, auth=gh_auth)
            @info("Pull request created: $(pr.html_url)")
        else
            println("Open the pull request by going to: ")
            println("https://github.com/$(fork.full_name)/pull/new/$(HTTP.escapeuri(branch_name))?expand=1")
        end
    end
end

function step7(state::WizardState)
    print(state.outs, "Your build script was:\n\n\t")
    println(state.outs, replace(state.history, "\n" => "\n\t"))

    printstyled(state.outs, "\t\t\t# Step 7: Deployment\n\n", bold=true)

    terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
    deploy_select = request(terminal,
        "How should we deploy this build recipe?",
        RadioMenu([
            "Prepare a pull request against the community buildtree, Yggdrasil",
            "Write to a local file",
            "Print to stdout",
        ])
    )
    println(state.outs)

    if deploy_select == 1
        yggdrasil_select = request(terminal,
            "Would you like to actually open the pull request or just prepare it?",
            RadioMenu([
                "Go ahead, open it",
                "No, just prepare it and let me look at it first",
            ])
        )
        yggdrasil_deploy(state, yggdrasil_select == 1)
    elseif deploy_select == 2
        filename = nonempty_line_prompt(
            "filename",
            "Enter path to build_tarballs.jl to write to:";
            ins=state.ins,
            outs=state.outs,
        )
        mkpath(dirname(abspath(filename)))
        open(filename, "w") do io
            print_build_tarballs(io, state)
        end
    else
        println(state.outs, "Your generated build_tarballs.jl:")
        println(state.outs, "\n```")
        print_build_tarballs(state.outs, state)
        println(state.outs, "```")
    end
end


text:
using TensorBoardLogger #import the TensorBoardLogger package
using Logging #import Logging package

logger = TBLogger("audiologs", tb_append) #create tensorboard logger

################log audio example: cos wave ################
samplerate = 44100
x = collect(1:samplerate*2)
x = cos.(440*pi*x/samplerate)
#using logger interface
with_logger(logger) do
    @info "cos/loggerinterface" x = TBAudio(x, samplerate)
end
#using explicit function interface
log_audio(logger, "cos/explicitinterface", x, samplerate, step = 0)

################log scalar example: sin+cos wave################
y = collect(1:samplerate*2)
y = sin.(440*pi*y/samplerate)
z = x+y
#using logger interface
with_logger(logger) do
    @info "sin+cos/loggerinterface" z = TBAudio(z, samplerate)
end
#using explicit function interface
log_audio(logger, "sin+cos/explicitinterface", z, samplerate, step = 0)


text:
using Main.EntropyM: entropy

while true

    println("\nText to analyse [leave blank to break]: ")
    human_input = readline()
    human_input == "" ? break :
    println("\nAlphabet [leave blank for default]: ")
    human_alphabet = readline()
    human_alphabet == "" ? human_alphabet = "abcdefghijklmnopqrstuvwxyz" :
    println()

    println("\nEntropy: ", entropy(uppercase(human_input), uppercase(human_alphabet)))

end # while

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
module DebugRequest

import ..Layer, ..request
using ..IOExtras

const live_mode = true

include("IODebug.jl")

"""
    request(DebugLayer, ::IO, ::Request, body) -> HTTP.Response

Wrap the `IO` stream in an `IODebug` stream and print Message data.
"""
abstract type DebugLayer{Next <:Layer} <: Layer{Next} end
export DebugLayer

function request(::Type{DebugLayer{Next}}, io::IO, req, body; kw...) where Next

    @static if live_mode
        return request(Next, IODebug(io), req, body; kw...)
    else
        iod = IODebug(io)
        try
            return request(Next, iod, req, body; kw...)
        finally
            show_log(stdout, iod)
        end
    end
end

end # module DebugRequest


text:
# Author: Jeremy Werner

# Start in the same directory as kNN.jl
include("kNN.jl")

using RDatasets

iris = dataset("datasets", "iris")

# PetalLength and PetalWidth will be used for 2D example
x_all = hcat(iris[:PetalLength], iris[:PetalWidth])

subset_n = 120
x_train = x_all[1:subset_n, 1:2]
x_test = x_all[(subset_n + 1):150, 1:2]
y_in = array(iris[:Species][1:subset_n])

real_y = array(iris[:Species][(subset_n+1):(size(iris)[1])])
n1 = length(real_y)

# Fitting for k=2,6
test_predictions_2 = kNN(x_train, y_in, 2, x_test)
test_predictions_6 = kNN(x_train, y_in, 6, x_test)

# Classification rates
class_2 = sum(real_y .== test_predictions_2) / n1
class_6 = sum(real_y .== test_predictions_6) / n1

# Higher dimensional test

x_all = hcat(iris[:PetalLength], iris[:PetalWidth], iris[:SepalLength], iris[:SepalWidth])

subset_n = 120
x_train = x_all[1:subset_n, 1:4]
x_test = x_all[(subset_n + 1):150, 1:4]
y_in = array(iris[:Species][1:subset_n])

test_predictions_3 = kNN(x_train, y_in, 3, x_test)
class_3_big = sum(real_y .== test_predictions_3) / n1


text:
# function to more accurately compute cos(x + xn)
# see https://github.com/heltonmc/Bessels.jl/pull/13
# written by @oscardssmith
function cos_sum(x, xn)
    s = x + xn
    n, r = reduce_pi02_med(s)
    lo = r.lo - ((s - x) - xn)
    hi = r.hi + lo
    y = Base.Math.DoubleFloat64(hi, r.hi-hi+lo)
    n = n&3
    if n == 0
        return Base.Math.cos_kernel(y)
    elseif n == 1
        return -Base.Math.sin_kernel(y)
    elseif n == 2
        return -Base.Math.cos_kernel(y)
    else
        return Base.Math.sin_kernel(y)
    end
end
# function to more accurately compute sin(x + xn)
function sin_sum(x, xn)
    s = x + xn
    n, r = reduce_pi02_med(s)
    lo = r.lo - ((s - x) - xn)
    hi = r.hi + lo
    y = Base.Math.DoubleFloat64(hi, r.hi-hi+lo)
    n = n&3
    if n == 0
        return Base.Math.sin_kernel(y)
    elseif n == 1
        return Base.Math.cos_kernel(y)
    elseif n == 2
        return -Base.Math.sin_kernel(y)
    else
        return -Base.Math.cos_kernel(y)
    end
end
@inline function reduce_pi02_med(x::Float64)
    pio2_1 = 1.57079632673412561417e+00

    fn = round(x*(2/pi))
    r  = muladd(-fn, pio2_1, x)
    w  = fn * 6.07710050650619224932e-11
    y = r-w
    return unsafe_trunc(Int, fn), Base.Math.DoubleFloat64(y, (r-y)-w)
end


text:
const YPCommandsList = ["start-game",
                        "stop-game",
                        "accept-card",
                        "reject-card"]

function YPServer(msg, currws)
    !inRoom(currws) && return
    querytype = msg["querytype"]

    if querytype == "start-game"
        startnewYP(currws)
    elseif querytype == "stop-game"
        stopYP(currws)
    elseif querytype == "accept-card"
        acceptCard(currws)
    elseif querytype == "reject-card"
        rejectCard(currws)
    end
    return
end

function startnewYP(currws)
    roomname = getRoomfromWS[currws]
    handle = getHandlefromWS[currws]

    if haskey(RoomGames, roomname)
        #Ongoing Game
        return errormsg(currws, "room-game-ongoing")
    end

    playerlist = getPlayers(roomname)

    if length(playerlist) <= 2
        return errormsg(currws, "room-insufficient-players")
    elseif length(playerlist) > 5
        # Too many players
        # No error for now
    elseif !(handle in playerlist)
        return errormsg(currws, "room-not-in-game-start")
    end

    RoomGames[roomname] = g = initialiseYP(playerlist)

    msg = Dict{String, Any}()
    msg["responsetype"] = "game-game-started"
    msg["game-status"] = gameState(g)
    return broadcastmsg(roomname, JSON.json(msg))
end

# Overloaded required method
inGameUsers(g::YesPleaseGame) = g.UserList

function stopYP(currws)
    roomname = getRoomfromWS[currws]
    handle = getHandlefromWS[currws]

    # If no games are running, end
    (!haskey(RoomGames, roomname)) && return
    # Stopping Player not currently playing
    !(handle in inGameUsers(RoomGames[roomname])) && return errormsg(currws, "room-not-in-game")

    removeGameRef(roomname)
    return
end

function acceptCard(currws::WebSocket)
    roomname = getRoomfromWS(currws)
    handle = getHandlefromWS(currws)
    g = RoomGames[roomname]
    # Player not currently playing
    !(handle in inGameUsers(g)) && return errormsg(currws, "room-not-in-game")

    exitcode = acceptCard(g, handle)
    if exitcode == "game-ended"
        return errormsg(currws, "game-ended")
    elseif exitcode == "wrong-player"
        return errormsg(currws, "game-wrong-player")
    elseif exitcode == "success"
        return broadcastmsg(roomname, "{ \"responsetype\" : \"game-card-accepted\" }")
    end
end

function rejectCard(currws::WebSocket)
    roomname = getRoomfromWS(currws)
    handle = getHandlefromWS(currws)
    g = RoomGames[roomname]
    # Player not currently playing
    !(handle in inGameUsers(g)) && return errormsg(currws, "room-not-in-game")

    exitcode = rejectCard(g, handle)
    if exitcode == "game-ended"
        return errormsg(currws, "game-ended")
    elseif exitcode == "wrong-player"
        return errormsg(currws, "game-wrong-player")
    elseif exitcode == "insufficient-tokens"
        return errormsg(currws, "game-insufficient-tokens")
    elseif exitcode == "success"
        return broadcastmsg(roomname, "{ \"responsetype\" : \"game-card-rejected\" }")
    end
end


text:
function sum(dv::AbstractDistributionVector{<:Normal}, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractVector{Bool} = Falses(length(dv)))
    length(dv) == length(isgapfilled) || error(
        "argument gapfilled must have the same length as dv ($(length(dv))" *
        "but was $(length(isgapfilled)).")
    if isa(ms, HandleMissingStrategy)
        # need to allocate anyway with subsetting
        nonmissing = findall(.!ismissing.(Ref(dv), axes(dv,1)))
        if !isempty(nonmissing)
            return(sum(@inbounds(dv[nonmissing]),
                isgapfilled = @inbounds(isgapfilled[nonmissing])))
        end
    end
    # uncorrelated, only sum diagonal
    Ssum = s = Ssumnonfilled = zero(eltype(nonmissingtype(eltype(dv))))
    nterm = 0
    for (i,d) in enumerate(dv)
        μ,σ = params(d)
        Ssum += μ
        if !isgapfilled[i]
            Ssumnonfilled += μ
            s += abs2(σ)
            nterm += 1
        end
    end
    nterm > 0 || error("Expected at least one nonmissing term, but mu = $μ")
    relerr = √s/Ssumnonfilled
    Normal(Ssum, Ssum * relerr)
end

function sum(dv::AbstractDistributionVector{D},
    acf::AutoCorrelationFunction, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractVector{Bool}=Falses(length(dv)),
    storage::AbstractVector{Union{Missing,DS}} =
       Vector{Union{Missing,eltype(D)}}(undef, length(dv))) where
    {D<:Normal, DS<:eltype(D)}
    # currently use Matrix-method, Maybe implement fast with loop
    corrM = Symmetric(cormatrix_for_acf(length(dv), coef(acf)))
    return(sum_normals(
        dv, corrM, ms;
        isgapfilled=isgapfilled, storage = storage))
end

function sum(dv::AbstractDistributionVector{D},
    corr::Symmetric, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractArray{Bool,1}=Falses(length(dv)),
    storage::AbstractVector{Union{Missing,DS}} =
       Vector{Union{Missing,eltype(D)}}(undef, length(dv))) where
    {D<:Normal, DS<:eltype(D)}
    sum_normals(dv, corr, ms; isgapfilled=isgapfilled, storage=storage)
end

function sum_normals(dv::AbstractDistributionVector{D},
    corr::Symmetric, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractArray{Bool,1} = Falses(length(dv)),
    storage::AbstractVector{Union{Missing,DS}} =
        Vector{Union{Missing,eltype(D)}}(undef, length(dv))) where
    {D<:Normal, DS<:eltype(D)}
    μ = params(dv, Val(1))
    σ = params(dv, Val(2))
    # var_sum (s) is the sum across all Sigma, i.e. σT * corr * σ
    # missings and gapfilled values do not count -> set to zero
    if isa(ms, HandleMissingStrategy)
        # check on missings in dv, σ may have more finite values
        storage .= ifelse.(ismissing.(Ref(dv),axes(dv,1)), zero(DS), σ) #coalesce(σ, zero(DS))
        Spure = disallowmissing(storage)
    else
        Spure = disallowmissing(σ)
    end
    Spure[isgapfilled] .= zero(DS)
    s = transpose(Spure) * corr * Spure
    Ssum = Ssumnonfilled = zero(DS)
    nterm = 0
    for (i,d) in enumerate(dv)
        ismissing(d) && continue
        μ,σ = params(d)
        Ssum += μ
        if !isgapfilled[i]
            Ssumnonfilled += μ
            nterm += 1
        end
    end
    nterm > 0 || error("Expected at least one nonmissing term, but mu = $μ")
    relerr = √s/Ssumnonfilled
    Normal(Ssum, Ssum * relerr)
end

mean(dv::AbstractDistributionVector{<:Normal},
    ms::MissingStrategy=PassMissing(); kwargs...) =
    mean_normals(dv, ms; kwargs...)
mean(dv::AbstractDistributionVector{<:Normal}, corr::Symmetric,
    ms::MissingStrategy=PassMissing(); kwargs...) =
    mean_normals(dv, corr, ms; kwargs...)
mean(dv::AbstractDistributionVector{<:Normal}, acf::AutoCorrelationFunction,
    ms::MissingStrategy=PassMissing(); kwargs...) =
    mean_normals(dv, acf, ms; kwargs...)

function mean_normals(dv::AbstractDistributionVector{<:Normal}, x...; kwargs...)
    ds = sum(dv, x...; kwargs...)
    n = count(i -> !ismissing(dv,i), axes(dv,1))
    Normal(ds.μ/n, ds.σ/n)
end


text:
@testset "done" begin
    c = done("v", 1)

    @test typeof(c) <: Continuation
    @test typeof(c) == Continuation{Tuple{String, Int64}}
    @test c.data == ("v", 1)
    @test c.state == Done
end

@testset "next" begin
    c = next("v", 1)

    @test typeof(c) <: Continuation
    @test typeof(c) == Continuation{Tuple{String, Int64}}
    @test c.data == ("v", 1)
    @test c.state == Next
end

@testset "isdone" begin
    @test isdone(done("v", 1)) == true
    @test isdone(next("v", 1)) == false
end

@testset "isnext" begin
    @test isnext(done("v", 1)) == false
    @test isnext(next("v", 1)) == true
end


text:
reponame = "QuantEcon/SimpleDifferentialOperators.jl"
path = "docs/notebooks"

# vanilla call
activate_github_path(reponame, path = path, activate = false)
@test "Project.toml" ∈ readdir()
@test "Manifest.toml" ∈ readdir()

# vanilla call repeated (should return 0)
@test_logs (:warn, "There's already a Project.toml in the current directory, and force = false.") activate_github_path(reponame, path = path, force = false, activate = false)

# clean
rm("Project.toml")
rm("Manifest.toml")

# error handling
@test_throws Exception activate_github_path(reponame * "flerg", path, force = false, activate = false)
@test_throws Exception activate_github_path(reponame, path = "src") # no TOML here


text:
module gazel
include("feed.jl")
include("cal_contr.jl")
include("control.jl")

"""
    __init__()

init function
"""
function __init__()
    #check_deps()

    # Get the name of the display
    display = _getDisplay()

    # setup constants
    eqn(i,p,s) =  i * 0.3
    i0 = 0.5
    s_span = (0.0,10.0)

    # setup problem
    prob = ODEProblem(eqn,i0,s_span)

    # solve the ODE and convert the descretized soln into an array
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, dt=0.1, save_everystep=true)
    arsol = Array(sol)

    # start the feed
    _start(display, arsol)

end

end # Gazel


text:
module IDRsSolver

#### IterativeSolver.jl reporting protocol
export idrs, syl, stein, idrs_core

type ConvergenceHistory{T, R}
    isconverged::Bool
    threshold::T
    mvps::Int
    residuals::R
end

function empty!(ch::ConvergenceHistory)
    ch.isconverged = false
    ch.mvps = 0
    empty!(ch.residuals)
    ch
end

function push!(ch::ConvergenceHistory, resnorm::Number)
    push!(ch.residuals, resnorm)
    ch
end

####

# define abstract versions of rand, dot, norm and axpy!

arand!(C) = rand!(C)
#adot(A,B) = vecdot(A, B)
adot(A,B) = dot(A[:], B[:])
anorm(A) = vecnorm(A)
aaxpy!(a, X, Y) = BLAS.axpy!(a, X, Y)

function omega(t, s, angle = .7)
    ns = anorm(s)
    nt = anorm(t)
    ts = adot(t,s)
    rho = abs(ts/(nt*ns))
    om = ts/(nt*nt)
    if rho < angle
        om = om*angle/rho
    end
    om
end

linsys_op(x, A) = A*x
idrs(A, b, x0 = zeros(b); s = 8, tol = sqrt(eps(anorm(b))), maxiter = length(b)^2) = idrs_core(linsys_op, (A,), b, x0; s = s, tol=tol, maxiter=maxiter)

stein_op(X, A, B) = X + A*X*B
stein(A, B, C, X0 = zeros(C); s = 8, tol = sqrt(eps(anorm(C))), maxiter = length(C)^2) = idrs_core(stein_op, (A,B), C, X0; s = s, tol=tol, maxiter=maxiter)

syl_op(X, A, B) = A*X + X*B
syl(A, B, C, X0 = zeros(C); s = 8, tol = sqrt(eps(anorm(C))), maxiter = length(C)^2) = idrs_core(syl_op, (A,B), C, X0; s = s, tol=tol, maxiter=maxiter)

function idrs_core{T}(op, args, C::T, X0 = zero(C); s = 8, tol = sqrt(eps(anorm(C))), maxiter = length(C)^2)

    X = X0
    R = C - op(X, args...)
    normR = anorm(R)
    res = [normR]
	iter = 0

    if normR <= tol           # Initial guess is a good enough solution
        return X0, ConvergenceHistory(res[end] < tol, tol, length(res), res)
    end

    Z = zero(C)

    P = T[arand!(copy(C)) for k in 1:s]
    U = T[copy(Z) for k in 1:s]
    G = T[copy(Z) for k in 1:s]
    Q = copy(Z)
    V = copy(Z)

    M = eye(eltype(C),s,s)
    f = zeros(eltype(C),s)
    c = zeros(eltype(C),s)

    om = 1.
    iter = 0
    while normR > tol && iter < maxiter
        for i in 1:s,
            f[i] = adot(P[i], R)
        end
        for k in 1:s

            # Solve small system and make v orthogonal to P

            c = \(LowerTriangular(M[k:s,k:s]),f[k:s])
            copy!(V, G[k])
            scale!(c[1], V)

            copy!(Q, U[k])
            scale!(c[1], Q)
            for i = k+1:s
                aaxpy!(c[i-k+1], G[i], V)
                aaxpy!(c[i-k+1], U[i], Q)
            end

            # Compute new U[:,k] and G[:,k], G[:,k] is in space G_j

            #V = R - V
            scale!(-1., V)
            aaxpy!(1., R, V)

            copy!(U[k], Q)
            aaxpy!(om, V, U[k])
            G[k] = op(U[k], args...)

            # Bi-orthogonalise the new basis vectors

            for i in 1:k-1
                alpha = adot(P[i],G[k])/M[i,i]
                aaxpy!(-alpha, G[i], G[k])
                aaxpy!(-alpha, U[i], U[k])
            end

            # New column of M = P'*G  (first k-1 entries are zero)

            for i in k:s
                M[i,k] = adot(P[i],G[k])
            end

            #  Make r orthogonal to q_i, i = 1..k

            beta = f[k]/M[k,k]
            aaxpy!(-beta, G[k], R)
            aaxpy!(beta, U[k], X)

            normR = anorm(R)
            res = [res; normR]
            iter += 1
            if normR < tol || iter > maxiter
                return X, ConvergenceHistory(normR < tol, tol, length(res), res)
            end
            if k < s
                f[k+1:s] = f[k+1:s] - beta*M[k+1:s,k]
            end

        end

        # Now we have sufficient vectors in G_j to compute residual in G_j+1
        # Note: r is already perpendicular to P so v = r

        copy!(V, R)
        Q = op(V, args...)
        om = omega(Q, R, 0.7)
        aaxpy!(-om, Q, R)
        aaxpy!(om, V, X)

        normR = anorm(R)
        iter += 1
        res = [res; normR]
    end
    if normR > tolc
        warn("idrs() did not converge")
    end
    return X, ConvergenceHistory(res[end]<tol, tol, length(res), res)
end

end # module


text:
using Distributed
addprocs(5)

using LinearAlgebra
using GeoPhyInv
using Test
using Misfits
using BenchmarkTools
using Profile

function initialize(fp)
	println(" *********************************")
	println(" *********** $(fp) ***************")
	println(" *********************************")
	include(fp)
end

folder="Poisson"
for t in ["testscript_RandomEigenfns", "adj_state_expt", "adj_state", "testdAdx", "forw", "test_born"]
	fp = joinpath(folder, string(t, ".jl"))
	initialize(fp)
end

folder="Interpolation"
for t in ["interp_tests"]
	fp = joinpath(folder, string(t, ".jl"))
	initialize(fp)
end

folder="Models"
for t in ["Models", "param_adj"]
	fp = joinpath(folder, string(t, ".jl"))
	initialize(fp)
end

folder="Fdtd"
for t in ["accuracy", "backprop", "rho_projection", "create_snaps"]
	fp = joinpath(folder, string(t, ".jl"))
	initialize(fp)
end

folder="FWI"
for t in ["gradient_accuracy", "born_map"]
	fp = joinpath(folder, string(t, ".jl"))
	initialize(fp)
end

#tests = ["Fdtd_accuracy", "Fdtd_backprop",
#	 # "Poisson_testscript_RandomEigenfns",
#	"Models", "Data", "FWI_grad_test"]
#
#println("Running tests:")
#
#end


text:
module LaverySpline

import PPInterpolation: PP, computeDEF
using JuMP
using Cbc
#using GLPK
export makeLaverySpline, makeQuinticLaverySpline

#TODO remove the various println, move to an extra param optional

function makeLaverySpline(x::AbstractArray{TX}, y::AbstractArray{T}) where {T,TX}
    pp = PP(3, T, TX, length(y))
    computeLaverySpline(pp, x, y)
    return pp
end

function computeLaverySpline(pp::PP{3,T,TX}, x::AbstractArray{TX}, z::AbstractArray{T}) where {T,TX}
    h = x[2:end] - x[1:end-1] #h1 = x2-x1
    dz = (z[2:end] - z[1:end-1]) ./ h
    size = length(z)
    nIntervals = 10
    firstSize = (size - 1) * nIntervals #first sum from 1 to size-1 (included) and interior sum from 1 to nInterval
    #1,1 1,2, ..., 1,nIntervals ...,size-1,1, size-1,2, ..., size-1,nIntervals
    eps = 1e-4
    m = Model(with_optimizer(Cbc.Optimizer))
    # m = Model(with_optimizer(GLPK.Optimizer, tm_lim = 60.0, msg_lev = GLPK.OFF))
    @variable(m, bv[1:size])
    @variable(m, s0[1:firstSize] >= 0) #first slack
    # @objective(m, Min, sum(s0))

    @variable(m, s1[1:size-2] >= 0) #second slack
    @objective(m, Min, sum(s0) + eps * sum(s1)) #s0 and s1 don't have the same size! => add s0 to s1 elements! or expand s0 and s1
    @constraint(m, s1 .>= bv[2:size-1])
    @constraint(m, s1 .>= -bv[2:size-1])
    deltatk = 1.0 / (nIntervals)

    for k = 1:nIntervals
        for i = 1:size-1
            dzi = dz[i]
            tk = -0.5 + (k - 0.5) / nIntervals
            @constraint(m, s0[k+(i-1)*nIntervals] >= ((-1 + 6 * tk) * bv[i] + (1 + 6 * tk) * bv[i+1] - 12 * dzi * tk) * deltatk)
            @constraint(m, s0[k+(i-1)*nIntervals] >= -((-1 + 6 * tk) * bv[i] + (1 + 6 * tk) * bv[i+1] - 12 * dzi * tk) * deltatk)
        end
    end
    @constraint(m, bv[1] == dz[1] - 0.5 * (bv[2] - dz[1]))
    @constraint(m, bv[size] == dz[size-1] - 0.5 * (bv[size-1] - dz[size-1]))
    optimize!(m)

    println("Objective: ", objective_value(m))

    b = zeros(size)

    for i = 1:size
        println("b $i: ", value(bv[i]))
        b[i] = value(bv[i])
    end
    #b[1] = dz[1]-0.5*(b[2]-dz[1]) #natural BC
    #b[end] = dz[end]-0.5*(b[end-1]-dz[end]) #natural BC

    c = pp.c
    for i = 1:n-1
        c[1, i] = (3 * dz[i] - b[i+1] - 2 * b[i]) / h[i]
        c[2, i] = (b[i+1] + b[i] - 2 * dz[i]) / (h[i]^2)
    end
    pp.a[1:end] = z
    pp.b[1:end] = b
    pp.x[1:end] = x
end

function makeQuinticLaverySpline(x::AbstractArray{TX}, z::AbstractArray{T}, eps1 = 1e-3, eps2 = 1e-3, order = 3) where {T,TX}
    h = x[2:end] - x[1:end-1] #h1 = x2-x1
    dz = (z[2:end] - z[1:end-1]) ./ h
    size = length(z)
    nIntervals = 10
    firstSize = (size - 1) * nIntervals #first sum from 1 to size-1 (included) and interior sum from 1 to nInterval
    #1,1 1,2, ..., 1,nIntervals ...,size-1,1, size-1,2, ..., size-1,nIntervals

    #first: compute slopes for a cubic.
    m = Model(with_optimizer(Cbc.Optimizer))
    # m = Model(with_optimizer(GLPK.Optimizer, tm_lim = 60.0, msg_lev = GLPK.OFF))
    @variable(m, bv[1:size])
    @variable(m, s0[1:firstSize] >= 0) #first slack
    # @objective(m, Min, sum(s0))

    @variable(m, s1[1:size-2] >= 0) #second slack
    @objective(m, Min, sum(s0) + eps1 * sum(s1)) #s0 and s1 don't have the same size! => add s0 to s1 elements! or expand s0 and s1
    @constraint(m, s1 .>= bv[2:size-1])
    @constraint(m, s1 .>= -bv[2:size-1])
    deltatk = 1.0 / (nIntervals)

    for k = 1:nIntervals
        for i = 1:size-1
            dzi = dz[i]
            tk = -0.5 + (k - 0.5) / nIntervals
            @constraint(m, s0[k+(i-1)*nIntervals] >= ((-1 + 6 * tk) * bv[i] + (1 + 6 * tk) * bv[i+1] - 12 * dzi * tk) * deltatk)
            @constraint(m, s0[k+(i-1)*nIntervals] >= -((-1 + 6 * tk) * bv[i] + (1 + 6 * tk) * bv[i+1] - 12 * dzi * tk) * deltatk)
        end
    end
    @constraint(m, bv[1] == dz[1] - 0.5 * (bv[2] - dz[1]))
    @constraint(m, bv[size] == dz[size-1] - 0.5 * (bv[size-1] - dz[size-1]))
    optimize!(m)

    println("Objective: ", objective_value(m))

    b = zeros(T, size)

    for i = 1:size
        # println("b $i: ", value(bv[i]))
        b[i] = value(bv[i])
    end

    #second: compute second derivative for quintic

    m = Model(with_optimizer(Cbc.Optimizer))
    # m = Model(with_optimizer(GLPK.Optimizer, tm_lim = 60.0, msg_lev = GLPK.OFF))
    @variable(m, cv[1:size])
    @variable(m, s0[1:firstSize] >= 0) #first slack
    @variable(m, s1[1:size-2] >= 0) #second slack
    @objective(m, Min, sum(s0) + eps2 * sum(s1)) #s0 and s1 don't have the same size! => add s0 to s1 elements! or expand s0 and s1
    @constraint(m, s1 .>= cv[2:size-1])
    @constraint(m, s1 .>= -cv[2:size-1])

    if order == 1
        for k = 1:nIntervals
            for i = 1:size-1
                deltatk = (x[i+1] - x[i]) / (nIntervals)
                hk = deltatk * (k - 0.5)
                @constraint(m, s0[k+(i-1)*nIntervals] >= (b[i] + hk * cv[i] + 3 * hk^2 * ((cv[i+1] - 3 * cv[i]) / (2 * h[i]) + 2 * (5 * dz[i] - 3 * b[i] - 2 * b[i+1]) / (h[i]^2)) + 4 * hk^3 * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 5 * (hk^4) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
                @constraint(m, -s0[k+(i-1)*nIntervals] <= (b[i] + hk * cv[i] + 3 * hk^2 * ((cv[i+1] - 3 * cv[i]) / (2 * h[i]) + 2 * (5 * dz[i] - 3 * b[i] - 2 * b[i+1]) / (h[i]^2)) + 4 * hk^3 * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 5 * (hk^4) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
            end
        end
    elseif order == 2
        #second der constraint
        for k = 1:nIntervals
            for i = 1:size-1
                deltatk = (x[i+1] - x[i]) / (nIntervals)
                hk = deltatk * (k - 0.5)
                @constraint(m, s0[k+(i-1)*nIntervals] >= (cv[i] + 6 * hk * ((cv[i+1] - 3 * cv[i]) / (2 * h[i]) + 2 * (5 * dz[i] - 3 * b[i] - 2 * b[i+1]) / (h[i]^2)) + 12 * hk^2 * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 20 * (hk^3) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
                @constraint(m, -s0[k+(i-1)*nIntervals] <= (cv[i] + 6 * hk * ((cv[i+1] - 3 * cv[i]) / (2 * h[i]) + 2 * (5 * dz[i] - 3 * b[i] - 2 * b[i+1]) / (h[i]^2)) + 12 * hk^2 * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 20 * (hk^3) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
            end
        end
        #third der
    elseif order == 3
        for k = 1:nIntervals
            for i = 1:size-1
                deltatk = (x[i+1] - x[i]) / (nIntervals)
                hk = deltatk * (k - 0.5)
                @constraint(m, s0[k+(i-1)*nIntervals] >= (6 * ((cv[i+1] - 3 * cv[i]) / (2 * h[i]) + 2 * (5 * dz[i] - 3 * b[i] - 2 * b[i+1]) / (h[i]^2)) + 24 * hk * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 60 * (hk^2) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
                @constraint(m, -s0[k+(i-1)*nIntervals] <= (6 * ((cv[i+1] - 3 * cv[i]) / (2 * h[i]) + 2 * (5 * dz[i] - 3 * b[i] - 2 * b[i+1]) / (h[i]^2)) + 24 * hk * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 60 * (hk^2) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
            end
        end
    elseif order == 4
        for k = 1:nIntervals
            for i = 1:size-1
                deltatk = (x[i+1] - x[i]) / (nIntervals)
                hk = deltatk * (k - 0.5)
                @constraint(m, s0[k+(i-1)*nIntervals] >= (24 * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 120 * (hk) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
                @constraint(m, -s0[k+(i-1)*nIntervals] <= (24 * ((3 * cv[i] - 2 * cv[i+1]) / (2 * h[i]^2) + (8 * b[i] + 7 * b[i+1] - 15 * dz[i]) / (h[i]^3)) + 120 * (hk) * ((cv[i+1] - cv[i]) / (2 * h[i]^3) + 3 * (2 * dz[i] - b[i+1] - b[i]) / (h[i]^4))) * deltatk)
            end
        end
    end
    @constraint(m, cv[1] == 0)
    @constraint(m, cv[size] == 0)
    optimize!(m)

    println("Objective: ", objective_value(m))

    y2 = zeros(T, size)

    for i = 1:size
        #        println("c $i: ", value(cv[i]))
        y2[i] = value(cv[i])
    end

    #b[1] = dz[1]-0.5*(b[2]-dz[1]) #natural BC
    #b[end] = dz[end]-0.5*(b[end-1]-dz[end]) #natural BC

    c = zeros(T, (4, size))
    computeDEF(c, b, y2, dz, h)
    return PP(5, copy(z), b, c, copy(x))
end
end


text:
"""
    $(SIGNATURES)

Calculate ΔG'⁰ for each `rxn_string` in `rxn_strings`. The `uncertainty_representation` can
be either `fullrank`, `cov`, or `sqrt`, and is used to represent the uncertainties. Specifically,
`cov` returns a full covariance matrix, `precision` returna the precision matrix (i.e.
the inverse of the covariance matrix), `sqrt` returns the square root of the covariance,
based on the uncertainty vectors, and `fullrank` returns a full-rank square root of the
covariance that is a compressed form of the sqrt result. If `minimize_norm` is  `true`,
an orthogonal projection is used to minimize the Euclidian norm of the result vector. If
`balance_warn` is true, emit a warning if any reaction is unbalanced, if `skip_unbalanced`
is true, return `nothing` if an unbalanced reaction is found.

# Example
```
equilibrator = Equilibrator()

ntp3 = "gtp + h2o = gdp + pi"
adk3 = "adp + gdp = amp + gtp"

μ, Q = standard_dg_prime_multi(eq, bigg.([ntp3, adk3]); uncertainty_representation="fullrank")
```
"""
function standard_dg_prime_multi(
    equilibrator::Equilibrator,
    rxn_strings::Vector{String};
    balance_warn = true,
    skip_unbalanced = false,
    uncertainty_representation = "cov",
    minimize_norm = false,
)
    rxns =
        [equilibrator.cc.parse_reaction_formula(rxn_string) for rxn_string in rxn_strings]

    if balance_warn && any(!rxn.is_balanced() for rxn in rxns)
        for (k, rxn) in enumerate(rxns)
            if !rxn.is_balanced
                @warn "Reaction at index $k is unbalanced."
            end
        end
        skip_unbalanced && return nothing
    end

    pymeans, pyQ = equilibrator.cc.standard_dg_prime_multi(
        rxns,
        uncertainty_representation = uncertainty_representation,
        minimize_norm = minimize_norm,
    )

    if uncertainty_representation == "cov"
        Q = pyQ.m_as("kJ^2/mol^2")u"kJ^2/mol^2"
    elseif uncertainty_representation == "fullrank"
        Q = pyQ.m_as("kJ/mol")u"kJ/mol"
    elseif uncertainty_representation == "precision"
        Q = pyQ.m_as("mol/kJ")u"mol/kJ"
    elseif uncertainty_representation == "sqrt"
        Q = pyQ.m_as("kJ/mol")u"kJ/mol"
    end
    return pymeans.m_as("kJ/mol")u"kJ/mol", Q
end


text:
######################################################################################
# Prints a single JSON file benchmark, or the comparison of two benchmarks, that were
# previously created by run_benchmark.jl
######################################################################################
include("benchmark_utils.jl")

if !(1 <= length(ARGS) <= 2)
    println("Specify one or two benchmark JSON files as input:\n\n$(basename(@__FILE__)) baseline_file [alternative_file]")
else
    length(ARGS) > 1 ? print_benchmark_comparison(ARGS[1], ARGS[2]) : print_benchmark(ARGS[1])
end

max_issues_repo_name: maxbennedich/julia-regression-analysis

max_forks_repo_name: maxbennedich/julia-regression-analysis


text:
using HypothesisTests, Base.Test

t = BinomialTest(26, 78)
@test pvalue(t) ≈ 0.004334880883507431
@test pvalue(t, tail=:left) ≈ 0.002167440441753716
@test pvalue(t, tail=:right) ≈ 0.9989844298129187
@test_ci_approx ci(t) (0.23058523962930383, 0.4491666887959782)
@test_ci_approx ci(t, tail=:left) (0.0, 0.4313047758370174)
@test_ci_approx ci(t, tail=:right) (0.2451709633730693, 1.0)
@test_ci_approx ci(t, method=:wald) (0.22871819521037956, 0.43794847145628707)
@test_ci_approx ci(t, tail=:left, method=:wald) (0.0, 0.42112912485444692)
@test_ci_approx ci(t, tail=:right, method=:wald) (0.24553754181221971, 1.0)
@test_ci_approx ci(t, method=:wilson) (0.23872670036358601, 0.44358590287381217)
@test_ci_approx ci(t, tail=:left, method=:wilson) (0.0, 0.42541288951088108)
@test_ci_approx ci(t, tail=:right, method=:wilson) (0.25242832328277831, 1.0)
@test_ci_approx ci(t, method=:jeffrey) (0.23626570247518358, 0.44251318323879296)
@test_ci_approx ci(t, tail=:left, method=:jeffrey) (0.0, 0.42466492683653623)
@test_ci_approx ci(t, tail=:right, method=:jeffrey) (0.25098836986261724, 1.0)
@test_ci_approx ci(t, method=:agresti_coull) (0.2384423809121706, 0.44387022232522744)
@test_ci_approx ci(t, tail=:left, method=:agresti_coull) (0.0, 0.42558712894362222)
@test_ci_approx ci(t, tail=:right, method=:agresti_coull) (0.25225408385003706, 1.0)
show(IOBuffer(), t)

t = BinomialTest([trues(6); falses(3)])
@test pvalue(t) ≈ 0.5078125000000002
@test pvalue(t, tail=:left) ≈ 0.91015625
@test pvalue(t, tail=:right) ≈ 0.2539062500000001
@test_ci_approx ci(t) (0.2992950562085405, 0.9251453685803082)
@test_ci_approx ci(t, tail=:left) (0.0, 0.9022531865607242)
@test_ci_approx ci(t, tail=:right) (0.3449413659437032, 1.0)
@test_ci_approx ci(t, method=:wald) (0.35868803903340479, 0.97464529429992841)
@test_ci_approx ci(t, tail=:left, method=:wald) (0.0,0.92513047859481645)
@test_ci_approx ci(t, tail=:right, method=:wald) (0.40820285473851681,1.0)
@test_ci_approx ci(t, method=:wilson) (0.35420213558039609,0.87941618161308899)
@test_ci_approx ci(t, tail=:left, method=:wilson) (0.0,0.85802909820500495)
@test_ci_approx ci(t, tail=:right, method=:wilson) (0.39825972868840931,1.0)
@test_ci_approx ci(t, method=:jeffrey) (0.34779179347226591,0.89578677833922582)
@test_ci_approx ci(t, tail=:left, method=:jeffrey) (0.0,0.86830830610561005)
@test_ci_approx ci(t, tail=:right, method=:jeffrey) (0.39604343455469687,1.0)
@test_ci_approx ci(t, method=:agresti_coull) (0.350905767251112,0.88271254994237336)
@test_ci_approx ci(t, tail=:left, method=:agresti_coull) (0.0,0.86049746046629294)
@test_ci_approx ci(t, tail=:right, method=:agresti_coull) (0.39579136642712159,1.0)
show(IOBuffer(), t)

x = [55, 58, 61, 61, 62, 62, 62, 63, 63, 64, 66, 68, 68, 69, 69, 69, 70, 71, 72, 72]
@test pvalue(SignTest(x)) ≈ 1.907348632812499e-6
@test pvalue(SignTest(x), tail=:left) ≈ 1.0
@test pvalue(SignTest(x), tail=:right) ≈ 9.536743164062495e-7
@test pvalue(SignTest(x, 70)) ≈ 0.004425048828125003
@test pvalue(SignTest(x, 70), tail=:left) ≈ 0.0022125244140625013
@test pvalue(SignTest(x, 70), tail=:right) ≈ 0.9996356964111328
@test_ci_approx ci(SignTest(x, 70)) (62, 69)
@test_ci_approx ci(SignTest(x, 70), 0.0002) (61, 71)
show(IOBuffer(), SignTest(x, 70))

x = [9, 2, 7, 5]
y = [7, 2, 6, 4]
@test pvalue(SignTest(x, y)) ≈ 0.25
@test_ci_approx ci(SignTest(x, y)) (0, 2)
show(IOBuffer(), SignTest(x, y))

# www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf
x = [-4.7, 3.7, 22.4, 13.6, 8.7, 9.1, -7.8, 10.8, 15.6, 23.5, 14.4, 20.2, 6.5, 10.1, -6.9]
@test pvalue(SignTest(x)) ≈ 0.03515625
@test pvalue(SignTest(x), tail=:left) ≈ 0.996307373046875
@test pvalue(SignTest(x), tail=:right) ≈ 0.017578125000000007
@test_ci_approx ci(SignTest(x), 1-0.99995) (-7.8, 23.5)
@test_ci_approx ci(SignTest(x), 1-0.9999) (-6.9, 22.4)
@test_ci_approx ci(SignTest(x), 1-0.993) (-4.7, 20.2)
@test_ci_approx ci(SignTest(x), 1-0.965) (3.7, 15.6)
@test_ci_approx ci(SignTest(x), 1-0.882) (6.5, 14.4)
@test_ci_approx ci(SignTest(x), 1-0.7) (8.7, 13.6)
@test_ci_approx ci(SignTest(x), 1-0.6) (9.1, 10.8)
show(IOBuffer(), SignTest(x))


text:
# -*- coding: utf-8 -*-
"""
    halfar_solution(t, r, θ)

Returns the evaluation of the Halfar solutions for the SIA equation.

Arguments:
    - t: time
    - r: radial distance. The solutions have polar symmetry around the center of origin
    - ν = (A, H₀, R₀)
"""
function halfar_solution(t, r, ν)

    # parameters of Halfar solutions
    A, h₀, r₀ = ν

    Γ = 2 * A * (ρ * g)^n / (n+2)
    τ₀ = (7/4)^3 * r₀^4 / ( 18 * Γ * h₀^7 )   # characteristic time

    if r₀ * (t/τ₀)^(1/18) <= r
        return 0.0
    else
        return h₀ * (τ₀/t)^(1/9) * ( 1 - ( (τ₀/t)^(1/18) * (r/r₀) )^(4/3) )^(3/7)
    end
end


text:
export Line

## Standard interval

# angle is π*a where a is (false==0) and (true==1)
# or ranges from (-1,1].  We use 1 as 1==true.

"""
    Line{a}(c)

represents the line at angle `a` in the complex plane, centred at `c`.
"""
struct Line{angle,T<:Number} <: SegmentDomain{T}
    center::T
    α::Float64
    β::Float64

    #TODO get this inner constructor working again.
    Line{angle,T}(c,α,β) where {angle,T} = new{angle,T}(c,α,β)
    Line{angle,T}(c) where {angle,T} = new{angle,T}(c,-1.,-1.)
    Line{angle,T}() where {angle,T} = new{angle,T}(zero(T),-1.,-1.)
end

const RealLine{T} = Union{Line{false,T},Line{true,T}}

Line{a}(c,α,β) where {a} = Line{a,typeof(c)}(c,α,β)
Line{a}(c::Number) where {a} = Line{a,typeof(c)}(c)
Line{a}() where {a} = Line{a,Float64}()

angle(d::Line{a}) where {a} = a*π

reverseorientation(d::Line{true}) = Line{false}(d.center,d.β,d.α)
reverseorientation(d::Line{false}) = Line{true}(d.center,d.β,d.α)
reverseorientation(d::Line{a}) where {a} = Line{a-1}(d.center,d.β,d.α)

# ensure the angle is always in (-1,1]
Line(c,a,α,β) = Line{mod(a/π-1,-2)+1,typeof(c)}(c,α,β)
Line(c,a) = Line(c,a,-1.,-1.)
# true is negative orientation, false is positive orientation
# this is because false==0 and we take angle=0
Line(b::Bool) = Line{b}()
Line() = Line(false)

isambiguous(d::Line)=isnan(d.center)
convert(::Type{Domain{T}},d::Line{a}) where {a,T<:Number} = Line{a,T}(d.center, d.α, d.β)
convert(::Type{Line{a,T}},::AnyDomain) where {a,T<:Number} = Line{a,T}(NaN)
convert(::Type{IT},::AnyDomain) where {IT<:Line}=Line(NaN,NaN)

## Map interval

##TODO non-1 alpha,beta

isempty(::Line) = false

function line_tocanonical(α,β,x)
    @assert α==β==-1. || α==β==-.5

    if α==β==-1.
        2x/(1+sqrt(1+4x^2))
    elseif α==β==-.5
        x/sqrt(1 + x^2)
    end
end

function line_tocanonicalD(α,β,x)
    @assert α==β==-1. || α==β==-.5

    if α==β==-1.
        2/(1+4x^2+sqrt(1+4x^2))
    elseif α==β==-0.5
        (1 + x^2)^(-3/2)
    end
end
function line_fromcanonical(α,β,x)
    #TODO: why is this consistent?
    if α==β==-1.
        x/(1-x^2)
    else
        x*(1 + x)^α*(1 - x)^β
    end
end
function line_fromcanonicalD(α,β,x)
    if α==β==-1.
        (1+x^2)/(1-x^2)^2
    else
        (1 - (β-α)x - (β+α+1)x^2)*(1+x)^(α-1)*(1-x)^(β-1)
    end
end

function line_invfromcanonicalD(α,β,x)
    if α==β==-1.
        (1-x^2)^2/(1+x^2)
    else
        1/(1 - (β-α)x - (β+α+1)x^2)*(1+x)^(1-α)*(1-x)^(1-β)
    end
end

tocanonical(d::Line,x) = line_tocanonical(d.α,d.β,cis(-angle(d)).*(x-d.center))
tocanonical(d::Line{false},x) = line_tocanonical(d.α,d.β,x-d.center)
tocanonical(d::Line{true},x) = line_tocanonical(d.α,d.β,d.center-x)

tocanonicalD(d::Line,x) = cis(-angle(d)).*line_tocanonicalD(d.α,d.β,cis(-angle(d)).*(x-d.center))
tocanonicalD(d::Line{false},x) = line_tocanonicalD(d.α,d.β,x-d.center)
tocanonicalD(d::Line{true},x) = -line_tocanonicalD(d.α,d.β,d.center-x)

fromcanonical(d::Line,x) = cis(angle(d))*line_fromcanonical(d.α,d.β,x)+d.center
fromcanonical(d::Line{false},x) = line_fromcanonical(d.α,d.β,x)+d.center
fromcanonical(d::Line{true},x) = -line_fromcanonical(d.α,d.β,x)+d.center

fromcanonicalD(d::Line,x) = cis(angle(d))*line_fromcanonicalD(d.α,d.β,x)
fromcanonicalD(d::Line{false},x) = line_fromcanonicalD(d.α,d.β,x)
fromcanonicalD(d::Line{true},x) = -line_fromcanonicalD(d.α,d.β,x)

invfromcanonicalD(d::Line,x) = cis(-angle(d))*line_invfromcanonicalD(d.α,d.β,x)
invfromcanonicalD(d::Line{false},x) = line_invfromcanonicalD(d.α,d.β,x)
invfromcanonicalD(d::Line{true},x) = -line_invfromcanonicalD(d.α,d.β,x)

==(d::Line{a},m::Line{a}) where {a} = d.center == m.center && d.β == m.β &&d.α == m.α

# algebra
*(c::Real,d::Line{false}) = Line{sign(c)>0 ? false : true}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
*(c::Real,d::Line{true}) = Line{sign(c)>0 ? true : false}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
*(c::Number,d::Line) = Line(isapprox(d.center,0) ? d.center : c*d.center,angle(d)+angle(c),d.α,d.β)
*(d::Line,c::Number) = c*d
for OP in (:+,:-)
    @eval begin
        $OP(c::Number,d::Line{a}) where {a} = Line{a}($OP(c,d.center),d.α,d.β)
        $OP(d::Line{a},c::Number) where {a} = Line{a}($OP(d.center,c),d.α,d.β)
    end
end

# algebra
arclength(d::Line) = Inf
leftendpoint(d::Line) = -Inf
rightendpoint(d::Line) = Inf
complexlength(d::Line) =Inf

## vectorized

function convert(::Type{Line},d::ClosedInterval)
    a,b=d.left,d.right
    @assert abs(a) == abs(b) == Inf

    if isa(a,Real) && isa(b,Real)
        if a==Inf
            @assert b==-Inf
            Line(true)
        else
            @assert a==-Inf&&b==Inf
            Line(false)
        end
    elseif abs(real(a)) < Inf
        @assert real(a)==real(b)
        @assert sign(imag(a))==-sign(imag(b))

        Line(real(b),angle(b))
    elseif isnan(real(a)) && isnan(real(b))  # hack for -im*Inf
        Line([imag(a)*im,imag(b)*im])
    elseif isnan(real(a))  # hack for -im*Inf
        Line([real(b)+imag(a)*im,b])
    elseif isnan(real(b))  # hack for -im*Inf
        Line([a,real(a)+imag(b)*im])
    elseif abs(imag(a)) < Inf
        @assert imag(a)==imag(b)
        @assert sign(real(a))==-sign(real(b))

        Line(imag(b),angle(b))
    else
        @assert angle(b) == -angle(a)

        Line(0.,angle(b))
    end
end

issubset(a::IntervalOrSegment, b::Line) = leftendpoint(a)∈b && rightendpoint(a)∈b
issubset(a::Ray{angle}, b::Line{angle}) where angle = leftendpoint(a) ∈ b
issubset(a::Ray{true}, b::Line{false}) = true
issubset(a::Ray{false}, b::Line{true}) = true

function intersect(a::Union{Interval,Segment,Ray},b::Line)
    @assert a ⊆ b
    a
end

function union(a::Union{Interval,Segment,Ray},b::Line)
    @assert a ⊆ b
    b
end

intersect(b::Line,a::Union{Interval,Segment,Ray}) = intersect(a,b)
union(b::Line,a::Union{Interval,Segment,Ray}) = union(a,b)

function setdiff(b::Line,a::Segment)
    @assert a ⊆ b
    if leftendpoint(a)>rightendpoint(a)
        b\reverseorientation(a)
    else
        Ray([leftendpoint(b),leftendpoint(a)]) ∪ Ray([rightendpoint(a),rightendpoint(b)])
    end
end