text:
fft!(::AbstractVector{T}, ::AbstractVector{T}, ::Int, ::Int, ::Direction, ::AbstractFFTType, ::CallGraph{T}, ::Int) where {T} = nothing

@inline function direction_sign(d::Direction)
    Int(d)
end

function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{U}, start_out::Int, start_in::Int, v::Direction, t::AbstractFFTType, idx::Int) where {T,U}
    fft!(out, in, start_out, start_in, v, t, g, idx)
end

"""
$(TYPEDSIGNATURES)
Cooley-Tukey composite FFT, with a pre-computed call graph

# Arguments
`out`: Output vector
`in`: Input vector
`start_out`: Index of the first element of the output vector
`start_in`: Index of the first element of the input vector
`d`: Direction of the transform
`g`: Call graph for this transform
`idx`: Index of the current transform in the call graph

"""
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, start_out::Int, start_in::Int, d::Direction, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T,U}
    root = g[idx]
    left_idx = idx + root.left
    right_idx = idx + root.right
    left = g[left_idx]
    right = g[right_idx]
    N  = root.sz
    N1 = left.sz
    N2 = right.sz
    s_in = root.s_in
    s_out = root.s_out

    w1 = convert(T, cispi(direction_sign(d)*2/N))
    wj1 = one(T)
    tmp = g.workspace[idx]
    @inbounds for j1 in 0:N1-1
        wk2 = wj1
        g(tmp, in, N2*j1+1, start_in + j1*s_in, d, right.type, right_idx)
        j1 > 0 && @inbounds for k2 in 1:N2-1
            tmp[N2*j1 + k2 + 1] *= wk2
            wk2 *= wj1
        end
        wj1 *= w1
    end

    @inbounds for k2 in 0:N2-1
        g(out, tmp, start_out + k2*s_out, k2+1, d, left.type, left_idx)
    end
end

"""
$(TYPEDSIGNATURES)
Discrete Fourier Transform, O(N^2) algorithm, in place.

# Arguments
`out`: Output vector
`in`: Input vector
`N`: Size of the transform
`start_out`: Index of the first element of the output vector
`stride_out`: Stride of the output vector
`start_in`: Index of the first element of the input vector
`stride_in`: Stride of the input vector
`d`: Direction of the transform

"""
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T}
    tmp = in[start_in]
    @inbounds for j in 1:N-1
        tmp += in[start_in + j*stride_in]
    end
    out[start_out] = tmp

    wk = wkn = w = convert(T, cispi(direction_sign(d)*2/N))
    @inbounds for d in 1:N-1
        tmp = in[start_in]
        @inbounds for k in 1:N-1
            tmp += wkn*in[start_in + k*stride_in]
            wkn *= wk
        end
        out[start_out + d*stride_out] = tmp
        wk *= w
        wkn = wk
    end
end

function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T<:Real}
    halfN = N÷2
    wk = wkn = w = convert(Complex{T}, cispi(direction_sign(d)*2/N))

    tmpBegin = tmpHalf = in[start_in]
    @inbounds for j in 1:N-1
        tmpBegin += in[start_in + stride_in*j]
        iseven(j) ? tmpHalf += in[start_in + stride_in*j] : tmpHalf -= in[start_in + stride_in*j]
    end
    out[start_out] = convert(Complex{T}, tmpBegin)
    iseven(N) && (out[start_out + stride_out*halfN] = convert(Complex{T}, tmpHalf))

    @inbounds for d in 1:halfN
        tmp = in[start_in]
        @inbounds for k in 1:N-1
            tmp += wkn*in[start_in + k*stride_in]
            wkn *= wk
        end
        out[start_out + d*stride_out] = tmp
        wk *= w
        wkn = wk
    end
    @inbounds for k in halfN+1:N-1
        out[start_out + stride_out*k] = conj(out[start_out + stride_out*(N-k)])
    end
end

function fft!(out::AbstractVector{T}, in::AbstractVector{U}, start_out::Int, start_in::Int, d::Direction, ::DFT, g::CallGraph{T}, idx::Int) where {T,U}
    root = g[idx]
    fft_dft!(out, in, root.sz, start_out, root.s_out, start_in, root.s_in, d)
end

"""
$(TYPEDSIGNATURES)
Power of 2 FFT, in place

# Arguments
`out`: Output vector
`in`: Input vector
`N`: Size of the transform
`start_out`: Index of the first element of the output vector
`stride_out`: Stride of the output vector
`start_in`: Index of the first element of the input vector
`stride_in`: Stride of the input vector
`d`: Direction of the transform

"""
function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T, U}
    if N == 2
        out[start_out]              = in[start_in] + in[start_in + stride_in]
        out[start_out + stride_out] = in[start_in] - in[start_in + stride_in]
        return
    end
    m = N ÷ 2

    fft_pow2!(out, in, m, start_out               , stride_out, start_in            , stride_in*2, d)
    fft_pow2!(out, in, m, start_out + m*stride_out, stride_out, start_in + stride_in, stride_in*2, d)

    w1 = convert(T, cispi(direction_sign(d)*2/N))
    wj = one(T)
    @inbounds for j in 0:m-1
        j1_out = start_out + j*stride_out
        j2_out = start_out + (j+m)*stride_out
        out_j    = out[j1_out]
        out[j1_out] = out_j + wj*out[j2_out]
        out[j2_out] = out_j - wj*out[j2_out]
        wj *= w1
    end
end

function fft!(out::AbstractVector{T}, in::AbstractVector{U}, start_out::Int, start_in::Int, d::Direction, ::Pow2FFT, g::CallGraph{T}, idx::Int) where {T,U}
    root = g[idx]
    N = root.sz
    s_in = root.s_in
    s_out = root.s_out
    fft_pow2!(out, in, N, start_out, s_out, start_in, s_in, d)
end

"""
$(TYPEDSIGNATURES)
Power of 4 FFT, in place

# Arguments
`out`: Output vector
`in`: Input vector
`N`: Size of the transform
`start_out`: Index of the first element of the output vector
`stride_out`: Stride of the output vector
`start_in`: Index of the first element of the input vector
`stride_in`: Stride of the input vector
`d`: Direction of the transform

"""
function fft_pow4!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T, U}
    ds = direction_sign(d)
    plusi = ds*1im
    minusi = ds*-1im
    if N == 4
        out[start_out + 0]            = in[start_in] + in[start_in + stride_in]        + in[start_in + 2*stride_in] + in[start_in + 3*stride_in]
        out[start_out +   stride_out] = in[start_in] + in[start_in + stride_in]*plusi  - in[start_in + 2*stride_in] + in[start_in + 3*stride_in]*minusi
        out[start_out + 2*stride_out] = in[start_in] - in[start_in + stride_in]        + in[start_in + 2*stride_in] - in[start_in + 3*stride_in]
        out[start_out + 3*stride_out] = in[start_in] + in[start_in + stride_in]*minusi - in[start_in + 2*stride_in] + in[start_in + 3*stride_in]*plusi
        return
    end
    m = N ÷ 4

    @muladd fft_pow4!(out, in, m, start_out                 , stride_out, start_in              , stride_in*4, d)
    @muladd fft_pow4!(out, in, m, start_out +   m*stride_out, stride_out, start_in +   stride_in, stride_in*4, d)
    @muladd fft_pow4!(out, in, m, start_out + 2*m*stride_out, stride_out, start_in + 2*stride_in, stride_in*4, d)
    @muladd fft_pow4!(out, in, m, start_out + 3*m*stride_out, stride_out, start_in + 3*stride_in, stride_in*4, d)

    w1 = convert(T, cispi(direction_sign(d)*2/N))
    wj = one(T)

    w1 = convert(T, cispi(ds*2/N))
    w2 = convert(T, cispi(ds*4/N))
    w3 = convert(T, cispi(ds*6/N))
    wk1 = wk2 = wk3 = one(T)

    @inbounds for k in 0:m-1
        @muladd k0 = start_out + k*stride_out
        @muladd k1 = start_out + (k+m)*stride_out
        @muladd k2 = start_out + (k+2*m)*stride_out
        @muladd k3 = start_out + (k+3*m)*stride_out
        y_k0, y_k1, y_k2, y_k3 = out[k0], out[k1], out[k2], out[k3]
        @muladd out[k0] = (y_k0 + y_k2*wk2) + (y_k1*wk1 + y_k3*wk2)
        @muladd out[k1] = (y_k0 - y_k2*wk2) + (y_k1*wk1 - y_k3*wk3) * plusi
        @muladd out[k2] = (y_k0 + y_k2*wk2) - (y_k1*wk1 + y_k3*wk3)
        @muladd out[k3] = (y_k0 - y_k2*wk2) + (y_k1*wk1 - y_k3*wk3) * minusi
        wk1 *= w1
        wk2 *= w2
        wk3 *= w3
    end
end

function fft!(out::AbstractVector{T}, in::AbstractVector{U}, start_out::Int, start_in::Int, d::Direction, ::Pow4FFT, g::CallGraph{T}, idx::Int) where {T,U}
    root = g[idx]
    N = root.sz
    s_in = root.s_in
    s_out = root.s_out
    fft_pow4!(out, in, N, start_out, s_out, start_in, s_in, d)
end

"""
$(TYPEDSIGNATURES)
Power of 3 FFT, in place

# Arguments
out: Output vector
in: Input vector
N: Size of the transform
start_out: Index of the first element of the output vector
stride_out: Stride of the output vector
start_in: Index of the first element of the input vector
stride_in: Stride of the input vector
d: Direction of the transform
plus120: Depending on direction, perform either ±120° rotation
minus120: Depending on direction, perform either ∓120° rotation

"""
function fft_pow3!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction, plus120::T, minus120::T) where {T, U}
    if N == 3
        @muladd out[start_out + 0]            = in[start_in] + in[start_in + stride_in]          + in[start_in + 2*stride_in]
        @muladd out[start_out +   stride_out] = in[start_in] + in[start_in + stride_in]*plus120  + in[start_in + 2*stride_in]*minus120
        @muladd out[start_out + 2*stride_out] = in[start_in] + in[start_in + stride_in]*minus120 + in[start_in + 2*stride_in]*plus120
        return
    end

    # Size of subproblem
    Nprime = N ÷ 3

    ds = direction_sign(d)

    # Dividing into subproblems
    fft_pow3!(out, in, Nprime, start_out, stride_out, start_in, stride_in*3, d, plus120, minus120)
    fft_pow3!(out, in, Nprime, start_out + Nprime*stride_out, stride_out, start_in + stride_in, stride_in*3, d, plus120, minus120)
    fft_pow3!(out, in, Nprime, start_out + 2*Nprime*stride_out, stride_out, start_in + 2*stride_in, stride_in*3, d, plus120, minus120)

    w1 = convert(T, cispi(ds*2/N))
    w2 = convert(T, cispi(ds*4/N))
    wk1 = wk2 = one(T)
    for k in 0:Nprime-1
        @muladd k0 = start_out + k*stride_out
        @muladd k1 = start_out + (k+Nprime)*stride_out
        @muladd k2 = start_out + (k+2*Nprime)*stride_out
        y_k0, y_k1, y_k2 = out[k0], out[k1], out[k2]
        @muladd out[k0] = y_k0 + y_k1*wk1 + y_k2*wk2
        @muladd out[k1] = y_k0 + y_k1*wk1*plus120 + y_k2*wk2*minus120
        @muladd out[k2] = y_k0 + y_k1*wk1*minus120 + y_k2*wk2*plus120
        wk1 *= w1
        wk2 *= w2
    end
end

function fft!(out::AbstractVector{T}, in::AbstractVector{U}, start_out::Int, start_in::Int, d::Direction, ::Pow3FFT, g::CallGraph{T}, idx::Int) where {T,U}
    root = g[idx]
    N = root.sz
    s_in = root.s_in
    s_out = root.s_out
    p_120 = convert(T, cispi(2/3))
    m_120 = convert(T, cispi(4/3))
    if d == FFT_FORWARD
        fft_pow3!(out, in, N, start_out, s_out, start_in, s_in, d, m_120, p_120)
    else
        fft_pow3!(out, in, N, start_out, s_out, start_in, s_in, d, p_120, m_120)
    end
end


text:
BLAS.dot(a::Vector{Complex{T}}, b::Vector{Complex{T}}) where T = BLAS.dotc(a,b)

"""
    cg(A, x::Vector, b::Vector; iterations::Int = 30, relTol = 1.e-3
      , solverInfo = nothing )

Simple conjugate gradient algorithm.
The system matrix contained in AbstractLinearTrafo MUST be symmetrical and
and positive definite

# Arguments
* `A`                           - system matrix
* `x::Vector`                   - solution vector with the initial guess
* `b::Vector`                   - data vector (right hand side of equation)
* `iterations::Int = 30`        - maximum number of iterations
* `relTol = 1.e-3`              - stopping criteria for the relativ residual change
* `solverInfo = nothing`        - `solverInfo` object used to store convergence metrics
"""
function cg(A
            , x::Vector
            , b::Vector
            ; iterations::Int = 30
            , relTol = 1.e-3
            , solverInfo = nothing)

  r = b-A*x
  p = r
  rsold = BLAS.dot(r,r)
  rsnew = rsold
  r0= norm(b)
  iter_conv = iterations

  for i=1:iterations
    Ap = A*p
    bla = BLAS.dot(p,Ap)

    alpha = rsold/ bla
    #x = x+alpha*p;
    BLAS.axpy!(alpha,p,x)
    # BLAS.axpy!(a,X,Y) overwrites Y with a*X + Y
    BLAS.axpy!(-alpha,Ap,r)
    #r = r-alpha*Ap
    rsnew = BLAS.dot(r,r)
    if sqrt(abs(rsnew))/r0<relTol
      iter_conv = i
      break
    end

    beta = rsnew/rsold
    p = r+beta*p
    rsold = rsnew
  end

  if solverInfo != nothing
    storeInfo(solverInfo,A,b,x;residual=rsnew)
  end
  return x
end

"""
    cg(A, x::Vector, b::Vector, M; iterations::Int = 30, relTol = 1.e-3
      , solverInfo = nothing)

Preconditionned conjugate gradient algorithm.
The system matrix contained in AbstractLinearTrafo MUST be symmetrical and
and positive definite

The arguments are the same as for the non-preconditionned case.
`M` denotes the precondionner.
"""
function cg(A
            , x::Vector
            , b::Vector
            , M
            ; iterations::Int = 30
            , relTol = 1.e-3
            , solverInfo = nothing)

  r = b-A*x
  z = M*r
  p = z
  rsold = BLAS.dot(z,r)
  rsnew = rsold
  r0= norm(b)
  iter_conv = iterations

  for i=1:iterations
    Ap = A*p
    bla = BLAS.dot(p,Ap)

    alpha = rsold/ bla
    #x = x+alpha*p;
    BLAS.axpy!(alpha,p,x)
    # BLAS.axpy!(a,X,Y) overwrites Y with a*X + Y
    # r = r-alpha*Ap
    BLAS.axpy!(-alpha,Ap,r)
    z = M*r

    rsnew = BLAS.dot(z,r)
    if sqrt(abs(rsnew))/r0<relTol
      iter_conv = i
      break
    end

    beta = rsnew/rsold
    p = z+beta*p
    rsold = rsnew
  end

  if solverInfo != nothing
    storeInfo(solverInfo,A,b,x;residual=rsnew)
  end

  return x
end

ext: jl


text:
module EmDee

include("vec3.jl")
include("lennard_jones.jl")
include("nonbonded.jl")
include("molecular_graphs.jl")
include("modelling.jl")

end


text:
## Construction

_szerr(N::Integer, sz::Integer) = throw(ArgumentError(
    "cannot construct a SizedString{$N} from a $sz byte input"))

const UIntOrChar = Union{Unsigned, AbstractChar}

struct SizedString{N} <: AbstractString
    data::NTuple{N,UInt8}

    SizedString{N}(t::NTuple{M,UIntOrChar}) where {N, M} = N == M ? new{N}(t) : _szerr(N, M)
    function SizedString{N}(s::String) where {N}
        N == sizeof(s) ? new{N}(NTuple{N,UInt8}(s)) : _szerr(N, sizeof(s))
    end
end

SizedString{N}(cs::UIntOrChar...) where {N, M} = SizedString{N}(cs)
SizedString{N}(v::Vector{<:UIntOrChar}) where {N} = SizedString{N}(v...)
SizedString{N}(fs::SizedString{N}) where {N} = fs

# Fallback
SizedString{N}(s) where {N} = SizedString{N}(String(s))

macro sized_str(s)
    :(SizedString{length($(esc(s)))}($(esc(s))))
end

## AbstractString Implementation

ncodeunits(s::SizedString{N}) where {N} = N

sizeof(s::SizedString) = sizeof(s.data)

length(s::SizedString) = ncodeunits(s)
length(S::Type{SizedString{N}}) where {N} = N

lastindex(s::SizedString{N}) where {N} = N

function iterate(s::SizedString{N}, i::Int = 1) where {N}
    i > N && return nothing
    return Char(s.data[i]), i+1
end

codeunit(::SizedString) = UInt8
@propagate_inbounds codeunit(s::SizedString, i::Integer) = s.data[i]
@propagate_inbounds getindex(s::SizedString, i::Integer)::Char = s.data[i]

isvalid(s::SizedString, i::Integer) = checkbounds(Bool, s, i)

## IO

function read(io::IO, T::Type{SizedString{N}}) where N
    return read!(io, Ref{T}())[]::T
end

function write(io::IO, s::SizedString{N}) where N
    return write(io, Ref(s))
end


text:
"""
    Trues = Ones{Bool, N, Axes} where {N, Axes}

Lazy version of `trues` with axes.
Typically created using `Trues(dims)` or `Trues(dims...)`

# Example
```jldoctest
julia> Trues(1,3)
1×3 Ones{Bool,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}} = true

julia> Trues((2,3))
2×3 Ones{Bool,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}} = true
```
"""
const Trues = Ones{Bool, N, Axes} where {N, Axes}

""" `Falses = Zeros{Bool, N, Axes}` (see `Trues`) """
const Falses = Zeros{Bool, N, Axes} where {N, Axes}

# y[mask] = x when mask isa Trues (cf y[:] = x)

function Base.to_indices(A::AbstractArray{T,N}, inds, I::Tuple{Trues{N}}) where {T,N}
    @boundscheck axes(A) == axes(I[1]) || Base.throw_boundserror(A, I[1])
    (vec(LinearIndices(A)),)
end


text:
export KroneckerOperator

##########
# KroneckerOperator gives the kronecker product of two 1D operators
#########

struct KroneckerOperator{S,V,DS,RS,DI,RI,T} <: Operator{T}
    ops::Tuple{S,V}
    domainspace::DS
    rangespace::RS
    domaintensorizer::DI
    rangetensorizer::RI
end

KroneckerOperator(A,B,ds::Space,rs::Space,di,ri) =
    KroneckerOperator{typeof(A),typeof(B),typeof(ds),typeof(rs),typeof(di),typeof(ri),
                        promote_type(eltype(A),eltype(B))}((A,B),ds,rs,di,ri)

KroneckerOperator(A,B,ds::Space,rs::Space) = KroneckerOperator(A,B,ds,rs,
                    CachedIterator(tensorizer(ds)),CachedIterator(tensorizer(rs)))
function KroneckerOperator(A,B)
    ds=domainspace(A)⊗domainspace(B)
    rs=rangespace(A)⊗rangespace(B)
    KroneckerOperator(A,B,ds,rs)
end
KroneckerOperator(A::UniformScaling,B::UniformScaling) =
    KroneckerOperator(ConstantOperator(A.λ),ConstantOperator(B.λ))
KroneckerOperator(A,B::UniformScaling) = KroneckerOperator(A,ConstantOperator(B.λ))
KroneckerOperator(A::UniformScaling,B) = KroneckerOperator(ConstantOperator(A.λ),B)
KroneckerOperator(A::Fun,B::Fun) = KroneckerOperator(Multiplication(A),Multiplication(B))
KroneckerOperator(A::UniformScaling,B::Fun) = KroneckerOperator(ConstantOperator(A.λ),Multiplication(B))
KroneckerOperator(A::Fun,B::UniformScaling) = KroneckerOperator(Multiplication(A),ConstantOperator(B.λ))
KroneckerOperator(A,B::Fun) = KroneckerOperator(A,Multiplication(B))
KroneckerOperator(A::Fun,B) = KroneckerOperator(Multiplication(A),B)

function promotedomainspace(K::KroneckerOperator,ds::TensorSpace)
    A=promotedomainspace(K.ops[1],ds.spaces[1])
    B=promotedomainspace(K.ops[2],ds.spaces[2])
    KroneckerOperator(A,B,ds,rangespace(A)⊗rangespace(B))
end

function promoterangespace(K::KroneckerOperator,rs::TensorSpace)
    A=promoterangespace(K.ops[1],rs.spaces[1])
    B=promoterangespace(K.ops[2],rs.spaces[2])
    KroneckerOperator(A,B,domainspace(K),rs)
end

function convert(::Type{Operator{T}},K::KroneckerOperator) where T<:Number
    if T == eltype(K)
        K
    else
        ops=Operator{T}(K.ops[1]),Operator{T}(K.ops[2])
        KroneckerOperator{typeof(ops[1]),typeof(ops[2]),typeof(K.domainspace),typeof(K.rangespace),
                            typeof(K.domaintensorizer),typeof(K.rangetensorizer),T}(ops,
              K.domainspace,K.rangespace,
              K.domaintensorizer,K.rangetensorizer)
    end
end

function colstart(A::KroneckerOperator,k::Integer)
    K=block(A.domaintensorizer,k)
    M = blockbandwidth(A,2)
    if isfinite(M)
        blockstart(A.rangetensorizer,max(Block(1),K-M))
    else
        blockstart(A.rangetensorizer,Block(1))
    end
end

function colstop(A::KroneckerOperator,k::Integer)
    k == 0 && return 0
    K=block(A.domaintensorizer,k)
    st=blockstop(A.rangetensorizer,blockcolstop(A,K))
    # zero indicates above dimension
    min(size(A,1),st)
end

function rowstart(A::KroneckerOperator,k::Integer)
    K=block(rangespace(A),k)
    K2 = Int(K)-blockbandwidth(A,1)
    K2 ≤ 1 && return 1
    ds = domainspace(A)
    K2 ≥ blocksize(ds,1) && return size(A,2)
    blockstart(ds,K2)
end

function rowstop(A::KroneckerOperator,k::Integer)
    K=block(rangespace(A),k)
    ds = domainspace(A)
    K2 = Int(K)+blockbandwidth(A,2)
    K2 ≥ blocksize(ds) && return size(A,2)
    st=blockstop(ds,K2)
    # zero indicates above dimension
    st==0 ? size(A,2) : min(size(A,2),st)
end

bandwidths(K::KroneckerOperator) = (ℵ₀,ℵ₀)

isblockbanded(K::KroneckerOperator) = all(isblockbanded,K.ops)
isbandedblockbanded(K::KroneckerOperator) =
    all(op->isbanded(op) && isinf(size(op,1)) && isinf(size(op,2)),K.ops)
israggedbelow(K::KroneckerOperator) = all(israggedbelow,K.ops)

blockbandwidths(K::KroneckerOperator) =
    (blockbandwidth(K.ops[1],1)+blockbandwidth(K.ops[2],1),
    blockbandwidth(K.ops[1],2)+blockbandwidth(K.ops[2],2))

# If each block were in turn BlockBandedMatrix, these would
# be the    bandwidths
subblock_blockbandwidths(K::KroneckerOperator) =
    (max(blockbandwidth(K.ops[1],1),blockbandwidth(K.ops[2],2)) ,
           max(blockbandwidth(K.ops[1],2),blockbandwidth(K.ops[2],1)))

# If each block were in turn BandedMatrix, these are the bandwidths
function subblockbandwidths(K::KroneckerOperator)
    isbandedblockbanded(K) || return (ℵ₀,ℵ₀)

    if all(hastrivialblocks,domainspace(K).spaces) &&
            all(hastrivialblocks,rangespace(K).spaces)
        subblock_blockbandwidths(K)
    else
        dt = domaintensorizer(K).iterator
        rt = rangetensorizer(K).iterator
        # assume block size is repeated and square
        @assert all(b->isa(b,AbstractFill),dt.blocks)
        @assert rt.blocks ≡ dt.blocks

        sb = subblock_blockbandwidths(K)
        # divide by the size of each block
        sb_sz = mapreduce(getindex_value,*,dt.blocks)
        # spread by sub block szie
        (sb[1]+1)*sb_sz-1,(sb[2]+1)*sb_sz-1
    end
end

const Wrappers = Union{ConversionWrapper,MultiplicationWrapper,DerivativeWrapper,LaplacianWrapper,
                       SpaceOperator,ConstantTimesOperator}

isbandedblockbanded(P::Union{PlusOperator,TimesOperator}) = all(isbandedblockbanded,P.ops)

blockbandwidths(P::PlusOperator) = mapreduce(blockbandwidths, (a,b) -> max.(a,b), P.ops)
subblockbandwidths(P::PlusOperator) = mapreduce(subblockbandwidths, (a,b) -> max.(a,b), P.ops)

blockbandwidths(P::TimesOperator) = mapreduce(blockbandwidths, (a,b) -> a .+ b, P.ops)
subblockbandwidths(P::TimesOperator) = mapreduce(subblockbandwidths, (a,b) -> a .+ b, P.ops)

domaintensorizer(R::Operator) = tensorizer(domainspace(R))
rangetensorizer(R::Operator) = tensorizer(rangespace(R))

domaintensorizer(P::PlusOperator) = domaintensorizer(P.ops[1])
rangetensorizer(P::PlusOperator) = rangetensorizer(P.ops[1])

domaintensorizer(P::TimesOperator) = domaintensorizer(P.ops[end])
rangetensorizer(P::TimesOperator) = rangetensorizer(P.ops[1])

for FUNC in (:blockbandwidths, :subblockbandwidths, :isbandedblockbanded,:domaintensorizer,:rangetensorizer)
    @eval $FUNC(K::Wrappers) = $FUNC(K.op)
end

domainspace(K::KroneckerOperator) = K.domainspace
rangespace(K::KroneckerOperator) = K.rangespace

domaintensorizer(K::KroneckerOperator) = K.domaintensorizer
rangetensorizer(K::KroneckerOperator) = K.rangetensorizer

# we suport 4-indexing with KroneckerOperator
# If A is K x J and B is N x M, then w
# index to match KO=reshape(kron(B,A),N,K,M,J)
# that is
# KO[k,n,j,m] = A[k,j]*B[n,m]
# TODO: arbitrary number of ops

getindex(KO::KroneckerOperator,k::Integer,n::Integer,j::Integer,m::Integer) =
    KO.ops[1][k,j]*KO.ops[2][n,m]

function getindex(KO::KroneckerOperator,kin::Integer,jin::Integer)
    j,m=KO.domaintensorizer[jin]
    k,n=KO.rangetensorizer[kin]
    KO[k,n,j,m]
end

function getindex(KO::KroneckerOperator,k::Integer)
    if size(KO,1) == 1
        KO[1,k]
    elseif size(KO,2) == 1
        KO[k,1]
    else
        throw(ArgumentError("[k] only defined for 1 x ∞ and ∞ x 1 operators"))
    end
end

*(A::KroneckerOperator,B::KroneckerOperator) =
    KroneckerOperator(A.ops[1]*B.ops[1],A.ops[2]*B.ops[2])

## Shorthand

⊗(A,B) = kron(A,B)

Base.kron(A::Operator,B::Operator) = KroneckerOperator(A,B)
Base.kron(A::Operator,B) = KroneckerOperator(A,B)
Base.kron(A,B::Operator) = KroneckerOperator(A,B)
Base.kron(A::AbstractVector{T},B::Operator) where {T<:Operator} =
    Operator{promote_type(eltype(T),eltype(B))}[kron(a,B) for a in A]
Base.kron(A::Operator,B::AbstractVector{T}) where {T<:Operator} =
    Operator{promote_type(eltype(T),eltype(A))}[kron(A,b) for b in B]
Base.kron(A::AbstractVector{T},B::UniformScaling) where {T<:Operator} =
    Operator{promote_type(eltype(T),eltype(B))}[kron(a,1.0B) for a in A]
Base.kron(A::UniformScaling,B::AbstractVector{T}) where {T<:Operator} =
    Operator{promote_type(eltype(T),eltype(A))}[kron(1.0A,b) for b in B]

## transpose

Base.transpose(K::KroneckerOperator)=KroneckerOperator(K.ops[2],K.ops[1])

for TYP in (:ConversionWrapper,:MultiplicationWrapper,:DerivativeWrapper,:IntegralWrapper,:LaplacianWrapper),
    FUNC in (:domaintensorizer,:rangetensorizer)
    @eval $FUNC(S::$TYP) = $FUNC(S.op)
end

Base.transpose(S::SpaceOperator) =
    SpaceOperator(transpose(S.op), transpose(domainspace(S)), transpose(rangespace(S)))
Base.transpose(S::ConstantTimesOperator) = sp.c*transpose(S.op)

### Calculus

#TODO: general dimension
function Derivative(S::TensorSpace{SV,DD},order::Vector{Int}) where {SV,DD<:EuclideanDomain{2}}
    @assert length(order)==2
    if order[1]==0
        Dy=Derivative(S.spaces[2],order[2])
        K=Operator(I,S.spaces[1])⊗Dy
        T=eltype(Dy)
    elseif order[2]==0
        Dx=Derivative(S.spaces[1],order[1])
        K=Dx⊗Operator(I,S.spaces[2])
        T=eltype(Dx)
    else
        Dx=Derivative(S.spaces[1],order[1])
        Dy=Derivative(S.spaces[2],order[2])
        K=Dx⊗Dy
        T=promote_type(eltype(Dx),eltype(Dy))
    end
    # try to work around type inference
    DerivativeWrapper{typeof(K),typeof(domainspace(K)),Vector{Int},T}(K,order)
end

DefiniteIntegral(S::TensorSpace) = DefiniteIntegralWrapper(mapreduce(DefiniteIntegral,⊗,S.spaces))

### Copy

# finds block lengths for a subrange
blocklengthrange(rt, B::Block) = [blocklength(rt,B)]
blocklengthrange(rt, B::BlockRange) = blocklength(rt,B)
function blocklengthrange(rt, kr)
    KR=block(rt,first(kr)):block(rt,last(kr))
    Klengths=Vector{Int}(length(KR))
    for ν in eachindex(KR)
        Klengths[ν]=blocklength(rt,KR[ν])
    end
    Klengths[1]+=blockstart(rt,KR[1])-kr[1]
    Klengths[end]+=kr[end]-blockstop(rt,KR[end])
    Klengths
end

function bandedblockbanded_convert!(ret, S::SubOperator, KO, rt, dt)
    pinds = parentindices(S)
    kr,jr = pinds

    kr1,jr1 = reindex(S,pinds,(1,1))

    Kshft = block(rt,kr1)-1
    Jshft = block(dt,jr1)-1

    for J=blockaxes(ret,2)
        jshft = (J==Block(1) ? jr1 : blockstart(dt,J+Jshft)) - 1
        for K=blockcolrange(ret,J)
            Bs = view(ret,K,J)
            Bspinds = parentindices(Bs)
            kshft = (K==Block(1) ? kr1 : blockstart(rt,K+Kshft)) - 1
            for ξ=1:size(Bs,2),κ=colrange(Bs,ξ)
                Bs[κ,ξ] = S[reindex(Bs,Bspinds,(κ,ξ))...]
            end
        end
    end

    ret
end

function default_BandedBlockBandedMatrix(S)
    KO = parent(S)
    rt=rangespace(KO)
    dt=domainspace(KO)
    ret = BandedBlockBandedMatrix(Zeros, S)
    bandedblockbanded_convert!(ret, S, parent(S), rt, dt)
end

BandedBlockBandedMatrix(S::SubOperator) = default_BandedBlockBandedMatrix(S)

const Trivial2DTensorizer = CachedIterator{Tuple{Int,Int},
                                             TrivialTensorizer{2}}

# This routine is an efficient version of KroneckerOperator for the case of
# tensor product of trivial blocks

function BandedBlockBandedMatrix(S::SubOperator{T,KroneckerOperator{SS,V,DS,RS,
                                     Trivial2DTensorizer,Trivial2DTensorizer,T},
                                     Tuple{BlockRange1,BlockRange1}}) where {SS,V,DS,RS,T}
    KR,JR = parentindices(S)
    KR_i, JR_i = Int.(KR), Int.(JR)

    KO=parent(S)

    ret = BandedBlockBandedMatrix(Zeros, S)

    A,B = KO.ops

    AA = convert(BandedMatrix, view(A, Block(1):last(KR),Block(1):last(JR)))::BandedMatrix{eltype(S)}
    Al,Au = bandwidths(AA)
    BB = convert(BandedMatrix, view(B, Block(1):last(KR),Block(1):last(JR)))::BandedMatrix{eltype(S)}
    Bl,Bu = bandwidths(BB)
    λ,μ = subblockbandwidths(ret)

    for J in blockaxes(ret,2), K in blockcolrange(ret,J)
        n,m = KR_i[Int(K)],JR_i[Int(J)]
        Bs = view(ret, K, J)
        l = min(Al,Bu+n-m,λ)
        u = min(Au,Bl+m-n,μ)
        @inbounds for j=1:m, k=max(1,j-u):min(n,j+l)
            a = AA[k,j]
            b = BB[n-k+1,m-j+1]
            c = a*b
            inbands_setindex!(Bs,c,k,j)
        end
    end
    ret
end

convert(::Type{BandedBlockBandedMatrix}, S::SubOperator{T,KroneckerOperator{SS,V,DS,RS,
                                     Trivial2DTensorizer,Trivial2DTensorizer,T},
                                     Tuple{BlockRange1,BlockRange1}}) where {SS,V,DS,RS,T} =
    BandedBlockBandedMatrix(S)

## TensorSpace operators

## Conversion

conversion_rule(a::TensorSpace,b::TensorSpace) = conversion_type(a.spaces[1],b.spaces[1])⊗conversion_type(a.spaces[2],b.spaces[2])
maxspace(a::TensorSpace,b::TensorSpace) = maxspace(a.spaces[1],b.spaces[1])⊗maxspace(a.spaces[2],b.spaces[2])

# TODO: we explicetly state type to avoid type inference bug in 0.4

ConcreteConversion(a::BivariateSpace,b::BivariateSpace) =
    ConcreteConversion{typeof(a),typeof(b),
                        promote_type(prectype(a),prectype(b))}(a,b)

Conversion(a::TensorSpace,b::TensorSpace) = ConversionWrapper(promote_type(prectype(a),prectype(b)),
                KroneckerOperator(Conversion(a.spaces[1],b.spaces[1]),Conversion(a.spaces[2],b.spaces[2])))

function Multiplication(f::Fun{TS},S::TensorSpace) where {TS<:TensorSpace}
    lr=LowRankFun(f)
    ops=map(kron,map(a->Multiplication(a,S.spaces[1]),lr.A),map(a->Multiplication(a,S.spaces[2]),lr.B))
    MultiplicationWrapper(f,+(ops...))
end

## Functionals
Evaluation(sp::TensorSpace,x::Vec) = EvaluationWrapper(sp,x,zeros(Int,length(x)),⊗(map(Evaluation,sp.spaces,x)...))
Evaluation(sp::TensorSpace,x::Tuple) = Evaluation(sp,Vec(x...))

# it's faster to build the operators to the last b
function mul_coefficients(A::SubOperator{T,KKO,Tuple{UnitRange{Int},UnitRange{Int}}}, b) where {T,KKO<:KroneckerOperator}
    P = parent(A)
    kr,jr = parentindices(A)
    dt,rt = domaintensorizer(P),rangetensorizer(P)
    KR,JR = Block(1):block(rt,last(kr)),Block(1):block(dt,last(jr))
    M = P[KR,JR]
    M*pad(b, size(M,2))
end


text:
module VIEstimators

    using Random, LinearAlgebra, StatsFuns, Zygote, ZygoteRules, Distributions
    import ChainRulesCore: frule, DoesNotExist

    # Variational inference-related things.
    include("util.jl")
    include("gaussian.jl")
    include("elbo_estimators.jl")
end # module


text:
# # Foundations of Bilevel Programming: Example 3
# This example is from the book _Foundations of Bilevel Programming_ by Stephan
# Dempe, Chapter 3.7, Page 59. [url](https://www.springer.com/gp/book/9781402006319)

# Model of the problem

# First level
# ```math
# \min \sum_{i\in I} x_i - z,\\
# \notag s.t.\\
# y^a_i\geq 0, \forall i \in I,\\
# y^a_i\leq 1, \forall i \in I,\\
# y^b_i\geq 0, \forall i \in I,\\
# y^b_i\leq 1, \forall i \in I,\\
# y^a_i + y^b_i = 1, \forall i \in I,\\
# z\geq 0,\\
# z\leq 1,\\
# y^a_1 + y^a_2 + y^a_3 \geq z,\\
# -y^b_1 - y^b_4 + y^a_3 \geq z,\\
# y^b_7 - y^b_6 + y^a_4 \geq z,\\
# y^a_5 + y^a_6 + y^a_7 \geq z.\\
# ```
# Second level
# ```math
# \min -\sum_{i\in I}x_i,\\
# \notag s.t.\\
# \sum x_i \geq 0,\\
# x_i \leq y^a_i, \forall i\in I,\\
# x_i \leq y^b_i, \forall i\in I,\\
# I = \{1,...,7\}
# ```

using BilevelJuMP
using Ipopt
using JuMP
using Test

model = BilevelModel(Ipopt.Optimizer, mode = BilevelJuMP.ProductMode(1e-9))

# Global variables

I = 7 # maximum literals
clauses = [[1,2,3],[-1,-4,3],[7,-6,4],[5,6,7]]
atol = 1e-6
# First we need to create all of the variables in the upper and lower problems:

# Upper level variables
@variable(Upper(model), ya[i=1:I])
@variable(Upper(model), yb[i=1:I])
@variable(Upper(model), z)

#Lower level variables
@variable(Lower(model), x[i=1:I])

# Then we can add the objective and constraints of the upper problem:

# Upper level objecive function
@objective(Upper(model), Min, sum(x[i] for i in 1:I) - z)

# Upper level constraints
@constraint(Upper(model), ca, z <= 1)
@constraint(Upper(model), cb, z >= 0)
@constraint(Upper(model), c1[i=1:I], ya[i] >= 0)
@constraint(Upper(model), c2[i=1:I], ya[i] <= 1)
@constraint(Upper(model), c3[i=1:I], yb[i] >= 0)
@constraint(Upper(model), c4[i=1:I], yb[i] <= 1)
@constraint(Upper(model), c5[i=1:I], ya[i] + yb[i] == 1)

# for c in clauses
@constraint(Upper(model), cc[k in eachindex(clauses)], sum(i > 0 ? ya[i] : yb[-i] for i in clauses[k]) >= z)

# Followed by the objective and constraints of the lower problem:

# Lower objective function
@objective(Lower(model), Min, -sum(x[i] for i in 1:I))

# Lower constraints
@constraint(Lower(model), b1[i=1:I], x[i] >= 0)
@constraint(Lower(model), b2[i=1:I], x[i] <= ya[i])
@constraint(Lower(model), b3[i=1:I], x[i] <= yb[i])

# Initial Starting conditions

JuMP.set_start_value.(x, 0)
JuMP.set_start_value.(ya, 1)
JuMP.set_start_value.(yb, 0)
JuMP.set_start_value(z, 1)
for i in 1:I
    JuMP.set_dual_start_value.(b1, 0)
    JuMP.set_dual_start_value.(b2, 0)
    JuMP.set_dual_start_value.(b3, -1)
end

# Now we can solve the problem and verify the solution again that reported by
# Dempe.

optimize!(model)
primal_status(model)
termination_status(model)

# Auto testing
@test objective_value(model) ≈ -1 atol=atol
@test value.(x) ≈ zeros(I) atol=atol
@test value.(ya) ≈ ones(I) atol=atol
@test value.(yb) ≈ zeros(I) atol=atol
@test value(z) ≈ 1 atol=atol
# @show dual.(b1) #≈ 6 atol=atol
# @show dual.(b2) #≈ 2 atol=atol
# @show dual.(b3) #≈ 2 atol=atol


text:
# general array for testing input to SpecialFunctions.jl
x = 1e-5:0.01:50.0

### Tests for besselk0
k0_SpecialFunctions = SpecialFunctions.besselk.(0, x)

k0_64 = besselk0.(Float64.(x))
k0_32 = besselk0.(Float32.(x))

# make sure output types match input types
@test k0_64[1] isa Float64
@test k0_32[1] isa Float32

# test against SpecialFunctions.jl
@test k0_64 ≈ k0_SpecialFunctions
@test k0_32 ≈ k0_SpecialFunctions

### Tests for besselk0x
k0x_SpecialFunctions = SpecialFunctions.besselkx.(0, x)

k0x_64 = besselk0x.(Float64.(x))
k0x_32 = besselk0x.(Float32.(x))

# make sure output types match input types
@test k0x_64[1] isa Float64
@test k0x_32[1] isa Float32

# test against SpecialFunctions.jl
@test k0x_64 ≈ k0x_SpecialFunctions
@test k0x_32 ≈ k0x_SpecialFunctions

### Tests for besselk1
k1_SpecialFunctions = SpecialFunctions.besselk.(1, x)

k1_64 = besselk1.(Float64.(x))
k1_32 = besselk1.(Float32.(x))

# make sure output types match input types
@test k1_64[1] isa Float64
@test k1_32[1] isa Float32

# test against SpecialFunctions.jl
@test k1_64 ≈ k1_SpecialFunctions
@test k1_32 ≈ k1_SpecialFunctions

### Tests for besselk1x
k1x_SpecialFunctions = SpecialFunctions.besselkx.(1, x)

k1x_64 = besselk1x.(Float64.(x))
k1x_32 = besselk1x.(Float32.(x))

# make sure output types match input types
@test k1x_64[1] isa Float64
@test k1x_32[1] isa Float32

# test against SpecialFunctions.jl
@test k1x_64 ≈ k1x_SpecialFunctions
@test k1x_32 ≈ k1x_SpecialFunctions

### Tests for besselk
@test besselk(0, 2.0) == besselk0(2.0)
@test besselk(1, 2.0) == besselk1(2.0)

@test besselk(5, 8.0) ≈ SpecialFunctions.besselk(5, 8.0)
@test besselk(5, 88.0) ≈ SpecialFunctions.besselk(5, 88.0)

@test besselk(100, 3.9) ≈ SpecialFunctions.besselk(100, 3.9)
@test besselk(100, 234.0) ≈ SpecialFunctions.besselk(100, 234.0)

# test small arguments and order
m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:700.0]
@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x]

# test medium arguments and order
m = 30:200; x = 5.0:5.0:100.0
t = Float64.([besselk(m, x) for m in m, x in x])
@test t ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x]

# test large orders
m = 200:5:1000; x = 400.0:10.0:1200.0
t = Float64.([besselk(m, x) for m in m, x in x])
@test t ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x]

# Float 32 tests for aysmptotic expansion
m = 20:5:200; x = 5.0f0:2.0f0:400.0f0
t = [besselk(m, x) for m in m, x in x]
@test t[10] isa Float32
@test t ≈ Float32.([SpecialFunctions.besselk(m, x) for m in m, x in x])

# test for low values and medium orders
m = 20:5:50; x = [1f-3, 1f-2, 1f-1, 1f0, 1.5f0, 2.0f0, 4.0f0]
t = [besselk(m, x) for m in m, x in x]
@test t[5] isa Float32
@test t ≈ Float32.([SpecialFunctions.besselk(m, x) for m in m, x in x])

@test iszero(besselk(20, 1000.0))
#@test isinf(besselk(250, 5.0))

### Tests for besselkx
@test besselkx(0, 12.0) == besselk0x(12.0)
@test besselkx(1, 89.0) == besselk1x(89.0)

@test besselkx(15, 82.123) ≈ SpecialFunctions.besselk(15, 82.123)*exp(82.123)
@test besselkx(105, 182.123) ≈ SpecialFunctions.besselk(105, 182.123)*exp(182.123)


text:
# ---
# title: 102. Binary Tree Level Order Traversal
# id: problem102
# author: Indigo
# date: 2021-06-02
# difficulty: Medium
# categories: Tree, Breadth-first Search
# link: <https://leetcode.com/problems/binary-tree-level-order-traversal/description/>
# hidden: true
# ---
#
# Given a binary tree, return the _level order_ traversal of its nodes' values.
# (ie, from left to right, level by level).
#
# For example:
# Given binary tree `[3,9,20,null,null,15,7]`,
#
#
#
#         3
#        / \
#       9  20
#         /  \
#        15   7
#
#
# return its level order traversal as:
#
#
#
#     [
#       [3],
#       [9,20],
#       [15,7]
#     ]
#
#
#
## @lc code=start
using LeetCode

function level_order(root::Union{Nothing, TreeNode{Int}})
    isnothing(root) && return Int[]
    q = TreeNode{Int}[]
    q2 = TreeNode{Int}[]
    res = Vector{Int}[]
    push!(q, root)
    while !isempty(q)
        push!(res, Int[])
        while !isempty(q)
            fst = popfirst!(q)
            push!(res[end], fst.val)
            !isnothing(fst.left) && push!(q2, fst.left)
            !isnothing(fst.right) && push!(q2, fst.right)
        end
        q, q2 = q2, q
    end
    return res
end
## @lc code=end


text:
#new linear function
my_f(x,y) = 2x+3y


text:
"""
rounded_capacity_inequality!(cmp, demand, capacity, edge_tail, edge_head, edge_x)

Generate the rounded capacity cuts and put them in `cmp::CnstrMgrPointer`.

Input:
- Customers are numbered 2, 3, ... n_customers + 1.
- Depot is numbered 1.

Returns lists of `S` and `RHS` in the form of
`x(S:S) <= |S| - k(S)`

`return list_S, list_RHS, is_integer_and_feasible, max_violation`
"""
function rounded_capacity_cuts!(
    cut_manager::CutManager,
    demand::Vector{Int},
    capacity::Int,
    edge_tail::Vector{Int},
    edge_head::Vector{Int},
    edge_x::Vector{Float64};
    max_n_cuts = 1000::Cint,
    integrality_tolerance = 1e-4::Cdouble
)
    n_customers, _demand, n_edges, _edge_tail, _edge_head, _edge_x = input_conversion(demand, edge_tail, edge_head, edge_x)

    new_cmp = init_CMP()

    _integer_and_feasible = Ref{Cchar}()
    _max_violation = Ref{Cdouble}()

    ccall(
        (:CAPSEP_SeparateCapCuts, LIB_CVRPSEP),
        Cvoid,
        (
            Cint, Ptr{Cint}, Cint,
            Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
            CnstrMgrPointer, Cint, Cdouble,
            Ref{Cchar}, Ref{Cdouble}, CnstrMgrPointer
        ),
        Cint(n_customers), Cint.(_demand), Cint(capacity),
        Cint(n_edges), Cint.(_edge_tail), Cint.(_edge_head), Cdouble.(_edge_x),
        cut_manager.cmp, Cint(max_n_cuts), Cdouble(integrality_tolerance),
        _integer_and_feasible, _max_violation, new_cmp
    )

    cut_manager.new_cuts = convert_CMP(new_cmp)
    jl_CMGR_MoveCnstr!(new_cmp, cut_manager.cmp)

    cut_manager.is_integer_and_feasible = Bool(_integer_and_feasible[])
    cut_manager.violation = Float64(_max_violation[])

    S = map(x -> x.IntList, cut_manager.new_cuts)
    # then .+ 1 for Julia customer numbering:
    #       2, 3, ... n_customers + 1
    for s in S
        s .+= 1
    end
    RHS = map(x -> x.RHS, cut_manager.new_cuts)

    return S, RHS
end


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
dims = MPI.Dims_create(nnodes, [0,0,0])

@test prod(dims) == nnodes

comm_cart = MPI.Cart_create(comm, dims; periodic=[0, 1, 0], reorder=1)
@test MPI.Comm_size(comm_cart) == nnodes

comm_cart2 = MPI.Cart_create(comm, dims; periodic=(false, true, false), reorder=true)
@test MPI.Comm_size(comm_cart2) == nnodes

comm_sub1 = MPI.Cart_sub(comm_cart, (false, true, true))
@test MPI.Comm_size(comm_sub1) == div(nnodes, dims[1])

comm_sub2 = MPI.Cart_sub(comm_cart, (true, false, true))
@test MPI.Comm_size(comm_sub2) == div(nnodes, dims[2])

comm_sub3 = MPI.Cart_sub(comm_cart, (true, true, false))
@test MPI.Comm_size(comm_sub3) == div(nnodes, dims[3])

comm_sub4 = MPI.Cart_sub(comm_cart, Cint[true, true, false])
@test MPI.Comm_size(comm_sub4) == MPI.Comm_size(comm_sub3)

comm_cart = comm_cart2 = comm_sub1 = comm_sub2 = comm_sub3 = comm_sub4 = nothing

GC.gc()

MPI.Finalize()
@test MPI.Finalized()


text:
using DualNumbers

immutable DualQuaternion{T<:Real} <: Number
  q0::Quaternion{T}
  qe::Quaternion{T}
  norm::Bool
end

DualQuaternion( q0::Quaternion, qe::Quaternion, n::Bool = false) =
  DualQuaternion( promote( q0, qe )..., n )

DualQuaternion( d1::Dual, d2::Dual, d3::Dual, d4::Dual, n::Bool = false) =
  DualQuaternion( Quaternion( d1.re, d2.re, d3.re, d4.re, n ),
                  Quaternion( d1.du, d2.du, d3.du, d4.du ), n )

DualQuaternion( x::Real ) = DualQuaternion( Quaternion( x ), Quaternion( zero(x) ), abs(x) == one(x) )

DualQuaternion( d::Dual ) = DualQuaternion( d, zero(d), zero(d), zero(d), abs(d) == one(d.re) )

DualQuaternion( q::Quaternion ) = DualQuaternion( q, zero(q), q.norm )

DualQuaternion( a::Vector ) = DualQuaternion( zero(Quaternion{typeof(a[1])}), Quaternion( a ) )

convert{T}(::Type{DualQuaternion{T}}, x::Real) =
    DualQuaternion( convert( Quaternion{T}, x ), convert(Quaternion{T},0) )

convert{T}(::Type{DualQuaternion{T}}, d::Dual) =
    DualQuaternion( convert( Dual{T}, d ), convert( Dual{T}, 0 ), convert( Dual{T}, 0 ), convert( Dual{T}, 0 ) )

convert{T}(::Type{DualQuaternion{T}}, q::Quaternion) =
    DualQuaternion( convert( Quaternion{T}, q ), convert( Quaternion{T}, 0 ), q.norm )

convert{T<:Real}(::Type{DualQuaternion{T}}, q::DualQuaternion{T}) = q

convert{T}(::Type{DualQuaternion{T}}, dq::DualQuaternion) =
    DualQuaternion( convert( Quaternion{T}, dq.q0 ), convert( Quaternion{T}, dq.qe ), dq.norm )

promote_rule{T<:Real}(::Type{DualQuaternion{T}}, ::Type{T}) = DualQuaternion{T}
promote_rule{T<:Real}(::Type{DualQuaternion}, ::Type{T}) = DualQuaternion
promote_rule{T<:Real,S<:Real}(::Type{DualQuaternion{T}}, ::Type{S}) = DualQuaternion{promote_type(T,S)}
promote_rule{T<:Real,S<:Real}(::Type{Quaternion{T}}, ::Type{DualQuaternion{S}}) = DualQuaternion{promote_type(T,S)}
promote_rule{T<:Real,S<:Real}(::Type{DualQuaternion{T}}, ::Type{DualQuaternion{S}}) = DualQuaternion{promote_type(T,S)}

dualquat( q1, q2, n=false ) = DualQuaternion( q1, q2, n )
dualquat( d1, d2, d3, d4, n=false ) = DualQuaternion( d1, d2, d3, d4, n )
dualquat( x ) = DualQuaternion( x )

function show(io::IO, dq::DualQuaternion)
  show( io, dq.q0 )
  print( io, " + ( " )
  show( io, dq.qe )
  print( io, " )du" )
end

Q0( dq::DualQuaternion ) = dq.q0
Qe( dq::DualQuaternion ) = dq.qe

(/)(dq::DualQuaternion, x::Real) = DualQuaternion( dq.q0 / x, dq.qe / x )

(/)(dq::DualQuaternion, d::Dual) =
  DualQuaternion( dual( dq.q0.s , dq.qe.s  ) / d,
                  dual( dq.q0.v1, dq.qe.v1 ) / d,
                  dual( dq.q0.v2, dq.qe.v2 ) / d,
                  dual( dq.q0.v3, dq.qe.v3 ) / d )

abs2( dq::DualQuaternion ) = dq.norm ? dual( one( dq.q0.s ) ) :
  dual( abs2( dq.q0 ) ,
        2.0 * ( dq.q0.s  * dq.qe.s  +
                dq.q0.v1 * dq.qe.v1 +
                dq.q0.v2 * dq.qe.v2 +
                dq.q0.v3 * dq.qe.v3 ) )

abs( dq::DualQuaternion ) = dq.norm ? dual( one( dq.q0.s ) ) : sqrt( abs2( dq ) )

conj( dq::DualQuaternion ) = DualQuaternion( conj( dq.q0 ), conj( dq.qe ), dq.norm )
dconj( dq::DualQuaternion ) = DualQuaternion( dq.q0, -dq.qe, dq.norm )

inv( dq::DualQuaternion ) = dq.norm ? conj( dq ) : conj( dq ) / abs2( dq )

function normalize( dq::DualQuaternion )
  if ( dq.norm )
    return dq
  end
  a = abs( dq )
  if abs( a ) > 0
    qa = dq / a
    dualquat( qa.q0, qa.qe, true )
  else
    dq
  end
end

function normalizea( dq::DualQuaternion )
  if ( dq.norm )
    return ( dq, one( dual ) )
  end
  a = abs( dq )
  if abs( a ) > 0
    qa = dq / a
    dualquat( qa.q0, qa.qe, true ), a
  else
    dq, zero( dual )
  end
end

(-)(dq::DualQuaternion) = DualQuaternion( -dq.q0, -dq.qe, dq.norm)

(+)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion( dq.q0 + dw.q0, dq.qe + dw.qe )
(-)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion( dq.q0 - dw.q0, dq.qe - dw.qe )
(*)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion( dq.q0 * dw.q0,
                                                              dq.q0 * dw.qe + dq.qe * dw.q0,
                                                              dq.norm && dw.norm )
(/)(dq::DualQuaternion, dw::DualQuaternion) = dq*inv(dw)

function angleaxis( dq::DualQuaternion )
  tq = dq.qe * conj( dq.q0 )
  t = [ 2.0 * tq.v1, 2.0 * tq.v2, 2.0 * tq.v3 ]
  q0s = dq.q0.s
  th0, s0 = angleaxis( dq.q0 )
  sq0 = quat( 0.0, s0 )
  if abs( abs( q0s ) - one( q0s ) ) == 0
    th = dual( th0, 0.5 * abs( quat( 0, t ) ) )
    th, dualquat( sq0 )
  else
    th = dual( th0, 0.5 * dot( t, s0 ) )
    s0c1 = cross( s0, t )
    tanth = tan( th0 )
    s0c2 = ( s0c1 / tanth + t ) * 0.5
    sqev = cross( s0c2, s0 )
    th, dualquat( sq0, quat( 0.0, sqev ) )
  end
end

function angle( dq::DualQuaternion )
  th, ax = angleaxis( dq )
  th
end

function axis( dq::DualQuaternion )
  th, ax = angleaxis( dq )
  ax
end

function exp( dq::DualQuaternion )
  se = dual( dq.q0.s, dq.qe.s )
  se = exp( se )
  dq = dualquat( quat( 0.0, imag( dq.q0 ) ), quat( 0.0, imag( dq.qe ) ) )
  dq, th = normalizea( dq )
  if dq.norm
    dualquat( se ) * ( dualquat( cos( th ) ) + dq * dualquat( sin( th ) ) )
  else
    dualquat( se )
  end
end

function log( dq::DualQuaternion )
  dq, a = normalizea( dq )
  sl = log( a )
  th, s = angleaxis( dq )
  s * dualquat( th ) + dualquat( sl )
end

(^)(dq::DualQuaternion, dw::DualQuaternion) = exp( dw * log( dq ) )

function sqrt( dq::DualQuaternion )
  exp( 0.5 * log( dq ) )
end

dualquatrand() = dualquat( quatrand(), quatrand() )
ndualquatrand() = normalize( dualquatrand() )


text:
using ComponentArrays
using LinearAlgebra
import ReachabilityAnalysis.IntervalArithmetic
using Plots
import Polyhedra
using ReachabilityAnalysis
using Symbolics
using UnPack

## Constants and Symobics
const g = 9.80665
vars = @variables rx, ry, vx, vy, m_prop, Isp, m_dry, T, t

## General Utilities
# Convert from Intervals and numbers to LazySets
to_set(a::AbstractArray) = reduce(×, to_set.(a))
to_set(i::IntervalArithmetic.Interval) = Interval(i)
to_set(s::LazySet) = s
to_set(x) = Singleton([x])

# Make a ComponentArray of the layers of a neural network
function to_componentarray(network)
    layers = network.layers
    pairs = map(enumerate(layers)) do (i, layer)
        s = Symbol("L", i)
        s => (W=layer.weights, b=layer.bias)
    end
    return ComponentArray(NamedTuple(pairs))
end

function set_parameters!(network, values)
    for (i, key) in enumerate(valkeys(values))
        layer = values[key]
        network.layers[i].weights .= layer.W
        network.layers[i].bias .= layer.b
    end
end

## ODE Functions
function make_ode_functions(controls)
    @unpack kp, kd = controls

    function lander_thrust_on!(dx, x, p, t)
        # Unpack variables
        rx, ry, vx, vy, m_prop, Isp, m_dry, T, t = x

        # Intermediate parameters
        m = m_dry + m_prop
        θ = -kd*vx - kp*rx
        Fx, Fy = T .* (sin(θ), cos(θ))

        # State derivatives
        dx[1] = vx              # rx
        dx[2] = vy              # ry
        dx[3] = Fx/m            # vx
        dx[4] = -g + Fy/m       # vy
        dx[5] = -T/(Isp * g)    # m_prop
        dx[6:8] .= zero(x[1])   # Isp, m_dry, T
        dx[end] = one(x[1])      # t
        return nothing
    end

    function lander_thrust_off!(dx, x, p, t)
        # Unpack variables
        rx, ry, vx, vy, m_prop, Isp, m_dry, T, t = x
        _1 = one(x[1]) # Needed because numbers can't be converted into Taylor?

        # State derivatives
        dx .= zero(x[1])
        dx[1] = vx       # rx
        dx[2] = vy       # ry
        dx[4] = -_1*g    # vy
        dx[end] = _1     # t
        return nothing
    end

    function lander_terminated!(dx, x, p, t)
        dx .= zero(x[1])
        dx[end] = one(x[1])
        return nothing
    end

    return lander_thrust_on!, lander_thrust_off!, lander_terminated!
end

## Constraints
angle_sat = deg2rad(45)
# u_constraints = HPolyhedron([α ≥ -angle_sat, α ≤ angle_sat], controls)
above_ground = HPolyhedron([ry > 0, m_prop ≥ 0.1], vars)
burnt_out = HPolyhedron([ry > 0, m_prop ≤ 0.1], vars)
below_ground = HalfSpace(ry ≤ 0, vars)

## Modes
controls = ComponentArray(kp=0.01, kd=0.05)
f_on!, f_off!, f_terminated! = make_ode_functions(controls)

# Mode 1: Thrust on
thrust_on = @system(x'=f_on!(x), x ∈ above_ground, dims=9)

# Mode 2: Thrust off
thrust_off = @system(x'=f_off!(x), x ∈ burnt_out, dims=9)

# Mode 3: Terminated (hit ground)
terminated = @system(x'=f_terminated!(x), x ∈ below_ground, dims=9)

modes = [thrust_on, thrust_off, terminated]

## Mode Transitions
# State automoaton
automaton = LightAutomaton(3)

hits_ground = @map(x -> x, dim:9, x ∈ below_ground) #ConstrainedIdentityMap(length(vars), below_ground)
burns_out = @map(x -> x, dim:9, x ∈ HalfSpace(m_prop ≤ 0.1, vars)) #ConstrainedIdentityMap(length(vars), burnt_out)
resetmaps = Any[]

# Transition 1: thrust on -> thrust off
add_transition!(automaton, 1, 2, 1)
push!(resetmaps, burns_out)

# Transition 2: thrust on -> terminated
add_transition!(automaton, 1, 3, 2)
push!(resetmaps, hits_ground)

# Transition 3: thrust off -> terminated
add_transition!(automaton, 2, 3, 3)
push!(resetmaps, hits_ground)

## System Definition
H = HybridSystem(; automaton, modes, resetmaps=identity.(resetmaps))

## Initial Conditions
x0 = Any[
    0.0,        # rx
    120..150,   # ry
    0 ± 20,     # vx
    -45 ± 5,    # vy
    10.0,       # m_prop
    70 ± 3,     # Isp
    60 ± 5,     # m_dry
    900..950, # T
    0.0,        # t
] |> to_set |> concretize

ic = [(1, x0)]

prob = InitialValueProblem(H, ic)

##
boxdirs = BoxDirections(9)
alg = TMJets(; orderT=9, disjointness=BoxEnclosure(), maxsteps=10000)
@time sol = solve(prob, alg;
    T = 10.0,
    # max_jumps = 1,
    # # intersect_source_invariant = false,
    # # intersection_method = TemplateHullIntersection(boxdirs),
    # clustering_method = LazyClustering(1),
    # disjointness_method = BoxEnclosure(),
);
reach_sol = sol[1];

##
plt1 = plot(reach_sol; vars=(1,2), color=:goldenrod, label="Flowpipe", title="Position Set Propagation", xlabel="Cross Range (m)", ylabel="Altitude (m)", aspect_ratio=1)
plot!(reach_sol[1]; vars=(1,2), color=:red, label="Initial Uncertainty")
plot!(reach_sol[end]; vars=(1,2), color=:green, label="Final Uncertainty")

##
reqs = (-0.5..0.5) × (-10..0)
plot(reach_sol; vars=(3,4), label="Flowpipe", title="Velocity Set Propagation", xlabel="Horizontal Velocity (m/s)", ylabel="Vertical Velocity (m/s)", aspect_ratio=1)
plot!(reach_sol[1]; vars=(3,4), label="Initial Uncertainty")
plot!(reach_sol[end]; vars=(3,4), label="Final Uncertainty")
plot!(reqs, label="Requirements")

##
plot(sol[1], vars=(0,5))

##
plt = plot(1; xlims=xlims(plt1), ylims=ylims(plt1))
anim = @animate for u in reach_sol
    plot!(plt, u; vars=(1,2), legend=false, color=1, label="Flowpipe", title="Position Set Propagation", xlabel="Cross Range (m)", ylabel="Altitude (m)", aspect_ratio=1, xlims=xlims(plt1), ylims=ylims(plt1))
end
gf = gif(anim, "setprop_anim.gif", fps=7)
# plot(reach_sol; vars=(1,2), label="Flowpipe", title="Position Set Propagation", xlabel="Cross Range (m)", ylabel="Altitude (m)", aspect_ratio=1)
# plot!(reach_sol[1]; vars=(1,2), label="Initial Uncertainty")
# plot!(reach_sol[end]; vars=(1,2), label="Final Uncertainty")


text:
"""
Gage the geometric center (mean value) of an input array.

Arguments
---------
:Type{GeometricCenter}
    a Julia Type object and must be ``GeometricCenter``

array::Array
    a floating point array
"""
function gage(::Type{GeometricCenter}, array::Array)
    if issubtype(typeof(array[1]), AbstractArray)
        return gage_array_of_arrays(GeometricCenter, array)
    else
        if length(size(array)) == 1 # check array dimension
            return array
        else
            return mean(array, 1)
        end
    end
end


text:
include("../src/HartreeFock.jl")

using .HartreeFock
using BenchmarkTools

function H2O(r::Real, θ::Real)
    Molecule("H2O",
        [
            Atom("O", vec3(0, 0, 0)),
            Atom("H", vec3(r*cosd(θ/2),  r*sind(θ/2), 0)),
            Atom("H", vec3(r*cosd(θ/2), -r*sind(θ/2), 0))
        ]
    )
end

function H2O_pes(r_range::AbstractArray, θ_range::AbstractArray)
    Nr = length(r_range)
    Nt = length(θ_range)
    total_energy = zeros(Float64, Nr, Nt)
    Threads.@threads for i in 1:Nr
    # for i in 1:Nr
        r = r_range[i] / a0
        for j in  1:Nt
            θ = θ_range[j]
            e = hartree_fock(H2O(r, θ))
            total_energy[i, j] = e
        end
    end
    i, j = argmin(total_energy).I
    println("minimal energy of H2O is $(total_energy[i,j]*hartree)eV, with the bond length = $(r_range[i])Å, bond angle = $(θ_range[j])ᵒ")
end

# display(@benchmark H2O_pes(0.9:0.01:1.1, 85:115))
H2O_pes(0.9:0.01:1.1, 85:115)


text:
using DashBootstrapComponents, DashHtmlComponents

toast = html_div([
    dbc_button(
        "Open toast",
        id="simple-toast-toggle",
        color="primary",
        className="mb-3",
        n_clicks=0,
    ),
    dbc_toast(
        [html_p("This is the content of the toast", className="mb-0")],
        id="simple-toast",
        header="This is the header",
        icon="primary",
        dismissable=true,
    ),
]);

callback!(
    app,
    Output("simple-toast", "is_open"),
    Input("simple-toast-toggle", "n_clicks"),
) do n
    if n > 0
        return true
    end
    return no_update()
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module Registry

# for this module
import Pkg
import Pkg.TOML
using InteractiveUtils

# for testings decoding of metadata:
import MLJBase: Found, Continuous, Finite, Infinite
import MLJBase: OrderedFactor, Count, Multiclass, Binary

const srcdir = dirname(@__FILE__) # the directory containing this file
const environment_path = joinpath(srcdir, "..")

## METHODS TO GENERATE METADATA AND WRITE TO ARCHIVE

function finaltypes(T::Type)
    s = InteractiveUtils.subtypes(T)
    if isempty(s)
        return [T, ]
    else
        return reduce(vcat, [finaltypes(S) for S in s])
    end
end

const project_toml = joinpath(srcdir, "../Project.toml")
const packages = map(Symbol,
                     keys(TOML.parsefile(project_toml)["deps"])|>collect)
filter!(packages) do pkg
    !(pkg in (:InteractiveUtils, :Pkg, :MLJModelInterface))
end

const package_import_commands =  [:(import $pkg) for pkg in packages]

macro update()
    mod = __module__
    _update(mod, false)
end

macro update(ex)
    mod = __module__
    test_env_only = eval(ex)
    test_env_only isa Bool || "b in @update(b) must be Bool. "
    _update(mod, test_env_only)
end

function _update(mod, test_env_only)

    test_env_only && @info "Testing registry environment only. "

    program1 = quote
        @info "Packages to be searched for model implementations:"
        for pkg in $packages
            println(pkg)
        end
        using Pkg
        Pkg.activate($environment_path)
        @info "resolving registry environment..."
        Pkg.resolve()
    end

    program2 = quote

        @info "Instantiating registry environment..."
        Pkg.instantiate()

        @info "Loading registered packages..."
        import MLJBase
        import MLJModels
        using Pkg.TOML

        # import the packages
        $(Registry.package_import_commands...)

        @info "Generating model metadata..."

        modeltypes = MLJModels.Registry.finaltypes(MLJBase.Model)
        filter!(modeltypes) do T
            !(T in [MLJBase.ProbabilisticNetwork,
                    MLJBase.DeterministicNetwork,
                    MLJBase.UnsupervisedNetwork,
                    MLJBase.Interval,
                    MLJBase.Static])
        end

        # generate and write to file the model metadata:
        packages = string.(MLJModels.Registry.packages)
        meta_given_package = Dict()
        for pkg in packages
            meta_given_package[pkg] = Dict()
        end
        for M in modeltypes
            _info = MLJBase.info_dict(M)
            pkg = _info[:package_name]
            if !(pkg in ["unknown",])
                modelname = _info[:name]
                haskey(meta_given_package, pkg) ||
                    error("Bad `package_name` trait for $M: "*
                          "$pkg not a registered package. ")
                meta_given_package[pkg][modelname] = _info
                println(M, "\u2714 ")
            else
                println(M, "\u2717 (package_name is :unknown)")
            end
        end
        print("\r")

        open(joinpath(MLJModels.Registry.srcdir, "../Metadata.toml"), "w") do file
            TOML.print(file, MLJModels.encode_dic(meta_given_package))
        end

        # generate and write to file list of models for each package:
        models_given_pkg = Dict()
        for pkg in packages
            models_given_pkg[pkg] = collect(keys(meta_given_package[pkg]))
        end
        open(joinpath(MLJModels.Registry.srcdir, "../Models.toml"), "w") do file
            TOML.print(file, models_given_pkg)
        end

        :(println("Local Metadata.toml updated."))

    end

    mod.eval(program1)
    test_env_only || mod.eval(program2)

    println("You can now check the registry by running "*
            "`MLJModels.check_registry().\n\n"*
            "You can safely ignore \"conflicting import\" warnings. ")

    true
end

end # module


text:
export IncidenceList
export incidence_list

IncidenceVector = SparseVector

type IncidenceList
    nnodes::Int         # number of nodes
    ntimestamps::Int   # number of timestamps
    edges::Vector{IncidenceVector{Int}}
end

nodes(g::IncidenceList) = collect(1:g.nnodes)
num_nodes(g::IncidenceList) = g.nnodes
timestamps(g::IncidenceList) = collect(1:g.ntimestamps)
num_timestamps(g::IncidenceList) = g.ntimestamps
function edges(g::IncidenceList)
    nn = g.nnodes
    edges = g.edges
    elists = Tuple{Int, Int, Int}[]
    for e in edges
        v1, v2 = e.nzind
        i = mod(v1, nn)
        j = mod(v2, nn)
        t = round(Int, (v1-i)/nn + 1)
        push!(elists, (i,j,t)) # (i,j,t) node i to node j at time stamp t
    end
    elists
end

"""
`incidence_list(n)`

Create an incidence matrix with `n` nodes
"""
incidence_list(n::Int) = IncidenceList(n, 1, IncidenceVector{Int}[])

function add_edge!(g::IncidenceList, i::Int, j::Int, t::Int)
    nn = g.nnodes
    nt = g.ntimestamps
    if t > nt
        g.ntimestamps = t
    end
    i = i + (t-1)*nn
    j = j + (t-1)*nn
    i <= t*nn && j <= t*nn || throw(DimensionMismatch(""))

    push!(g.edges, sparsevec([i,j], [1,1], t*nn))
    g
end


text:
export VRGradient, LinearVRG, UniformVRG

################################################################################
# Variance Reduced Gradient Storage
################################################################################
abstract type AbstractVRGradient{N} end

Base.length(vrg::AbstractVRGradient{N}) where N = N

Base.firstindex(vrg::AbstractVRGradient) = 1
Base.lastindex(vrg::AbstractVRGradient) = length(vrg)
Base.iterate(vrg::AbstractVRGradient) = (i = firstindex(vrg); iterate(vrg,i))
Base.iterate(vrg::AbstractVRGradient, i) =
	i > lastindex(vrg) ? nothing : (vrg[i], i+1)

Base.getindex(vrg::AbstractVRGradient, i::Int) = not_implemented()

reset!(vrg::AbstractVRGradient) = not_implemented()

grad!(res, vrg::AbstractVRGradient, x, i) = not_implemented()
vrgrad!(res, vrg::AbstractVRGradient, x, i, innv_weight) = not_implemented()
innov!(res, vrg::AbstractVRGradient, x, i) = not_implemented()
grad_approx!(res, vrg::AbstractVRGradient) = not_implemented()

store_from_point!(vrg::AbstractVRGradient, x, ids) = not_implemented()
store_from_grad!(vrg::AbstractVRGradient, grad, state) = not_implemented()
store_from_innov!(vrg::AbstractVRGradient, innov, state) = not_implemented()

vrgrad_store!(res, vrg::AbstractVRGradient, x, i, innv_weight)=not_implemented()

not_implemented() = error("Not Implemented.")

# Convenience remappings
store_from_point!(vrg::AbstractVRGradient, x) =
	store_from_point!(vrg,x,1:length(vrg))

function grad!(res, vrg::AbstractVRGradient, x)
	grad!(res, vrg, x, 1)
	res_i = similar(res)
	for i = 2:length(vrg)
		grad!(res_i, vrg, x, i)
		res .+= res_i
	end
	res ./= length(vrg)
	return res
end

################################################################################
struct VRGradient{T<:Real, Vg<:AbstractArray{T}, GF, N} <: AbstractVRGradient{N}
	grad!::GF
	ysum::Vg
	y::Vector{Vg}

	function VRGradient(grad!, buf, n)
		ysum = similar(buf)
		y = [similar(buf) for _ = 1:n]
		vrg = new{eltype(buf), typeof(buf), typeof(grad!), n}(grad!, ysum, y)
		reset!(vrg)
		return vrg
	end
end

Base.getindex(vrg::VRGradient, i::Int) = vrg.y[i]

function reset!(vrg::VRGradient)
	vrg.ysum .= zero(eltype(vrg.ysum))
	for y in vrg.y
		y .= zero(eltype(y))
	end
end

function grad!(res, vrg::VRGradient, x, i)
	vrg.grad!(res, x, i)
	return i
end

function vrgrad!(res, vrg::VRGradient, x, i, innv_weight)
	grad!(res,vrg,x,i)
	res .-= vrg[i]
	res .*= innv_weight
	res .+= vrg.ysum./length(vrg)
end

function innov!(res, vrg::VRGradient, x, i)
	grad!(res,vrg,x,i)
	res .-= vrg[i]
	return i
end

grad_approx!(res, vrg::VRGradient) = (res .= vrg.ysum./length(vrg))

function store_from_point!(vrg::VRGradient, x, ids)
	for i in ids
		vrg.ysum .-= vrg[i]
		vrg.grad!(vrg[i], x, i)
		vrg.ysum .+= vrg[i]
	end
end

function store_from_grad!(vrg::VRGradient, grad, state)
	i = state
	vrg.ysum .-= vrg[i]
	vrg[i] .= grad
	vrg.ysum .+= grad
end

function store_from_innov!(vrg::VRGradient, innov, state)
	i = state
	vrg.ysum .+= innov
	vrg[i] .+= innov
end

function vrgrad_store!(res, vrg::VRGradient, x, i, innv_weight)
	grad_approx!(res, vrg)
	res .-= innv_weight.*vrg.y[i]
	vrg.ysum .-= vrg.y[i]
	grad!(vrg.y[i],vrg,x,i)
	res .+= innv_weight.*vrg.y[i]
	vrg.ysum .+= vrg.y[i]
end

################################################################################
struct LinearVRG{
		T<:Real, Vg<:AbstractArray{T}, D<:AbstractArray{T},
		Vd<:AbstractArray{D}, GF, N
		} <: AbstractVRGradient{N}
	grad::GF
	ysum::Vg
	y::Vector{T}
	data::Vd

	function LinearVRG(grad, data, buf)
		y = Vector{eltype(buf)}(undef, length(data))
		vrg = new{
				eltype(buf), typeof(buf), eltype(data),
				typeof(data), typeof(grad), length(data)
				}(
				grad, buf, y, data
				)
		reset!(vrg)
		return vrg
	end
end

Base.getindex(vrg::LinearVRG, i::Int) = vrg.data[i]*vrg.y[i]

function reset!(vrg::LinearVRG)
	vrg.ysum .= zero(eltype(vrg.ysum))
	vrg.y .= zero(eltype(vrg.y))
end

sgrad(vrg, x, i) = vrg.grad(dot(vrg.data[i],x),i)
function grad!(res, vrg::LinearVRG, x, i)
	sg = sgrad(vrg,x,i)
	res .= vrg.data[i].*sg
	return (sg,i)
end

function vrgrad!(res, vrg::LinearVRG, x, i, innv_weight)
	grad_approx!(res, vrg)
	innv = sgrad(vrg, x, i) - vrg.y[i]
	res .+= vrg.data[i].*innv.*innv_weight
end

function innov!(res, vrg::LinearVRG, x, i)
	innv = sgrad(vrg, x, i) - vrg.y[i]
	res .= vrg.data[i].*innv
	return (innv,i)
end

grad_approx!(res, vrg::LinearVRG) = (res .= vrg.ysum./length(vrg))

function store_from_point!(vrg::LinearVRG, x, ids)
	for i in ids
		innv = sgrad(vrg, x, i) - vrg.y[i]
		vrg.ysum .+= vrg.data[i].*innv
		vrg.y[i] += innv
	end
end

function store_from_grad!(vrg::LinearVRG, grad, state)
	(sg,i) = state
	innv = sg - vrg.y[i]
	vrg.ysum .+= vrg.data[i].*innv
	vrg.y[i] += innv
end

function store_from_innov!(vrg::LinearVRG, innov, state)
	(innov_s,i) = state
	vrg.ysum .+= innov
	vrg.y[i] += innov_s
end

function vrgrad_store!(res, vrg::LinearVRG, x, i, innv_weight)
	grad_approx!(res, vrg)
	innv = sgrad(vrg, x, i) - vrg.y[i]
	res .+= vrg.data[i].*innv.*innv_weight
	vrg.ysum .+= vrg.data[i].*innv
	vrg.y[i] += innv
end

################################################################################
struct UniformVRG{
		T<:Real, Vg<:AbstractArray{T}, GF, N
		} <: AbstractVRGradient{N}
	grad!::GF
	ymean::Vg
	xmem::Vg
	buf::Vg
	function UniformVRG(grad, buf, n)
		vrg = new{eltype(buf), typeof(buf), typeof(grad), n}(
			grad, buf, similar(buf), similar(buf))
		reset!(vrg)
		return vrg
	end
end

Base.getindex(vrg::UniformVRG, i::Int) =
	(vrg.grad!(vrg.buf,vrg.xmem,i); vrg.buf)

function reset!(vrg::UniformVRG)
	vrg.xmem .= zero(eltype(vrg.xmem))
	store_from_point!(vrg, vrg.xmem)
end

function grad!(res, vrg::UniformVRG, x, i)
	vrg.grad!(res, x, i)
	return nothing
end

function vrgrad!(res, vrg::UniformVRG, x, i, innv_weight)
	grad_approx!(res, vrg)
	vrg.grad!(vrg.buf, x, i)
	res .+= innv_weight.*vrg.buf
	vrg.grad!(vrg.buf, vrg.xmem, i)
	res .-= innv_weight.*vrg.buf
end

function innov!(res, vrg::UniformVRG, x, i)
	vrg.grad!(res, x, i)
	vrg.grad!(vrg.buf, vrg.xmem, i)
	res .-= vrg.buf
	return nothing
end

grad_approx!(res, vrg::UniformVRG) = (res .= vrg.ymean)

function store_from_point!(vrg::UniformVRG, x)
	vrg.xmem .= x
	vrg.grad!(vrg.ymean, x, 1)
	for i = 2:length(vrg)
		vrg.grad!(vrg.buf, x, i)
		vrg.ymean .+= vrg.buf
	end
	vrg.ymean ./= length(vrg)
end

invalid() = error("UniformVRG only support uniform store operations.")

store_from_point!(vrg::UniformVRG, x, ids) = invalid()
store_from_grad!(vrg::UniformVRG, grad, state) = invalid()
store_from_innov!(vrg::UniformVRG, innov, state) = invalid()
vrgrad_store!(res, vrg::UniformVRG, x, i, innv_weight) = invalid()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""CustomResourceValidation is a list of validation methods for CustomResources.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation(;
        openAPIV3Schema=nothing,
    )

    - openAPIV3Schema::IoK8sApiextensionsApiserverPkgApisApiextensionsV1JSONSchemaProps : openAPIV3Schema is the OpenAPI v3 schema to use for validation and pruning.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation <: SwaggerModel
    openAPIV3Schema::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1JSONSchemaProps } # spec name: openAPIV3Schema

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation(;openAPIV3Schema=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation, Symbol("openAPIV3Schema"), openAPIV3Schema)
        setfield!(o, Symbol("openAPIV3Schema"), openAPIV3Schema)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation = Dict{Symbol,Symbol}(Symbol("openAPIV3Schema")=>Symbol("openAPIV3Schema"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation = Dict{Symbol,String}(Symbol("openAPIV3Schema")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1JSONSchemaProps")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation }, name::Symbol, val)
end


text:
@testset "serial ego" begin
    # EI Min
    f = EfficientGlobalOptimization.Branin_Forrester()
    yₒ, xₒ = optimum(f)
    lb = lowerbounds(f)
    ub = upperbounds(f)
    x0 = sampling(lb, ub, 7)
    y0 = transpose(map(f, eachcol(x0)))
    krg = Kriging(x0, y0, MeanConst(0.0), SEArd([0.0,0.0], 3.0))
    kbd = (fill(-5.0, 3), fill(5.0, 3))
    mbd = ([-1.0], [1.0])
    tuner = KrigingTuner(kbd, mbd; global_every = 8, local_every = 2)
    ego = EGO(krg, tuner, lb, ub)
    try
        for i = 1:30
            x, fx = acquire(ego)
            y = f(x)
            append!(ego, x, y)
        end
    catch e
        if isa(e, NoFeasibleInfill)
            @info e.msg
        else
            rethrow(e)
        end
    end
    yo, xo = optimum(ego)
    @test yo[1] ≈ yₒ rtol = 0.001
end

@testset "parallel ego" begin

    f = EfficientGlobalOptimization.Branin_Forrester()
    yₒ, xₒ = optimum(f)
    lb = lowerbounds(f)
    ub = upperbounds(f)
    x0 = sampling(lb, ub, 7)
    y0 = transpose(map(f, eachcol(x0)))
    krg = Kriging(x0, y0, MeanConst(0.0), SEArd([0.0,0.0], 3.0))
    kbd = (fill(-5.0, 3), fill(5.0, 3))
    mbd = ([-1.0], [1.0])
    tuner = KrigingTuner(kbd, mbd; global_every = 8, local_every = 2)
    ego = EGO(krg, tuner, lb, ub)
    try
        for j = 1:6
            newx = []
            for i = 1:5
                x, fx = acquire(ego, newx...)
                push!(newx, x)
            end
            if isempty(newx)
                break
            end
            newy = map(f, newx)
            append!(ego, hcat(newx...), newy)
        end
    catch e
        if isa(e, NoFeasibleInfill)
            @info e.msg
        else
            rethrow(e)
        end
    end
    yo, xo = optimum(ego)
    @test yo[1] ≈ yₒ rtol = 0.001
end


text:
# partial port to Julia from Python -- python scoure:  https://github.com/worldveil/deuces/blob/master/deuces/card.py

module PokerCard

"""
We represent cards as 32-bit integers so there is no object instantiation -
they are just ints. Most of the bits are used and have a specific meaning.
See below:
                                PokerCard:
                    bitrank     suit rank   prime
              +--------+--------+--------+--------+
              |xxxbbbbb|bbbbbbbb|cdhsrrrr|xxpppppp|
              +--------+--------+--------+--------+

    1) p = prime number of rank (deuce = 2, trey = 3m four = 5, ..., ace = 41)
    2) r = rank of card (deuce = 0, trey = 1, four = 2, five = 3, ... ace = 12)
    3) cdhs = suit of the card (bit turned on based on suit of the card)
    4) b = bit it turned on depending on the rank of the card # 13 total
    5) x = unused bits

This will allow us to do very important things like:
    - Make a unique prime product for each hand
    - Detect flushes
    - Detect straights
    It is also quite performant.

"""

STR_RANKS = "23456789TJQKA"
INT_RANKS = 1:13
PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]

# conversion from string => int
CHAR_RANK_TO_INT_RANK = Dict(zip(STR_RANKS, INT_RANKS))
CHAR_SUIT_TO_INT_SUIT = Dict(
    's' => 1, # spades
    'h' => 2, # hearts
    'd' => 4, # diamonds
    'c' => 8, # clubs
)

INT_SUIT_TO_CHAR_SUIT = "xshxdxxxc"

# pretty printing
PRETTY_SUITS = Dict(
    1 =>  string(Char(0x2660)), # spades
    2 =>  string(Char(0x2764)),  # hearts
    4 =>  string(Char(0x2666)), # diamonds
    8 =>  string(Char(0x2663)) # clubs
)

function newCard(card_string::String)
    """
    Converts Card string to binary integer representation, inspired by:
    http: //www dot suffecool dot net/poker/evaluator dot html

    # do not go to the site, it contains malware these days allegedly
    """

    rank_char = card_string[1]
    suit_char = card_string[2]

    rank_int = CHAR_RANK_TO_INT_RANK[rank_char]
    suit_int = CHAR_SUIT_TO_INT_SUIT[suit_char]
    rank_prime = PRIMES[rank_int]

    # now we can compose the unique 32 bit int
    bitrank  = 1 << (rank_int-1) << 16
    suit = suit_int << 12
    rank = (rank_int) << 8

    return bitrank | suit | rank | rank_prime
end

function get_rank_int(card_int::Integer)
    return (card_int >> 8) & 0xF
end

function get_suit_int(card_int::Integer)
    return (card_int >> 12) & 0xF
end

function int_to_str(card_int::Integer)
    rank_int = get_rank_int(card_int)
    suit_int = get_suit_int(card_int)
    return STR_RANKS[rank_int] *  INT_SUIT_TO_CHAR_SUIT[suit_int+1]
end

function get_bitrank_int(card_int::Integer)
    return (card_int >> 16) & 0x1FFF
end

function hand_to_binary(cards_strs::Array)
    """
    Expects a list of cards as strings and returns a list
    of integers of same length corresponding to those strings.
    """

    hand_ints = []

    for card_str in cards_strs
        push!(hand_ints, newCard(card_str))
    end
    return hand_ints
end

function prime_product_from_hand(card_ints::Array)
    product = 1
    for c in card_ints
        product *= (c & 0xFF)
    end

    return product
end

function prime_product_from_rankbits(rankbits)
    """
    Returns the prime product using the bitrank (b)
    bits of the hand. Each 1 in the sequence is converted
    to the correct prime and multiplied in.
    Params:
        rankbits = a single 32-bit (only 13-bits set) integer representing
                the ranks of 5 _different_ ranked cards
                (5 of 13 bits are set)
    Primarily used for evaulating flushes and straights,
    two occasions where we know the ranks are *ALL* different.
    Assumes that the input is in form (set bits):
                          rankbits
                    +--------+--------+
                    |xxxbbbbb|bbbbbbbb|
                    +--------+--------+
    """

    product = 1
    for i in 0:12
        # if the ith bit is set
        # println(bitstring(rankbits%UInt32))
        val = rankbits & (1 << i)
        if sum(val) > 0
            product *= PRIMES[i+1]
        end
    end
    return product
end

function int_to_pretty_str(card_int::Integer)
    suit_int = get_suit_int(card_int)
    rank_int = get_rank_int(card_int)

    r = STR_RANKS[rank_int]
    s = PRETTY_SUITS[suit_int]

    return "" * r*s * ""
end

end # PokerCard

max_stars_repo_name: HigginTown/Quads

max_issues_repo_name: HigginTown/Quads

max_forks_repo_name: HigginTown/Quads


text:
using SimpleSegmentedStacks
using Test

function wastestack(call, n)
    if n == 0
        return n
    end
    return call() do
        wastestack(call, n - 1)
    end + n
end

@noinline naivecall(f) = f()

function wastestack_in_task(call, n)
    task = @task wastestack(call, n)
    yield(task)
    try
        return fetch(task)
    catch err
        if @isdefined TaskFailedException
            @assert err isa TaskFailedException
        end
        return task.result
    end
end

@testset begin
    @test wastestack_in_task(SimpleSegmentedStacks.call, 2^18) isa Int
    @test wastestack(SimpleSegmentedStacks.call, 2^18) isa Int
    @test wastestack_in_task(naivecall, 2^18) isa StackOverflowError
end


text:
module Phases

# Imports
import WaterData
import GeometricalPredicates, VoronoiDelaunay, Dierckx # JLD issues
using BasicUnits
  const bar = BasicUnits.bar
  const cm = BasicUnits.cm
using DataFrames
using Ogre
import GeometricalPredicates: geta, getb, getc, getx, gety
import WaterData: lognorm12, unlognorm, findtriangle, barycoords, Polygon, BoundingBox
using RecipesBase, PlotUtils, LaTeXStrings
using Parameters

# Exports
export findphase, PlanetPhases

# EOS
piecewise = WaterData.load_piecewise_eoses()
fe = piecewise["fe"]
mgsio3 = piecewise["mgsio3"]
h2o = WaterData.load_full_eos()["gridPlusIdeal"]

@enum ModelType InteriorOnly WholePlanet

# Planet structure
# include with default values 70% core, core = 1/3 iron, 2/3 silicate
@with_kw type PlanetModel
  M::Float64 = M_earth
  f_iron::Float64 = 0.7/3
  f_silicate::Float64 = 0.7 - f_iron
  f_water::Float64 = 1 - f_iron - f_silicate
  Tirr::Float64 = 300.
  ɛ::Float64 = 1e-14
  κ = nothing
  γ::Float64 = 0.01
  Psurf::Float64 = 100bar
  model_type::ModelType = WholePlanet
end

# distinguishing the core & water layers
# here, "core" refers to both iron+silicate
isiron(M, Miron) = (M < Miron)
issilicate(M, Mcore, Miron) = (Miron <= M < Mcore)
iswater(M, Mcore) = (M >= Mcore)
function findphase_cored(M, Mcore, Miron, P, T)
  if iswater(M, Mcore)
    findphase(P, T)
  elseif issilicate(M, Mcore, Miron)
    Silicate
  elseif isiron(M, Mcore)
    Iron
  else
    error("Couldn't get phase")
  end
end
function structure(pm::PlanetModel)
  @unpack M, f_iron, f_silicate, Tirr, ɛ, κ, γ, Psurf, model_type = pm
  if model_type == WholePlanet
    structure, radius = Ogre.interior(M, f_iron, f_silicate, ɛ, Tirr, κ, γ, Psurf)
  elseif model_type == InteriorOnly
    structure, radius = error("Interior-only models are not *quite* there yet")
  end
  Ts = Ogre.temperature(structure) |> reverse
  Ps = Ogre.pressure(structure) |> reverse
  Rs = Ogre.radius(structure) |> reverse
  Ms = Ogre.mass(structure) |> reverse
  clamp!(Rs, 0, Inf)

  return Ms, Rs, Ps, Ts
end
function radius(pm::PlanetModel)
  Ms, Rs, Ps, Ts = structure(pm)
  return maximum(Rs)
end
function phases(pm::PlanetModel)
  Ms, Rs, Ps, Ts = structure(pm)
  Mcore = pm.M * (pm.f_iron + pm.f_silicate)
  Miron = pm.M * pm.f_iron
  phases = map(Ms, Ps, Ts) do M, P, T
    findphase_cored(M, Mcore, Miron, P, T)
  end
  return PlanetPhases(Rs, phases)
end
function finddensity_cored(M, Mcore, Miron, P, T)
  if iswater(M, Mcore)
    h2o(P, T)
  elseif issilicate(M, Mcore, Miron)
    mgsio3(P, T)
  elseif isiron(M, Mcore)
    fe(P, T)
  else
    error("Couldn't get density")
  end
end
function densities(pm::PlanetModel)
  Ms, Rs, Ps, Ts = structure(pm)
  Mcore = pm.M * (pm.f_iron + pm.f_silicate)
  Miron = pm.M * pm.f_iron
  densities = map(Ms, Ps, Ts) do M, P, T
    finddensity_cored(M, Mcore, Miron, P, T)
  end

  return PlanetInterior(Rs, densities)
end

"Phases of water (or core material)"
@enum Phase Gas Liquid SupercriticalFluid I II III V VI VII VIII X Plasma Superionic Unknown Silicate Iron Nothing
defaultcolor = :Red
"Mappings from phases to plot colors (pretty)"
const phase_colors = Dict(
  Gas => :AliceBlue,
  Liquid => :SpringGreen,
  SupercriticalFluid => :DeepSkyBlue,
  Nothing => :White,
  I => defaultcolor,
  II => defaultcolor,
  III => defaultcolor,
  V => defaultcolor,
  VI => defaultcolor,
  VII => :SeaGreen,
  VIII => defaultcolor,
  X => :DarkSlateBlue,
  Plasma => :Coral,
  Superionic => :Gold,
  Unknown => :Fuchsia,
  Silicate => "#4A291B",
  Iron => :Black )

"Mappings from strings to a Phase"
const phasemappings = Dict(
  "idealgas"=>Gas,
  "iapws"=>Unknown, "french"=>Unknown, "fallback"=>Unknown,
  "fluid"=>SupercriticalFluid,
  "I"=>I,
  "II"=>II,
  "III"=>III,
  "V"=>V,
  "VI"=>VI,
  "VII"=>VII, "ice VII"=>VII, "sugimura"=>VII,
  "VIII"=>VIII,
  "X"=>X, "ice X"=>X,
  "superionic"=>Superionic,
  "plasma"=>Plasma )
Base.convert(::Type{Phase}, s::String) = phasemappings[s]

"Critical pressure of water"
Pcrit = WaterData.Pc
"Critical temperature of water"
Tcrit = WaterData.Tc
# regions for water phases
dummy_eos = WaterData.ConstantEOS(1.)
liquid = let
  bounds = WaterData.BoundingBox(0, Inf, 0, Tcrit)
WaterData.BoundedEOS(dummy_eos, bounds)
end
gas = let
bounds = WaterData.BoundingBox(0, Pcrit, Tcrit, Inf)
WaterData.BoundedEOS(dummy_eos, bounds)
end
supercritical = let
bounds = WaterData.BoundingBox(Pcrit, Inf, Tcrit, Inf)
WaterData.BoundedEOS(dummy_eos, bounds)
end

"EOSes with defined regions and associated phases"
immutable EOSCollectionPhases
  eoses::Vector{WaterData.EOS}
  phases::Vector{Phase}
end
eostables = let
  eosdata = WaterData.load_full_eos()
  eosraw = eosdata["raw"]
  eosgrid = eosdata["grid"]
  idealgas = eosdata["gridPlusIdeal"].eoses[1]

  eoses = copy(eosraw.eoses)[[2, 4, 5, 6, 7, 8, 10, 11]]

  eos_strings = ["sugimura", "I", "III", "V", "VI", "VII", "X", "X"]
  eos_phases = [phasemappings[p] for p in eos_strings]

  # add ideal gas
  push!(eoses, idealgas)
  push!(eos_phases, Gas)

  # add other phases
  push!(eoses, liquid)
  push!(eos_phases, Liquid)
  push!(eoses, supercritical)
  push!(eos_phases, SupercriticalFluid)
  push!(eoses, gas)
  push!(eos_phases, Gas)

  EOSCollectionPhases(eoses, eos_phases)
end
function Base.getindex(c::EOSCollectionPhases, P, T)::Phase
  for (eos, phase) in zip(c.eoses, c.phases)
    if (P, T) in eos
      return phase
    end
  end
  return Unknown
end

"An unstructured EOS table with associated phases for each point"
immutable TablePhases <: WaterData.TabularEOS
  P::Vector{Float64}
  T::Vector{Float64}
  ρ::Vector{Float64}
  phase::Vector{String}
  tess::WaterData.Tessellation

  function TablePhases(P, T, ρ, phase)
    u = WaterData.unique_indices(zip(P,T))
    new(P[u], T[u], ρ[u], phase[u], WaterData.get_tessellation(P[u], T[u], uselog=true))
  end
end
function Base.in(P, T, eos::TablePhases)
  if !in(P, T, BoundingBox(eos))
    return false
  else
    xn, yn = lognorm12(P, T, eos)
    return in(xn, yn, eos.tess)
  end
end
function closestpoint(p::TablePhases, P, T)
  Prange = extrema(p.P)
  Trange = extrema(p.T)
  PP = lognorm12(P, Prange...)
  TT = lognorm12(T, Trange...)
  tri = findtriangle(p.tess, PP, TT)
  i = indmax(barycoords(tri, PP, TT))
  closest = if i == 1
    geta(tri)
  elseif i == 2
    getb(tri)
  elseif i == 3
    getc(tri)
  end

  Pclose = unlognorm(getx(closest), Prange...)
  Tclose = unlognorm(gety(closest), Trange...)

  return (Pclose, Tclose)
end
function Base.getindex(p::TablePhases, P, T)::Phase
  if (P, T) in p
    (Pclose, Tclose) = closestpoint(p, P, T)
    for (Pi, Ti, phase) in zip(p.P, p.T, p.phase)
      if Pi ≈ Pclose && Ti ≈ Tclose
        return phase
      end
    end
  end
  return Unknown
end

"The French EOS and its associated phases"
frenchphases = let
  df = readtable(WaterData.config.rawdata * "/French.eos", allowcomments=true)
  P = collect(df[:pressure]) * bar * 1e3
  T = collect(df[:temperature]) * K
  ρ = collect(df[:density]) * g/cm^3
  phase = collect(df[:state])

  TablePhases(P, T, ρ, phase)
end

# run through all phase sources to find the appropriate phase code
function findphase(P, T)::Phase
  thephase = Unknown
  allphases = [frenchphases, eostables]
  for p in allphases
    if p[P, T] != Unknown
      thephase = p[P, T]
      break
    end
  end
  return thephase
end

function change_indices(x)
  where_changed = (x[1:end-1] .!= x[2:end])
  push!(where_changed, true) # ending counts as a change
  find(where_changed)
end

type PlanetPhases
  rs::Vector{Float64}
  ϕs::Vector{Phase}
end

type PlanetInterior
  rs::Vector{Float64}
  xs::Vector{Float64}
end

@recipe function plot(p::PlanetPhases)
  rs, ϕs = p.rs, p.ϕs
  unshift!(ϕs, ϕs[1])
  unshift!(rs, 0)
  colors = Symbol[]
  for (p, c) in phase_colors
    unshift!(rs, 0)
    unshift!(ϕs, p)
    unshift!(colors, c)
  end

  Θs = linspace(0, 2π)
  phasegrid = repeat(ϕs, inner=(1, length(Θs)))
  phasegrid_str = string.(phasegrid)
  Rchanges = rs[change_indices(ϕs)]
  Rp = maximum(rs)
  Rencl = ceil(Rp/R_earth)

  legend := false
  proj := :polar
  border := nothing
  yaxis := nothing
  # colorbar := false
  xaxis := nothing
  ylims --> (0, Rencl)

  # the phase structure
  @series begin
    seriestype := :heatmap
    color := cgrad(colors)

    Θs, rs/R_earth, phasegrid_str
  end

  # the borders
  @series begin
    seriestype := :path
    color := :black
    linewidth := 0

    Θs, Rchanges'/R_earth
  end

  # the outer border
  @series begin
    seriestype := :path
    color := :black
    linewidth := 0.5

    Θs, [Rp/R_earth]
  end

  # the inner grid
  @series begin
    seriestype := :path
    color := :grey
    linestyle := :dot
    linewidth := 0.5

    Θs, collect(1:Rencl-1)'
  end

  # the outer grid
  @series begin
    seriestype := :path
    color := :grey
    linestyle := :dot
    linewidth := 0.5

    Θs, collect(Rencl:10)'
  end
end

end # module Phases


text:
using DataFrames
import DataFrames: _names
import Compat.String

bestmime(io, x) =
    mimewritable(MIME"text/html"(), x) ? writemime(io, MIME"text/html"(), x) : write(io, string(x))

function Base.writemime(io::IO,
                        ::MIME"text/html",
                        df::DataFrame)
    n = size(df, 1)
    cnames = _names(df)
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for column_name in cnames
        write(io, "$column_name	")
    end
    write(io, "
")
    tty_rows, tty_cols = Base.tty_size()
    for row in 1:min(n, tty_rows)
        write(io, "
")
        write(io, "$row	")
        for column_name in cnames
            cell = df[row, column_name]
            write(io, "")
            bestmime(io, cell)
            write(io, "	")
        end
        write(io, "
")
    end
    if n > 20
        write(io, "
")
        write(io, "⋮	")
        for column_name in cnames
            write(io, "⋮	")
        end
        write(io, "
")
    end
    write(io, "
")

end
immutable Colored
    color::String
    value
end

Base.writemime(io::IO, ::MIME"text/html", x::Colored) =
    write(io, """$(x.value)""")

color_incomplete(mx, x, flip=false) = (x < mx) $ flip ? Colored("#a35", x) : Colored("#374", x)

tryfloat(x, d=0.0) = try float(x) catch d end
function color_progress(d)
    df = DataFrame(Names = d[:Names])
    fst = true
    for (n, c) in eachcol(d)
        if fst
            fst = false
            continue
        end
        mx = c[1]
        df[n] = [Colored("black", mx), map(x -> color_incomplete(tryfloat(mx), tryfloat(x)), c[2:end])]
    end
    df
end


text:
__precompile__()
module Roots
using Compat

export findroot, AbstractRootsMethod,
     NewtonMethod, InterpolationMethod

abstract type AbstractRootsMethod end

include("newton.jl")
include("interpolation.jl")

"""
    findroot(f, x₀, m=NewtonMethod())

Find a zero of `f` starting from the point `x₀`.

Available algorithms are `NewtonMethod` and  `InterpolationMethod` (necessitates `gnuplot`).

The derivative of `f` is computed by numerical discretization. Multivariate
functions are supported by `NewtonMethod`.

Returns a tuple `(ok, x, it, normf)`.

**Usage Example**
```
ok, x, it, normf = findroot(x->exp(x)-x^4, 1., NewtonMethod(atol=1e-5))
```
"""
findroot(f, x0) = findroot(f, x0, NewtonMethod())
findroot(f, x0::Real, m::AbstractRootsMethod) = findroot(f, Float64(x0), m)
findroot(f, x0::Float64, m::NewtonMethod) = newton(f, x0, m)
findroot(f, x0::Vector{Float64}, m::NewtonMethod) = newton(f, x0, m)
findroot(f, x0::Float64, m::InterpolationMethod) = findzero_interp(f, x0, m)

end #module

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
using Base.Test

using MDPs

function test1{P,R}(::Type{P}, ::Type{R}, s::Int, a::Int)
  transition = rand(P, s, s, a)
  # Make each row of ``P`` sum to 1
  transition = transition ./ sum(transition, 1)
  # Generate ``a`` vectors of lenth ``s`` of rewards between -1 and 1
  reward = one(R) .- 2rand(R, s, a)
  (transition, reward)
end
for T in (Float32, Float64)
  P, R = test1(T, T, 100, 3)
  @test is_square_stochastic(P)
end

function test2{P,R}(::Type{P}, ::Type{R}, s::Int, a::Int)
  transition = rand(P, s, s, a)
  # Make each row of ``P`` sum to 1
  for k = 1:a
    for j = 1:s
      S = zero(P)
      @simd for i = 1:s
        @inbounds S += transition[i, j, k]
      end
      @simd for i = 1:s
        @inbounds transition[i, j, k] /= S
      end
    end
  end
  # Generate ``a`` vectors of lenth ``s`` of rewards between -1 and 1
  reward = one(R) .- 2rand(R, s, a)
  (transition, reward)
end
for T in (Float32, Float64)
  P, R = test2(T, T, 100, 3)
  @test is_square_stochastic(P)
end

function test3{P,R}(::Type{P}, ::Type{R}, s::Int, a::Int)
  transition = rand(P, s, s, a)
  reward = rand(R, s, a)
  for k = 1:a
    for j = 1:s
      S = zero(P)
      @simd for i = 1:s
        @inbounds S += transition[i, j, k]
      end
      @simd for i = 1:s
        @inbounds transition[i, j, k] /= S
      end
      @inbounds reward[j, k] = one(R) - 2reward[j, k]
    end
  end
  (transition, reward)
end
for T in (Float32, Float64)
  P, R = test3(T, T, 100, 3)
  @test is_square_stochastic(P)
end

function test4{P,R}(::Type{P}, ::Type{R}, s::Int, a::Int)
  transition = Array(P, s, s, a)
  reward = Array(R, s, a)
  for k = 1:a
    for j = 1:s
      S = zero(P)
      @simd for i = 1:s
        r = rand(P)
        @inbounds transition[i, j, k] = r
        S += r
      end
      @simd for i = 1:s
        @inbounds transition[i, j, k] /= S
      end
      @inbounds reward[j, k] = one(R) - convert(R, 2)*rand(R)
    end
  end
  (transition, reward)
end
for T in (Float32, Float64)
  P, R = test4(T, T, 100, 3)
  @test is_square_stochastic(P)
end

print("Test 1: ")
@time for x = 1:10000
  test1(Float64, Float64, 100, 3)
end

print("Test 2: ")
@time for x = 1:10000
  test2(Float64, Float64, 100, 3)
end

print("Test 3: ")
@time for x = 1:10000
  test3(Float64, Float64, 100, 3)
end

print("Test 4: ")
@time for x = 1:10000
  test4(Float64, Float64, 100, 3)
end


text:
#=

# 301: 3D Laplace equation
([source code](SOURCE_URL))

=#

module Example301_Laplace3D

using VoronoiFVM,ExtendableGrids
using GridVisualize

## Flux function which describes the flux
## between neigboring control volumes
function g!(f,u0,edge)
    u=unknowns(edge,u0)
    f[1]=u[1,1]-u[1,2]
end

function main(;Plotter=nothing)
    nspecies=1
    ispec=1
    X=collect(0:0.2:1)
    grid=VoronoiFVM.Grid(X,X,X)
    physics=VoronoiFVM.Physics(num_species=nspecies,flux=g!)
    sys=VoronoiFVM.System(grid,physics)
    enable_species!(sys,ispec,[1])
    boundary_dirichlet!(sys,ispec,5,0.0)
    boundary_dirichlet!(sys,ispec,6,1.0)
    inival=unknowns(sys,inival=0)
    solution=unknowns(sys)
    solve!(solution,inival,sys)
    scalarplot(grid,solution[1,:],Plotter=Plotter,zplane=0.5, flevel=0.5)
    return solution[43]
end

## Called by unit test

function test()
    main() ≈ 0.2
end

end