text:
export RationalFunction
export poles, residues
export lowest_terms

## ----
"""
    AbstractRationalFunction{T,X,P}

Abstract type for holding ratios of polynomials of type `P{T,X}`.

Default methods for basic arithmetic operations are provided.

Numeric methods to cancel common factors, compute the poles, and return the residues are provided.

!!! Note:
    Requires `VERSION >= v"1.2.0"`
"""
abstract type AbstractRationalFunction{T,X,P} end

function Base.show(io::IO, pq::AbstractRationalFunction)
    p,q = pqs(pq)
    print(io,"(")
    printpoly(io, p)
    print(io, ") // (")
    printpoly(io, q)
    print(io, ")")
end

## helper to make a rational function of given type
function rational_function(::Type{R}, p::P, q::Q) where {R<:AbstractRationalFunction,
                                                         T,X,   P<:AbstractPolynomial{T,X},
                                                         S,   Q<:AbstractPolynomial{S,X}}
    constructorof(R)(promote(p,q)...)
end

## ---- conversion
function Base.convert(::Type{PQ}, pq′::PQ′) where {T,X,P,PQ <: AbstractRationalFunction{T,X,P},
                                                   T′,X′,P′,PQ′<:AbstractRationalFunction{T′,X′,P′} }
    !isconstant(pq′) && assert_same_variable(X,X′)
    p′,q′=pqs(pq′)
    𝑷 = isconstant(pq′) ? P :  promote_type(P, P′)
    p,q = convert(𝑷, p′), convert(𝑷, q′)
    rational_function(PQ, p, q)
end

# R ⊂ R[x] ⊂ R(x)
function Base.convert(::Type{PQ}, p::Number) where {PQ <: AbstractRationalFunction}
   P = eltype(PQ)
   rational_function(PQ, p * one(P), one(P))
end

function Base.convert(::Type{PQ}, q::LaurentPolynomial) where {PQ <: AbstractRationalFunction}
    m = firstindex(q)
    if m >= 0
        p = convert(Polynomial, q)
        return convert(PQ, p)
    else
        z = variable(q)
        zᵐ = z^(-m)
        p = convert(Polynomial, zᵐ * q)
        return rational_function(PQ, p, convert(Polynomial, zᵐ))
    end

end

function Base.convert(::Type{PQ}, p::P) where {PQ <: AbstractRationalFunction, P<:AbstractPolynomial}
    T′ = _eltype(_eltype((PQ)))
    T =  T′ == nothing ? eltype(p) : T′
    X = indeterminate(PQ, p)

    𝐩 = convert(Polynomial{T,X}, p)
    rational_function(PQ, 𝐩, one(𝐩))
end

function Base.convert(::Type{P}, pq::PQ) where {P<:AbstractPolynomial, PQ<:AbstractRationalFunction}
    p,q = pqs(pq)
    isconstant(q) || throw(ArgumentError("Can't convert rational function with non-constant denominator to a polynomial."))
    convert(P, p) / constantterm(q)
end

function Base.convert(::Type{S}, pq::PQ) where {S<:Number, T,X,P,PQ<:AbstractRationalFunction}
    !isconstant(pq) && throw(ArgumentError("Can't convert non-constant rational function to a number"))
    S(pq(0))
end

# promotion rule to promote things upwards
function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,PQ <: AbstractRationalFunction{T,X,P},
                                                           T′,X′,P′,PQ′<:AbstractRationalFunction{T′,X′,P′} }
    assert_same_variable(X,X′)
    PQ_, PQ′_ = constructorof(PQ), constructorof(PQ′)
    𝑷𝑸 = PQ_ == PQ′ ? PQ_ : RationalFunction
    𝑷 = constructorof(typeof(variable(P)+variable(P′)))
    𝑻 = promote_type(T,T′)
    𝑷𝑸{𝑻,X,𝑷{𝑻,X}}
end
Base.promote_rule(::Type{PQ}, ::Type{P}) where {PQ <: AbstractRationalFunction, P<:AbstractPolynomial} = PQ
function Base.promote_rule(::Type{PQ}, ::Type{S}) where {T,X, P<:AbstractPolynomial{T,X}, PQ <: AbstractRationalFunction{T,X,P}, S<:Number}
    R = promote_type(S,T)
    P′ = constructorof(P){R,X}
    constructorof(PQ){R,X,P′}
end

## Look like rational numbers
# The p//q constructor is reserved for the `RationalFunction` type, but these should all be defined
# for other possible types.
function Base.://(p::PQ,q::PQ′) where {PQ <: AbstractRationalFunction, PQ′ <: AbstractRationalFunction}
    p0,p1 = p
    q0,q1 = q
    rational_function(promote_type(PQ, PQ′), p0*q1, p1*q0)
end

function Base.://(p::Union{Number,AbstractPolynomial},q::PQ) where {PQ <: AbstractRationalFunction}
    q0,q1 = q
    rational_function(PQ, p*q1, q0)
end
function Base.://(p::PQ, q::Union{Number,AbstractPolynomial}) where {PQ <: AbstractRationalFunction}
    p0, p1 = p
    rational_function(PQ, p0, p1*q)
end

Base.://(p::AbstractPolynomial,q::Number) = p // (q*one(p))
Base.://(p::Number, q::AbstractPolynomial) = (p*one(q)) // q

function Base.copy(pq::PQ) where {PQ <: AbstractRationalFunction}
    p,q = pqs(pq)
    rational_function(PQ, p, q)
end

## ----

# requires these field names
Base.numerator(pq::AbstractRationalFunction) = pq.num
Base.denominator(pq::AbstractRationalFunction) = pq.den

# Treat a RationalFunction as a tuple (num=p, den=q)
Base.length(pq::AbstractRationalFunction) = 2
function Base.iterate(pq::AbstractRationalFunction, state=nothing)
    state == nothing && return (numerator(pq), 1)
    state == 1 && return (denominator(pq), 2)
    nothing
end
Base.collect(pq::AbstractRationalFunction{T,X,P}) where {T,X,P} = collect(P, pq)
Base.broadcastable(pq::AbstractRationalFunction) = Ref(pq)

Base.eltype(pq::Type{<:AbstractRationalFunction{T,X,P}}) where {T,X,P} = P
Base.eltype(pq::Type{<:AbstractRationalFunction{T,X}}) where {T,X} = Polynomial{T,X}
Base.eltype(pq::Type{<:AbstractRationalFunction{T}}) where {T} = Polynomial{T,:x}
Base.eltype(pq::Type{<:AbstractRationalFunction}) = Polynomial{Float64,:x}
_eltype(pq::Type{<:AbstractRationalFunction{T,X,P}}) where {T,X,P} = P
_eltype(pq::Type{<:AbstractRationalFunction{T,X}}) where {T,X} = Polynomial{T,X}
_eltype(pq::Type{<:AbstractRationalFunction{T}}) where {T} = Polynomial{T}
_eltype(pq::Type{<:AbstractRationalFunction})  = Polynomial

"""
    pqs(pq)

Return `(p,q)`, where `pq=p/q`, as polynomials.
"""
pqs(pq::AbstractRationalFunction) = (numerator(pq), denominator(pq))

## ----

Base.size(F::AbstractRationalFunction) = ()
function Base.isinf(pq::AbstractRationalFunction)
    p,q=pqs(pq)
    iszero(q) && !iszero(p)
end
function Base.isnan(pq::AbstractRationalFunction)
    p,q= pqs(pq)
    iszero(p) && iszero(q)
end
function Base.iszero(pq::AbstractRationalFunction)
    p,q= pqs(pq)
    iszero(p) && !iszero(q)
end
function Base.isone(pq::AbstractRationalFunction)
    p,q = pqs(pq)
    isconstant(p) && isconstant(q) && p == q
end

## -----

_indeterminate(::Type{<:AbstractRationalFunction}) = nothing
_indeterminate(::Type{PQ}) where {T,X,PQ<:AbstractRationalFunction{T,X}} = X
indeterminate(pq::PQ) where {PQ<:AbstractRationalFunction} = indeterminate(PQ)
function indeterminate(::Type{PQ}, var=:x) where {PQ<:AbstractRationalFunction}
    X′ = _indeterminate(PQ)
    X = X′ == nothing ? Symbol(var) : X′
    X
end
function indeterminate(PP::Type{P}, p::AbstractPolynomial{T,Y}) where {P <: AbstractRationalFunction, T,Y}
    indeterminate(PP, Y)
end
function isconstant(pq::AbstractRationalFunction; kwargs...)
    p,q = pqs(lowest_terms(pq, kwargs...))
    isconstant(p) && isconstant(q)
end
isconstant(::Number) = true

function constantterm(pq::AbstractRationalFunction; kwargs...)
    p,q = pqs(pq)
    isconstant(pq) && return constantterm(p)/constantterm(q)
    throw(ArgumentError("No constant term defined for a non-constant polynomial"))
end

function Base.zero(pq::R) where {R <: AbstractRationalFunction}
    p,q = pqs(pq)
    rational_function(R, zero(p), one(q))
end
function Base.one(pq::R) where {R <: AbstractRationalFunction}
    p,q = pqs(pq)
    rational_function(R, one(p), one(q))
end
function variable(pq::R) where {R <: AbstractRationalFunction}
    p,q = pqs(pq)
    rational_function(R, variable(p), one(q))
end

function variable(::Type{PQ}) where {PQ <: AbstractRationalFunction}
    x = variable(eltype(PQ))
    rational_function(PQ, x, one(x))
end

# use degree as largest degree of p,q after reduction
function degree(pq::AbstractRationalFunction)
    pq′ = lowest_terms(pq)
    maximum(degree.(pqs(pq′)))
end

# Evaluation
function eval_rationalfunction(x, pq::AbstractRationalFunction{T}) where {T}
    num, den = pqs(pq)
    dn, dd = degree(num), degree(den)
    md = min(dn, dd)
    result = num(x)/den(x)
    md < 0 && return result
    while md >= 0
        !isnan(result) && return result
        num,den = derivative(num), derivative(den)
        result = num(x)/den(x)
        md -= 1
    end

    x*NaN
end

# equality
import Base: ==
function ==(p::AbstractRationalFunction{T,X,P}, q::AbstractRationalFunction{S,Y,Q}) where {T,X,P,S,Y,Q}
    isconstant(p) && isconstant(q) && p(0) == q(0) && return true
    X == Y || return false
    p₀, p₁ = pqs(p)
    q₀, q₁ = pqs(q)
    p₀ * q₁ == q₀ * p₁ || return false
end

function ==(p::AbstractRationalFunction{T,X,P}, q::Union{AbstractPolynomial, Number}) where {T,X,P}
    ==(promote(p,q)...)
end

function ==(p::Union{AbstractPolynomial, Number}, q::AbstractRationalFunction{T,X,P}) where {T,X,P}
   ==(promote(p,q)...)
end

function Base.isapprox(pq₁::PQ₁, pq₂::PQ₂,
                  rtol::Real = sqrt(eps(float(real(promote_type(T,S))))),
                  atol::Real = zero(float(real(promote_type(T,S))))) where {T,X,P,PQ₁<:AbstractRationalFunction{T,X,P},
                                                                     S,Y,Q,PQ₂<:AbstractRationalFunction{S,Y,Q}}

    p₁,q₁ = pqs(pq₁)
    p₂,q₂ = pqs(pq₂)

    isapprox(p₁*q₂, q₁*p₂; rtol=rtol, atol=atol)
end

# Arithmetic
function Base.:-(pq::PQ) where {PQ <: AbstractRationalFunction}
    p, q = copy(pq)
    rational_function(PQ, -p, q)
end

Base.:+(p::Number, q::AbstractRationalFunction) = q + p
Base.:+(p::AbstractRationalFunction,  q::Number) = p + q*one(p)
Base.:+(p::AbstractPolynomial, q::AbstractRationalFunction) = q + p
Base.:+(p::AbstractRationalFunction,  q::AbstractPolynomial) = p + (q//one(q))
Base.:+(p::AbstractRationalFunction, q::AbstractRationalFunction) = sum(promote(p,q))
# type should implement this
function Base.:+(p::R, q::R) where {T,X,P,R <: AbstractRationalFunction{T,X,P}}
    p0,p1 = pqs(p)
    q0,q1 = pqs(q)
    rational_function(R, p0*q1 + p1*q0, p1*q1)
end

Base.:-(p::Number, q::AbstractRationalFunction) = -q +  p
Base.:-(p::AbstractRationalFunction,  q::Number) = p - q*one(p)
Base.:-(p::AbstractPolynomial, q::AbstractRationalFunction) = -q + p
Base.:-(p::PQ,  q::AbstractPolynomial) where {PQ <: AbstractRationalFunction} = p - rational_function(PQ,q, one(q))
function Base.:-(p::AbstractRationalFunction, q::AbstractRationalFunction)
    p′, q′ = promote(p,q)
    p′ - q′
end
# type should implement this
function Base.:-(p::R, q::R) where {T,X,P,R <: AbstractRationalFunction{T,X,P}}
    p0,p1 = pqs(p)
    q0,q1 = pqs(q)
    rational_function(R, p0*q1 - p1*q0, p1*q1)
end

function Base.:*(p::Number, q::R) where {T, X, R <: AbstractRationalFunction{T,X}}
    q0,q1 = pqs(q)
    rational_function(R, (p*q0), q1)
end
function Base.:*(p::R,  q::Number) where {R <:AbstractRationalFunction}
    p0,p1 = pqs(p)
    rational_function(R, p0*q, p1)
end
function Base.:*(p::AbstractPolynomial, q::R) where {R <: AbstractRationalFunction}
    rational_function(R, p, one(p)) * q
end
Base.:*(p::R,  q::AbstractPolynomial) where {R <: AbstractRationalFunction} = p * rational_function(R,q, one(q))
# type should implement this
Base.:*(p::AbstractRationalFunction, q::AbstractRationalFunction) = prod(promote(p,q))
function Base.:*(p::R, q::R) where {T,X,P,R <: AbstractRationalFunction{T,X,P}}
    p0,p1 = pqs(p)
    q0,q1 = pqs(q)
    rational_function(R, p0*q0, p1*q1)
end

function Base.:/(p::Number, q::R) where {R <: AbstractRationalFunction}
    q0,q1 = pqs(q)
    rational_function(R,p*q1, q0)
end
function Base.:/(p::R,  q::Number) where {R <: AbstractRationalFunction}
    p0,p1 = pqs(p)
    rational_function(R, p0, (p1*q))
end
Base.:/(p::AbstractPolynomial, q::PQ) where {PQ <: AbstractRationalFunction} = rational_function(PQ, p,one(p)) / q
function Base.:/(p::PQ,  q::AbstractPolynomial) where {PQ <: AbstractRationalFunction}
    p0,p1 = pqs(p)
    rational_function(PQ,p0, p1*q)
end
function Base.:/(p::AbstractRationalFunction, q::AbstractRationalFunction)
    p′,q′ = promote(p,q)
    p′ / q′
end
# type should implement this
function Base.:/(p::PQ, q::PQ) where {T,X,P,PQ <: AbstractRationalFunction{T,X,P}}
    p0,p1 = pqs(p)
    q0,q1 = pqs(q)
    rational_function(PQ, p0*q1, p1*q0)
end

function Base.:^(pq::P, n::Int) where {P <: AbstractRationalFunction}
    p,q = pqs(pq)
    rational_function(P, p^n, q^n)
end

function Base.inv(pq::P) where {P <: AbstractRationalFunction}
    p,q = pqs(pq)
    rational_function(P, q, p)
end

# conj, transpose... TODO

## derivative and integrals
function derivative(pq::P, n::Int=1) where {P <: AbstractRationalFunction}
    n <= 0 && return pq
    while n >= 1
        p,q = pqs(pq)
        pq = rational_function(P, (derivative(p)*q - p * derivative(q)), q^2)
        n -= 1
    end
    pq
end

function integrate(pq::P) where {P <: AbstractRationalFunction}
    p,q = pqs(pq)
    isconstant(q) && return rational_function(P, integrate(p), q)
    # XXX could work here, e.g.:
    # d,r = partial_fraction
    # ∫d + Σ∫r for each residue (issue with logs)
    throw(ArgumentError("Can only integrate rational functions with constant denominators"))
end

## ----

"""
    divrem(pq::AbstractRationalFunction; method=:numerical, kargs...)

Return `d,r` with `p/q = d + r/q` where `degree(numerator(r)) < degree(denominator(q))`, `d` a Polynomial, `r` a `AbstractRationalFunction`.

* `method`: passed to `gcd`
* `kwargs...`: passed to `gcd`
"""
function Base.divrem(pq::PQ; method=:numerical, kwargs...) where {PQ <: AbstractRationalFunction}
    p,q = pqs(pq)
    degree(p) < degree(q) && return (zero(p), pq)

    d,r = divrem(p,q)
    (d, rational_function(PQ, r, q))

end

"""
    lowest_terms(pq::AbstractRationalFunction, method=:numerical)

Find GCD of `(p,q)`, `u`, and return `(p÷u)//(q÷u)`. Commonly referred to as lowest terms.

* `method`: passed to `gcd(p,q)`
* `kwargs`: passed to `gcd(p,q)`

By default, `AbstractRationalFunction` types do not cancel common factors. This method will numerically cancel common factors, returning the normal form, canonicalized here by `q[end]=1`. The result and original may be considered equivalent as rational expressions, but different when seen as functions of the indeterminate.

"""
function lowest_terms(pq::PQ; method=:numerical, kwargs...) where {T,X,
                                                                   P<:StandardBasisPolynomial{T,X},
                                                                   PQ<:AbstractRationalFunction{T,X,P}}
    p,q = pqs(pq)
    u,v,w = uvw(p,q; method=method, kwargs...)
    rational_function(PQ, v/w[end], w/w[end])
end

## ---- zeros, poles, ...
"""
    poles(pq::AbstractRationalFunction; method=:numerical, kwargs...)

For a rational function `p/q`, first reduces to normal form, then finds the roots and multiplicities of the resulting denominator.

"""
function poles(pq::AbstractRationalFunction; method=:numerical,  kwargs...)
    pq′ = lowest_terms(pq; method=method, kwargs...)
    den = denominator(pq′)
    mr = Multroot.multroot(den)
    (zs=mr.values, multiplicities = mr.multiplicities)
end

"""
    roots(pq::AbstractRationalFunction; kwargs...)

Return the `zeros` of the rational function (after cancelling commong factors, the `zeros` are the roots of the numerator.

"""
function roots(pq::AbstractRationalFunction; method=:numerical,  kwargs...)
    pq′ = lowest_terms(pq; method=method, kwargs...)
    den = numerator(pq′)
    mr = Multroot.multroot(den)
    (zs=mr.values, multiplicities = mr.multiplicities)
end

"""
    residues(pq::AbstractRationalFunction; method=:numerical,  kwargs...)

If `p/q =d + r/q`, returns `d` and the residues of a rational fraction `r/q`.

First expresses `p/q =d + r/q` with `r` of lower degree than `q` through `divrem`.
Then finds the poles of `r/q`.
For a pole, `λj` of multiplicity `k` there are `k` residues,
`rⱼ[k]/(z-λⱼ)^k`, `rⱼ[k-1]/(z-λⱼ)^(k-1)`, `rⱼ[k-2]/(z-λⱼ)^(k-2)`, …, `rⱼ[1]/(z-λⱼ)`.
The residues are found using this formula:
`1/j! * dʲ/dsʲ (F(s)(s - λⱼ)^k` evaluated at `λⱼ` ([5-28](https://stanford.edu/~boyd/ee102/rational.pdf)).

## Example

(From page 5-33 of above pdf)

```jldoctest
julia> s = variable(Polynomial, :s)
Polynomial(1.0*s)

julia> pq = (-s^2 + s + 1) // ((s-1) * (s+1)^2)
(1.0 + 1.0*s - 1.0*s^2) // (-1.0 - 1.0*s + 1.0*s^2 + 1.0*s^3)

julia> d,r = residues(pq);

julia> d
Polynomial(0.0)

julia> r
Dict{Float64, Vector{Float64}} with 2 entries:
  -1.0 => [-1.25, 0.5]
  1.0  => [0.25]

julia> iszero(d)
true

julia> z = variable(pq)
(1.0*s) // (1.0)

julia> for (λ, rs) ∈ r # reconstruct p/q from output of `residues`
           for (i,rᵢ) ∈ enumerate(rs)
               d += rᵢ/(z-λ)^i
           end
       end

julia> p′, q′ = lowest_terms(d);

julia> q′ ≈ (s-1) * (s+1)^2 # works, as q is monic
true

julia> p′ ≈ (-s^2 + s + 1)
true
```

!!! Note:
    There are several areas where numerical issues can arise. The `divrem`, the identification of multiple roots (`multroot`), the evaluation of the derivatives, ...

"""
function residues(pq::AbstractRationalFunction; method=:numerical,  kwargs...)

    d,r′ = divrem(pq)
    r = lowest_terms(r′; method=method, kwargs...)
    b,a = pqs(r)
    a′ = derivative(a)

    residues = Any[]
    mr = Multroot.multroot(a)

    for (λₖ, mₖ) ∈ zip(mr.values, mr.multiplicities)

        if mₖ == 1
            push!(residues, λₖ => [b(λₖ)/a′(λₖ)])
        else
            # returns rₖ₁, rₖ₂,...rₖₘ where rₖ,i/(s-λₖ)ⁱ is part of the decomposition
            s = variable(a)
            F = lowest_terms(r*(s-λₖ)^mₖ)
            rs = [F(λₖ)]
            j! = 1
            for j ∈ 1:mₖ-1
                dF = lowest_terms(derivative(F))
                pushfirst!(rs, 1/j! * dF(λₖ))
                j! *= (j+1)
            end
            push!(residues, λₖ => rs)
        end
    end
    d, Dict(residues...)
end

"""
    partial_fraction(pq::AbstractRationalFunction; method=:numerical, kwargs...)

For a rational function `p/q = d + r/q`, with `degree(r) < degree(q)` returns `d` and
the terms that comprise `r/q`: For each pole with multiplicity returns
`rⱼ[k]/(z-λⱼ)^k`, `rⱼ[k-1]/(z-λⱼ)^(k-1)`, `rⱼ[k-2]/(z-λⱼ)^(k-2)`, …, `rⱼ[1]/(z-λⱼ)`.

Should be if `p/q` is in normal form and `d,r=partial_fraction(p//q)` that
`d + sum(r) - p//q ≈ 0`

"""
function partial_fraction(pq::AbstractRationalFunction; method=:numerical, kwargs...)
    d,r = residues(pq; method=method, kwargs...)
    s = variable(pq)
    d, partial_fraction(Val(:residue), r, s)
end

function partial_fraction(::Val{:residue}, r, s::PQ) where {PQ}
    terms = []
    for (λₖ,rsₖ) ∈ r
        for  (rⱼ,j) ∈ zip(rsₖ, length(rsₖ):-1:1)
            push!(terms, rⱼ/(s-λₖ)^j)
        end
    end
    terms

end


text:
import PotentialUQ
using LinearAlgebra: I
using Random
using Turing
Random.seed!(1234);

n = 100
β = zeros(10)
A = 1e-1*randn(n, 10)
A[1:10, 1:10] += I(10)
b = randn(n)

Q = Matrix(1e-3*I(n))

snap = PotentialUQ.SNAP(A, b, β)
dsnap = PotentialUQ.SNAPDistribution(snap, Q)

samples = PotentialUQ.sample(dsnap; num_samples = 1000)

println("MAP Solution, Sampling ", dsnap.x.β)

# summaries = describe(samples)
# show(Base.stdout, "text/plain", summaries[1])
# show(Base.stdout, "text/plain", summaries[2])

ext: jl


text:
using DifferentialEquations
using Distributions
using Plots
using DelimitedFiles

#-----------------------------------
#specify the signaling topology
#-----------------------------------

topology = [0 1 1
            0 0 1
            0 0 0]

nCell = size(topology, 1)

#-----------------------------------
#specify the number of iterations
#-----------------------------------

nits = 10

#-----------------------------------
#specify the value of A0
#-----------------------------------

A0 = 0.7

#-----------------------------------
#specify the value of lambda
#-----------------------------------

lambda = 1.0
kappa = 1.0

#-------------------------------------------------
#specify the wait time distribution and pulse time
#-------------------------------------------------

mean = 50.0

pulse = 5.0

#-----------------------------------
# internal ODE system
#-----------------------------------

function GataPu1(du,u,p,t)
    a = p[1:5]
    b = p[6:8]
    g = p[9:17]
    L = p[18:20] #L = [A,B,C]
    du[1] = -b[1]*u[1] + (a[1]*L[1] + a[2]*u[1])/(1 + g[1]*L[1] + g[2]*u[1] + g[3]*u[1]*u[2])
    du[2] = -b[2]*u[2] + (a[3]*L[2] + a[4]*u[2])/(1 + g[4]*L[2] + g[5]*u[2] + g[6]*u[1]*u[2] + g[7]*u[1]*u[3])
    du[3] = -b[3]*u[3] + (a[5]*u[1])/(1 + g[8]*u[1] + g[9]*L[3])
end

#-----------------------------------
# parameter vector p
#-----------------------------------

p = zeros(0)
a = [1.0,0.25,1.0,0.25,0.01]
b = [0.01,0.01,0.01]
g = [1.0,0.25,1.0,1.0,0.25,1.0,0.13,0.01,10]
L = [1.0,0.5,0] #L = [A,B,C] = [A,0.5,0]
append!(p,a)
append!(p,b)
append!(p,g)
append!(p,L)

#-----------------------------------
# wait time distribution
#-----------------------------------

function rand_exponential(mean)
    if mean <= 0.0
        error("mean must be positive")
    end
    -mean*log(rand())
end

Nml = Normal(50.0, 10.0)

function wait_time(mean)
    round(rand_exponential(mean))
end

#-----------------------------------
#Simulation Code
#-----------------------------------

CellValues = [[20.0, 1.0, 20.0, A0] for i=1:nCell, n=1:1]
CellIts = zeros(nCell)
CellWaitTimes = [wait_time(mean) for i=1:nCell, n=1:1]
CellState = ["wait" for i=1:nCell, n=1:1]
PulseCoeffs = zeros(Float64, nCell, nCell)

### Define empty vectors for G,P,X concentrations and times for each cell
for x in [Symbol("Cell"*string(i)) for i=1:nCell]
    @eval $x = Any[]
end
for x in [Symbol("Time"*string(i)) for i=1:nCell]
    @eval $x = Any[]
end

total_time = 0.0

while min(CellIts...) < 1000 + 40*nCell
    t = min(CellWaitTimes...)

    for i = 1:nCell
        u0 = [CellValues[i][1], CellValues[i][2], CellValues[i][3]]
        p[18] = CellValues[i][4]
        tspan = (0.0, t)
        prob = ODEProblem(GataPu1,u0,tspan,p)
        sol = solve(prob)
        append!(eval(Symbol("Cell"*string(i))), sol.u)
        append!(eval(Symbol("Time"*string(i))), sol.t .+ total_time)
        CellValues[i][1] = sol[length(sol)][1]
        CellValues[i][2] = sol[length(sol)][2]
        CellValues[i][3] = sol[length(sol)][3]
    end

    global CellWaitTimes = CellWaitTimes .- t
    global total_time = total_time + t

    for i = 1:nCell
        if CellWaitTimes[i] == 0
            if CellState[i] == "wait"
                for m = 1:nCell
                    if topology[i,m] == 1
                        PulseCoeffs[i,m] = (CellValues[i][1]/(lambda + CellValues[i][2]))
                        CellValues[m][4] = CellValues[m][4] * PulseCoeffs[i,m]
                    elseif topology[i,m] == -1
                        PulseCoeffs[i,m] = (CellValues[i][2]/(kappa + CellValues[i][1]))
                        CellValues[m][4] = CellValues[m][4] * PulseCoeffs[i,m]
                    end
                end
                CellWaitTimes[i] = pulse
                CellState[i] = "pulse"
                CellIts[i] = CellIts[1] + 1
            else
                for m = 1:nCell
                    if topology[i,m] == 1
                        CellValues[m][4] = CellValues[m][4] / PulseCoeffs[i,m]
                    elseif topology[i,m] == -1
                        CellValues[m][4] = CellValues[m][4] / PulseCoeffs[i,m]
                    end
                end

                CellWaitTimes[i] = wait_time(mean)
                CellState[i] = "wait"
            end
        end
    end
end

for i = 1:nCell
    writedlm("../results/Trajectory_Cell_$i.txt", eval(Symbol("Cell"*string(i))), ", ")
end

for i = 1:nCell
    writedlm("../results/Trajectory_Time_$i.txt", eval(Symbol("Time"*string(i))), ", ")
end

max_stars_repo_name: alavendelm/Cell-Cell-Communication

max_issues_repo_name: alavendelm/Cell-Cell-Communication

max_forks_repo_name: alavendelm/Cell-Cell-Communication


text:
colorsets=[[(1,0,0),(0,1,0),(0,0,1),(1/2,0,1/2),(1/4,1/2,1/4),(1,1,1)],[(1,0,0),(0,1,0),(0,0,1),(1/2,1/2,0),(1/2,0,1/2),(0,1/2,1/2),(0.15,0.7,0.15)]]
setcolorset=(game,set::Int=2)->game.colors=colorsets[set] #modify game.color to set a temporary custom color

units=Dict("standard"=>Dict(:name=>"standard", :ir=>2, :pl=>[2])) #change :pl to [1,2,3] to place anywhere

units["tunnel"]=Dict(:name=>"tunnel", :ir=>0, :pl=>[1]) #these by default cannot be placed under/over surrounded territory without permission, such permissions are not hard coded but may rather result in a rejected/modified merge request
units["bridge"]=Dict(:name=>"bridge", :ir=>0, :pl=>[3])
units["queen"]=Dict(:name=>"queen", :canspawn=>true, :graphic=>[(0.5,0.5),(-0.5,0.5),(0,-0.5)])
units["vein"]=Dict(:name=>"vein", :ir=>0, :pl=>[2])
#units["universal"]=Dict(:name=>"universal", :canspawn=>true, :ir=>3, :pl=>[1,2,3],:graphic=>[(0.5,-0.5),(-0.5,-0.5),(0,0.5)]) #this one has some weird scoring issues, sometimes adds score then subtracts it the next move. Anyways not really meant to be generally used since it overrides the need to connect (like playing with only (super)queens), more like a GM unit
#units[:storage] :canstore=>true, :storemax=> #for storing (duh!)

#white unit is a spawn for all colors whereas any unit with color (1,1,1) "is" every color, so a white spawn is a spawn for all and "white" as a unit spec can be removed. White superpowers are being reduced, "units" facilitating multiple color connections can be built with tunnel and bridge units. White now has metaphysical superpowers, perfect storage
#units["white"]=Dict(:harvest=>white_harvest, :name=>"white", :inclusive=>true, :passover=>true, :baselife=>18, :canspawn=>true)
#build white unitsets, like molecules with colored branches
function standard_costfun(game,unit) #costs have been removed, our hive will dispatch all units proposed that create agreeable formations
	distance=nearestwhite(game,unit.loc,unit.loc[3]==2)
	cost=distance*unitparams[2]
	cost*=sum(unitparams[1]) #only if sum>1? How do we create such units? Are smaller units good? Do they take less space?
	return cost
end
function standard_live!(game,unit,lifemap)
	spreadlife!(game,unit,lifemap)
end
function standard_harvest(game,unit,ledger)
	unitharvest(game,unit,ledger)
end
standard_harvest! =standard_harvest #for savefile backwardscompat


text:
using SparseArrays

function edge_values(space, m)
    rs = refspace(space)
    supp = geometry(space)
    edgs = skeleton(supp,1)
    Conn = copy(transpose(connectivity(edgs, supp, identity)))
    Vals = zeros(scalartype(space), length(supp), length(edgs))
    for (i,sh) in enumerate(space.fns[m])
        tet_id = sh.cellid
        tet = chart(supp, cells(supp)[tet_id])
        for k in nzrange(Conn, tet_id)
            edg_id = rowvals(Conn)[k]
            # loc_id = abs(nonzeros(Conn)[k])
            # edge = CompScienceMeshes.edges(tet)[loc_id]
            edge = chart(edgs, cells(edgs)[edg_id])
            ctr = center(edge)
            tgt = tangents(ctr,1)
            ctr1 = neighborhood(tet, carttobary(tet, cartesian(ctr)))
            vals = rs(ctr1)
            # @show vals
            val = vals[sh.refid].value
            Vals[tet_id, edg_id] += dot(sh.coeff * val, tgt)
        end
    end
    return Vals
end

function check_edge_values(EV)
    for i in axes(EV,2)
        ev = [x for x in EV[:,i] if abs(x) > 1e-8]
        uev = unique(x -> round(x, digits=5), ev)
        if !(length(uev) <= 1)
            @show i
            @show uev
            error("stop")
        end
    end
end


text:
module Sequoia

using Compat, NullableArrays

export  # modules
        CSV, SQLite3,

        # types
        C, F, Cursor, Table, NullableTable,

        # consts
        sqlite3_lib,

        # methods
        ncols,
        reset!,
        step!

include("../deps/deps.jl")
include("cursor.jl")
# in
include("in/sqlite3.jl"); # using SQLite3
include("in/csv.jl")
# out
include("out/table.jl")
include("out/nullabletable.jl")
# include("out/matrix.jl")

end # module


text:
using Colors
include("helpers/listing.jl")

isinstalled(pkg) = try Pkg.installed(pkg) != nothing catch e false end

function main(window)
    push!(window.assets, "widgets")
    push!(window.assets, "codemirror")
    push!(window.assets, "tex")

    vbox(
        title(3, "Typography"),
        vskip(1em),
        md"Escher provides primitives that directly map to CSS font styling properties, as well as higher-level functions which form a standard typographic scale you can use to give your documents a consistent aesthetic.",
        vskip(1em),
        h1("High-level Functions"),
        md"""These use styles from $(Pkg.dir("Escher"))/assets/font.css by default""", # It would be better to have the served asset directory instead

        vskip(1em),
        h2("Block Quotes"),
        md"Block quotes is a style used when quoting large sections of text. They are also available within markdown by prepending a section with `>` ",
        vskip(1em),
        listing("""vbox(
          blockquote(md" `map` is a function that that applies a given **function** to each element of a **list**, returning a **list of results**."),
          md"e.g. `map(x->x*x, 1:4)` evaluates to \$(string(map(x->x*x, 1:4)))",
        )
        """),

        vskip(1em),
        h2("Titles and Headings"),
        "Useful things you might use if writing an article like this one",
        vskip(1em),
        listing("""vbox(
  title(4,"Title 4"),
  map(n->title(n,"Title \$n"), 3:-1:1 )...,
  heading(1,"Heading 1"),
  map(n->heading(n,"Heading \$n"), 2:4 )...,
  "Body Text",
  h1( "h1(x) is short for header(1,x)" ),
  h2( "h2" ),
  h3( "h3" ),
  h4( "h4" ),
  md\"\"\"# Heading 1 in Markdown
  ## Heading 2 in Markdown
  ### Heading 3 in Markdown
  #### Heading 4 in Markdown\"\"\",
)"""),

        h2("Code"),
        md"""To show code with syntax highlighting, you can use the `codemirror` function. Codemirror requires the `"codemirror"` asset.

        Use `push!(window.assets, "codemirror")` to do this.""",
        vskip(1em),
        listing("""
            push!(window.assets, "codemirror")

            codemirror(\"\"\"
                function foo()
                42
                end
                \"\"\"
            )
            """,
            codemirror("""
                function foo()
                    42
                end
                """
            )
        ),
        md"""Note that you also need the `"codemirror"` asset to be loaded even if you are writing code inside `md`.""",

        vskip(1em),
        h2("LaTeX"),
        md"LaTeX strings can be rendered with the `tex` function.",
        vskip(1em),
        listing("""
        tex("(a+b)^2=a^2+b^2+2ab")"""),

        vskip(1em),
        h1("Low-level Functions"),
        "Sometimes you want more explicit control over your document",

        vskip(1em),
        h2("Weights and Types"),
        md"Many classes of fonts are available, often with many weights. By default, Escher uses the [*Source Sans Pro*](http://www.google.com/fonts/specimen/Source+Sans+Pro) (sans-serif) and [*Source Code Pro*](http://www.google.com/fonts/specimen/Source+Code+Pro) (monospaced) font families for great-looking and legible type.",
md"`fontweight` can take integer multiples of 100 from 100 to 900, or values such as `bold`, `bolder`, and `lighter` ",
md"`fonttype` arguments include `serif`, `sansserif`, `slabserif`, and `monospace`.",
        vskip(1em),
        listing("""hbox(intersperse(hskip(1em),[
  vbox(map(x->fontweight(x,"Sans \$x."), 100:200:700 )...),
  vbox(map(x->fontweight(x,"Mono \$x.") |> fonttype(monospace), 300:200:700 )...),
  vbox(map(x->fontweight(x,"Serif \$x.") |> fonttype(serif), 500:200:700 )...),
  fontweight(bold, "Bold flavored text")
]))"""),

        vskip(1em),
        h2("Styling and Alignment"),
        "For those time where your squiggly lines are leaning, or if you want to write on the other side",
        md"Valid alignments are `raggedright`, `raggedleft`, `justifytext`, and `centertext`",
        md"note that `fontcase` only changes the style, and not the content of the text. Valid arguments are `ucase` and `lcase`",
        vskip(1em),
        listing("""vbox(
  fontstyle(italic,"italics is a CSS class used for styling text"),
  emph("emphasis is an HTML tag for adding semantic emphasis to a word or phrase"),
  textalign(centertext,fontcase(ucase, "Loudness equals power." )),
  fontcase(lcase, "Speak softly, and carry a BIG stick." ) |> textalign(raggedleft),
  lineheight(80px,"This text has a high lineheight"),
)"""),

        vskip(1em),
        h2("Colors and Families"),
        md"You can also change the color of text with *fontcolor* and [Colors.jl](https://github.com/JuliaGraphics/Colors.jl). You may change the family with *fontfamily*",
        vskip(1em),
        listing("""using Colors
hbox(
  map(h->fontcolor(LCHab(50,200,h),"nya "),20:20:320)...,
  fontcolor(colorant"magenta","ncat :3")
) |> fontfamily("Ubuntu")"""),

        vskip(1em),
        h2("Size"),
        md"```fontsize``` accepts many units of length such as `pt` and `px`, as well as keywords `((x)x)large|small` and `medium`",
        listing("""vbox(
  Escher.fontsize(xxlarge,"Big Text"),
  Escher.fontsize(36pt,"Bigger Text"),
  Escher.fontsize(48px,"Bigger Text"),
)"""),
        vskip(1em),

        title(3, "Embellishment"),
        vskip(1em),

md"A border may be styled by many functions. `borderstyle` may accept `solid`, `dashed`, `dotted`, and `noborder`. `fillcolor` and `bordercolor` can accept any color defined in colors.jl, and `borderwidth` may accept anything of type Length such as `em` or `px`.",
        listing("""greybox = container(12em,2.5em) |> fillcolor(colorant"grey") |> borderwidth(0.4em)

vbox(
 greybox,
 greybox |> bordercolor(colorant"orange") |> borderstyle(solid),
 greybox |> borderstyle(solid),
 greybox |> borderstyle(dotted),
 greybox |> borderstyle(dashed),
)"""),

        vskip(1em),
md"You may specify a border all at once with `border`. It takes a list of sides (optionally), a style, a width, and a color.",
        listing("""image(
    "https://upload.wikimedia.org/wikipedia/commons/thumb/5/5f/Utah_teapot_simple_2.png/220px-Utah_teapot_simple_2.png",
    alt="Teapot",
) |> border([right,top,left,bottom],dashed,2px,colorant"#2D2E30")"""),

        vskip(1em),
        md"`hline` and `vline` are very thin, bordered elements that are inteded for use as a separator. It is quite useful within a flow: `vline` within a `hbox` and `hline` in a `vbox`. As with any other bordered element, they can be styled accordingly",
        listing("""vbox(
          paper(4,hbox(
            [
              flex("Column A"),
              vline(),
              flex("Column B")
            ]
          )),
          vskip(4em),
          paper(4,vbox(
            [
              "Row 1",
              hline() |> borderwidth(3px),
              "Row 2"
            ]
          )),
        )"""),

    ) |> pad(2em)
end


text:
export NumericParser, DateParser, DateTimeParser
using TextParse
import Dates
import Dates: DateFormat

@deprecate Numeric(x...) NumericParser(x...)

DateParser(format::AbstractString...; locale="english")     = DateParser(Dates.DateFormat.(format, locale)...)
DateTimeParser(format::AbstractString...; locale="english") = DateTimeParser(Dates.DateFormat.(format, locale)...)
DateParser(format::DateFormat...)     = Either(parser.(TextParse.DateTimeToken.(Dates.Date,format))...; simplify=true)
DateTimeParser(format::DateFormat...) = Either(parser.(TextParse.DateTimeToken.(Dates.DateTime,format))...; simplify=true)

"""
    DateParser(format::DateFormat...)
    DateTimeParser(format::DateFormat...)

Create a parser matching either one format
using `TextParse.DateTimeToken` for `Dates.Date` and `Dates.DateTime` respectively.

    DateParser(format::AbstractString...; locale="english")
    DateTimeParser(format::AbstractString...; locale="english")

Convenience functions for above using `Dates.DateFormat.(format, locale)`.
"""
DateParser, DateTimeParser

import TextParse: tryparsenext

result_type(x::AbstractToken) = result_type(typeof(x))
result_type(::Type{<:AbstractToken{T}}) where T = T

struct AbstractTokenParser{P<:AbstractToken,T} <: LeafParser{NCodeunitsState{T},T}
    parser::P
    function AbstractTokenParser(p::AbstractToken)
        new{typeof(p), result_type(p)}(p)
    end
    function AbstractTokenParser{T}(a...; kw...) where T
        p = T(a...; kw...)
        new{typeof(p), result_type(p)}(p)
    end
end
"""
    NumericParser(x...) = parser(TextParse.Numeric(x...))
"""
NumericParser{T} = AbstractTokenParser{T} where {T <: TextParse.Numeric}
NumericParser(x...) = parser(TextParse.Numeric(x...))

parser(x::AbstractToken) = AbstractTokenParser(x)

regex_string(::TextParse.Numeric{<:Integer}) = "-?[[:digit:]]+"
function _printnode(io::IO, x::AbstractTokenParser)
    print_constructor(io, x)
end

print_constructor(io::IO, x::AbstractTokenParser) =
    print(io, x.parser)

print_constructor(io::IO, x::AbstractTokenParser{<:TextParse.DateTimeToken}) =
    print(io, x.parser.format)

iterate_state(parser::AbstractTokenParser, sequence, till, before_i, next_i, state::NCodeunitsState) =
    nothing

function iterate_state(parser::AbstractTokenParser, sequence, till, before_i, next_i, state::Nothing, opts=TextParse.default_opts)
    r,next_i_ = tryparsenext(parser.parser, sequence, next_i, till,opts)
    if isnull(r)
        nothing
    else
        NCodeunitsState(next_i,next_i_,get(r))
    end
end

"""
    TextParse.tryparsenext(x::CombinedParser,str,i,till,opts=TextParse.default_opts)

`TextParse.jl` integrates with `CombinedParsers.jl` both ways.
> `tryparsenext` returns a tuple `(result, nextpos)` where `result` is of type `Nullable{T}`, `Nullable{T}()` if parsing failed, non-null containing the parsed value if it succeeded.
> If parsing succeeded, `nextpos` is the position the next token, if any, starts at. If parsing failed, `nextpos` is the position at which the
parsing failed.

```jldoctest
julia> using TextParse

julia> p = ("Number:" * Repeat(' ') * TextParse.Numeric(Int))[3]
🗄 Sequence[3]
├─ Number\\:
├─ \\ *  |> Repeat
└─ <Int64>
::Int64

julia> parse(p, "Number:    42")
42

julia> TextParse.tryparsenext(p, "Number:    42")
(Nullable{Int64}(42), 14)
```

"""
function TextParse.tryparsenext(x::CombinedParser,str,i,till,opts=TextParse.default_opts)
    s = iterate_state(x,str,till,i,nothing)
    if s === nothing
        Nullable{result_type(x)}(),i
    else
        Nullable(get(x,str,till,tuple_pos(s),i,tuple_state(s))),tuple_pos(s)
    end
end

"""
    regex_prefix(x)

Prefix printed in parser tree node.
"""
regex_prefix(x::AbstractToken) = ""
"""
    regex_suffix(x)

Suffix printed in parser tree node.
"""
regex_suffix(x::AbstractToken) = ""
"""
    regex_inner(x::AbstractToken)

Regex representation of `x`.
See [`regex_string`](@ref)
"""
regex_inner(x::AbstractToken) = "$(typeof(x))"
regex_inner(::TextParse.Numeric{T}) where T = "$(T)"


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

using DifferentialEquations

function main()

function micro_1(du, u, parms, time)
    # PARAMETER DEFS
    # β transmition rate
    # r net population growth rate
    # μ hosts' natural mortality rate
    # Κ population size
    # α disease induced mortality rate

    β, r, μ, K, α = parms
    dS = r*(1-S/K)*S - β*S*I
    dI = β*S*I-(μ+α)*I
    du = [dS,dI]
end

# +
# PARAMETER DEFS
# w and m are used to define the other parameters allometrically

w = 1;
m = 10;
β = 0.0247*m*w^0.44;
r = 0.6*w^-0.27;
μ = 0.4*w^-0.26;
K = 16.2*w^-0.7;
α = (m-1)*μ;
# -

parms = [β,r,μ,K,α];
init = [K,1.];
tspan = (0.0,10.0);

sir_prob = ODEProblem(micro_1,init,tspan,parms)

sir_sol = solve(sir_prob);

end


text:
realtype(x::DataType)           = x
realtype{T}(::Type{Complex{T}}) = T

complextype(T::DataType)           = Complex{T}
complextype{T}(::Type{Complex{T}}) = Complex{T}

limitprecision(x)                     = x
limitprecision(x::Float64)            = Float32(x)
limitprecision(x::Complex128)         = Complex64(x)
limitprecision(x::Vector{Float64})    = convert(Vector{Float32}, x)
limitprecision(x::Vector{Complex128}) = convert(Vector{Complex64}, x)
limitprecision{T}(::Type{T})          = T
limitprecision(::Type{Float64})       = Float32
limitprecision(::Type{Complex128})    = Complex64

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
# Intersect
# ---------------------------------------------------------------------------------------------
#
# Intersection of a line segment and polygon in 2D
function intersect(l::LineSegment2D{T}, poly::Polygon{N, 2, T}) where {N,T}
    # Create the line segments that make up the polygon and intersect each one
    # until 2 unique points have been found
    p₁ = nan(Point2D{T})
    npoints = 0x0000
    for i ∈ 1:N-1
        hit, point = l ∩ LineSegment2D(poly[i], poly[i + 1])
        if hit
            if npoints === 0x0000
                npoints = 0x0001
                p₁ = point
            elseif !(p₁ ≈ point)
                return true, SVector(p₁, point)
            end
        end
    end
    hit, point = l ∩ LineSegment2D(poly[N], poly[1])
    if hit
        if npoints === 0x0000
            npoints = 0x0001
            p₁ = point
        elseif !(p₁ ≈ point)
            return true, SVector(p₁, point)
        end
    end
    return false, SVector(p₁, point)
end

# Möller, T., & Trumbore, B. (1997). Fast, minimum storage ray-triangle intersection.
function intersect(l::LineSegment3D{T}, tri::Triangle3D{T}) where {T}
    p = Point3D{T}(0, 0, 0)
    𝗲₁ = tri[2] - tri[1]
    𝗲₂ = tri[3] - tri[1]
    𝗱 = l.𝘂
    𝗽 = 𝗱 × 𝗲₂
    det = 𝗽 ⋅ 𝗲₁
    (det > -1e-8 && det < 1e-8) && return (false, p)
    inv_det = 1/det
    𝘁 = l.𝘅₁ - tri[1]
    u = (𝘁 ⋅ 𝗽)*inv_det
    (u < 0 || u > 1) && return (false, p)
    𝗾 = 𝘁 × 𝗲₁
    v = (𝗾 ⋅ 𝗱)*inv_det
    (v < 0 || u + v > 1) && return (false, p)
    t = (𝗾 ⋅ 𝗲₂)*inv_det
    (t < 0 || t > 1) && return (false, p)
    return (true, l(t))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
__precompile__()

module BasisMatrices

#=
TODO: still need to write fund, minterp
TODO: also need splidop, lindop
TODO: funeval fails for scalar input and does weird thing for 1-element
      vector input

TODO: rethink the `order::Matrix` thing. Maybe a convention we can adopt is
the following:

- If calling evalbase on a `*Params` type, we will say that order must be of
type `Integer` or `T<:Integer, Vector{T}`. The implementation for the `Vector`
case will just loop over each element and then call the `Integer` method,
which will have all the code
- If calling `evalbase` on a `Basis{1}` type, the semantics will be the same
as for the `*Params` case
- If calling `evalbase` on a `Basis{N>1}` type, then `order` can either be a
`T<:Integer Vector{T}` or `T<:Integer Matrix{T}`. If it is a vector, we call
the `Int` evalbase for the `*Params` along each dimension, then combine as
necessary. If it is a `Matrix` then we say that each column is the orders that
should be computed for each dimension and we farm it out that way.

So method signatures would be

```
evalbase(p::ChebParam, order::Int) => Matrix{Float64}
evalbase(p::SplineParam, order::Int) => SparseMatrixCSC{Float64}
evalbase(p::LinParam, order::Int) => SparseMatrixCSC{Float64}

evalbase(p::ChebParam, order::Vector{Int}) => Vector{Matrix{Float64}}
evalbase(p::SplineParam, order::Vector{Int}) => Vector{SparseMatrixCSC{Float64}}
evalbase(p::LinParam, order::Vector{Int}) => Vector{SparseMatrixCSC{Float64}}

evalbase(b::Basis{1}, order::Vector{Int}) => .... can't remember
evalbase(b::Basis{N}, order::Matrix{Int}) => Vector{... can't remembers}
```
=#

#=
Note that each subtype of `BT<:BasisFamily` (with associated `PT<:BasisParam`)
will define the following constructor methods:

```julia
# basis constructors
Basis(::BT, args...)
Basis(::PT)

# node constructor
nodes(::PT)
```

=#

import Base: ==, *, \

using QuantEcon: gridmake, gridmake!, ckron, fix, fix!

import Compat

# types
export BasisFamily, Cheb, Lin, Spline, Basis,
       BasisParams, ChebParams, LinParams, SplineParams,
       AbstractBasisMatrixRep, Tensor, Expanded, Direct,
       BasisMatrix, Interpoland, SplineSparse, RowKron

# functions
export nodes, get_coefs, funfitxy, funfitf, funeval, evalbase,
       derivative_op, row_kron, evaluate, fit!, update_coefs!,
       complete_polynomial, complete_polynomial!, n_complete

include("util.jl")
include("spline_sparse.jl")
include("basis.jl")
include("basis_structure.jl")
include("interp.jl")

# include the rest of the Julian API
include("cheb.jl")
include("spline.jl")
include("lin.jl")

# include comlpete
include("complete.jl")

# deprecations
@deprecate BasisStructure BasisMatrix

end # module


text:
#!/usr/bin/env julia

using VibronicToolkit

using ArgParse
using PyPlot

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "--conf"
        metavar = "FILE"
        help = "path to config file"
        required = true
    "--num-links"
        metavar = "P"
        help = "number of Trotter links"
        arg_type = Int
        required = true
    "--sampling-conf"
        metavar = "FILE"
        help = "path to sampling config file"
    "--sampling-beta"
        metavar = "T"
        help = "sampling reciprocal temperature"
        arg_type = Float64
    "--extent"
        metavar = "q1_min,q1_max,q2_min,q2_max"
        help = "comma-separated extent values"
        required = true
    "--out-path"
        metavar = "FILE"
        help = "path to output file"
        required = true
    "--contour"
        help = "draw contours instead of smooth gradients"
        action = :store_true
end
c = parse_args(ARGS, s, as_symbols=true)

sys = read(c[:conf], DenseSystem)
P = c[:num_links]
if c[:sampling_conf] !== nothing && c[:sampling_beta] !== nothing
    sampling_sys = read(c[:sampling_conf], DiagonalSystem)
    sampling_beta = c[:sampling_beta]
elseif c[:sampling_conf] !== nothing || c[:sampling_beta] !== nothing
    error("Both --sampling-conf and --sampling-beta must be given.")
end
extent = Tuple(parse.(Float64, split(c[:extent], ",")))
out_path = c[:out_path]

# Draw PES.
pes = ground_pes(sys, extent)
if c[:contour]
    cf = contourf(pes, origin="lower", extent=extent, cmap="magma_r")
    contour(pes, origin="lower", extent=extent, colors="k", linestyles="solid", linewidths=0.2)
else
    cf = imshow(pes, origin="lower", extent=extent, cmap="magma_r")
end
colorbar(cf)

if c[:sampling_conf] !== nothing
    # Draw sampling distribution.
    ds = sampling_sys.ds
    ws = weights(sampling_sys, sampling_beta)
    colors = map(plt["matplotlib"]["colors"]["hsv_to_rgb"], [(0.4, w, 0.8) for w in ws])
    scatter(ds[1,:], ds[2,:], marker="+", c=colors)
    for s in 1:size(ds, 2)
        sigma_1 = path_mean_std(sampling_sys, sampling_beta, P, s, 1)
        sigma_2 = path_mean_std(sampling_sys, sampling_beta, P, s, 2)
        ellipse_f = plt["matplotlib"]["patches"]["Ellipse"]
        ellipse = ellipse_f(xy=(ds[1,s], ds[2,s]), width=2sigma_1, height=2sigma_2, edgecolor=colors[s], fc="None", linestyle="dotted")
        gca()["add_patch"](ellipse)
        text(ds[1,s]+0.5, ds[2,s], "$(s)", ha="center", va="center")
    end
end

xlabel(L"q_1")
ylabel(L"q_2")

savefig(out_path, bbox_inches="tight")


text:
#__precompile__()
module CoolProp
using Compat

errcode = Ref{Clong}(0)
const buffer_length = 20000
message_buffer = Array(UInt8, buffer_length)

const inputs_to_get_global_param_string = ["version", "gitrevision", "errstring", "warnstring", "FluidsList", "incompressible_list_pure", "incompressible_list_solution", "mixture_binary_pairs_list", "parameter_list", "predefined_mixtures", "HOME", "cubic_fluids_schema"]

include("CoolPropHighLevel.jl")
include("CoolPropConfig.jl")
include("CoolPropInformation.jl")

# ---------------------------------
# Getter and setter for debug level
# ---------------------------------

"""
    get_debug_level()

Get the debug level.

# Return value
Level The level of the verbosity for the debugging output (0-10) 0: no debgging output
"""
function get_debug_level()
    ccall( (:get_debug_level, "CoolProp"), Cint, () )
end

"""
    set_debug_level(level::Integer)

Set the debug level.

# Arguments
* `level::Integer`: The level of the verbosity for the debugging output (0-10) 0: no debgging output
"""
function set_debug_level(level::Integer) # change ::Int to ::Integer to make set_debug_level(get_debug_level()) works on different machine
    ccall( (:set_debug_level, "CoolProp"), Void, (Cint,), level)
end

include("CoolPropHA.jl")
include("CoolPropLowLevel.jl")

for symorigin = [:PropsSI, :PhaseSI, :K2F, :F2K, :HAPropsSI, :AbstractState_factory, :AbstractState_free, :AbstractState_set_fractions, :AbstractState_update, :AbstractState_keyed_output, :AbstractState_output, :AbstractState_specify_phase, :AbstractState_unspecify_phase, :AbstractState_update_and_common_out, :AbstractState_update_and_1_out, :AbstractState_update_and_5_out, :AbstractState_set_binary_interaction_double, :AbstractState_set_cubic_alpha_C, :AbstractState_set_fluid_parameter_double, :AbstractState_first_saturation_deriv, :AbstractState_first_partial_deriv, :AbstractState_build_phase_envelope, :AbstractState_build_spinodal, :AbstractState_all_critical_points, :AbstractState_get_phase_envelope_data, :AbstractState_get_spinodal_data]
    sym = Symbol(lowercase(string(symorigin)))
    @eval const $sym = $symorigin
    @eval export $sym, $symorigin
end
const set_reference_stateS = set_reference_state
const set_reference_stateD = set_reference_state
const set_config_string = set_config
export set_reference_stateS, set_reference_stateD, set_reference_state
export get_global_param_string, get_parameter_information_string, get_fluid_param_string, get_param_index, get_input_pair_index, set_config
export saturation_ancillary, set_departure_functions, set_config_string, cair_sat
end #module


text:
module QuantumControlBase

include("conditionalthreads.jl")   # submodule ConditionalThreads

include("controlproblem.jl")
export ControlProblem, Objective, MinimalObjective, WeightedObjective

include("propagate.jl")
export propagate_objective, propagate_objectives

include("controls.jl")
export discretize, discretize_on_midpoints
export get_control_parameters, getcontrols
export get_tlist_midpoints

include("gradgen.jl")

include("liouvillian.jl")
export liouvillian

include("infohook.jl")
export chain_infohooks

include("shapes.jl")               # submodule Shapes
include("functionals.jl")          # submodule Functionals
include("weyl_chamber.jl")         # submodule WeylChamber
include("testutils.jl")            # submodule TestUtils
include("saving.jl")               # submodule Saving

include("optimize.jl")
export optimize, @optimize_or_load, optimization_savename, load_optimization
export default_optimization_savename_kwargs

using .TestUtils: optimize_with_dummy_method
optimize(problem, method::Val{:dummymethod}) = optimize_with_dummy_method(problem)

end


text:
function dynamic_loss(x, n_iter)
    for i=1:n_iter
        x = 2 .* x
    end
    return sum(x)
end

function dynamic_loss2(x)
    if sum(x) >= 0
        x = 2 .* x
    else
        x = 3 .* x
    end
    return sum(x)
end

@testset "grad: dynamic" begin
    x = rand(4)
    _, g1 = grad(dynamic_loss, x, 1; dynamic=true)
    _, g2 = grad(dynamic_loss, x, 2; dynamic=true)
    _, g3 = grad(dynamic_loss, x, 3; dynamic=true)
    @test g1[1] == fill(2, 4)
    @test g2[1] == fill(4, 4)
    @test g3[1] == fill(8, 4)

    _, g1 = grad(dynamic_loss2, fill(1.0, 4); dynamic=true)
    _, g2 = grad(dynamic_loss2, fill(-1, 4); dynamic=true)
    @test g1[1] == fill(2, 4)
    @test g2[1] == fill(3, 4)
end


text:
function molecule(::Molecule"H₂O")
  return """
  O 1.2091536548 1.7664118189 -0.0171613972
  H 2.1984800075 1.7977100627  0.0121161719
  H 0.9197881882 2.4580185570  0.6297938830
  """
end
molecule(m::Molecule"water") = molecule(alias(m))
alias(::Molecule"water") = Molecule"H₂O"()


text:
module WeatherWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const Class = @jimport java.lang.Class
const DataContext = @jimport org.orekit.data.DataContext
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const Frame = @jimport org.orekit.frames.Frame
const Geoid = @jimport org.orekit.models.earth.Geoid
const GlobalPressureTemperature2Model = @jimport org.orekit.models.earth.weather.GlobalPressureTemperature2Model
const GlobalPressureTemperatureModel = @jimport org.orekit.models.earth.weather.GlobalPressureTemperatureModel
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object
const TimeScale = @jimport org.orekit.time.TimeScale
const WeatherModel = @jimport org.orekit.models.earth.weather.WeatherModel

export GlobalPressureTemperature2Model
export GlobalPressureTemperatureModel
export WeatherModel
export equals
export get_a
export get_class
export get_pressure
export get_temperature
export get_water_vapor_pressure
export hash_code
export notify
export notify_all
export to_string
export wait
export weather_parameters

include("global_pressure_temperature2_model.jl")
include("global_pressure_temperature_model.jl")
include("weather_model.jl")

end


text:
# Build the dependencies
Mod = @eval module Anon end
Mod.include("build_ZMQ.v4.3.1.jl")
Mod.include("build_libhelicsSharedLib.v2.3.0.jl")


text:
include("../construction.jl")

function compute_obdm(wf::WaveFunction, flavor::Integer, hilbert_space::Array{FullState,1})::Array{DType,2}
    """ Constructs the one-body density matrix for a given wavefunction.
    """
    n_basis::Integer = find_n_basis(hilbert_space)
    lookup_table, inv_lookup_table = make_lookup_table(hilbert_space)

    obdm = zeros(DType, (n_basis, n_basis))
    for n = 1:length(lookup_table)
        state = lookup_table[n]

        # Choose up or down component.
        (state, rest) = f_to_s(state)[[flavor,3-flavor]]
        for k = 1:n_basis
            temp = annihilate(state, k)
            if !isnothing(temp)
                for i = 1:n_basis
                    new_state = create(temp, i)
                    if !isnothing(new_state)
                        full_state = (flavor == 1) ? s_to_f(new_state, rest) : s_to_f(rest, new_state)

                        # This is the index of the final state.
                        final_state = get(inv_lookup_table, full_state, nothing)
                        if !isnothing(final_state)
                            obdm[k,i] += sign(state, k, i) * wf[n]*conj(wf[final_state])
                        end
                    end
                end
            end
        end
    end

    return obdm
end

function compute_tbdm_up_down(wf::WaveFunction, lookup_table::LookupDict, inv_lookup_table::InvLookupDict, n_basis::Integer)::Array{DType,4}
    # The one-body density-matrix in the
    tbdm = zeros(DType, (n_basis, n_basis, n_basis, n_basis))
    for n = 1:length(lookup_table)
        s_up, s_down = f_to_s(lookup_table[n])
        for k = 1:n_basis
            temp_up = annihilate(s_up, k)
            if !isnothing(temp_up)
                for i = 1:n_basis
                    new_up = create(temp_up, i)
                    if !isnothing(new_up)
                        for l = 1:n_basis
                            temp_down = annihilate(s_down, l)
                            if !isnothing(temp_down)
                                for j = 1:n_basis
                                    new_down = create(temp_down, j)
                                    if !isnothing(new_down)
                                        final_state = inv_lookup_table[s_to_f(new_up,new_down)]
                                        tbdm[k,i,l,j] += sign(s_up,k,i)*sign(s_down,l,j) * wf[n]*conj(wf[final_state])
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    return tbdm
end


text:
function odd_iter_count(iteration, x, value)
    return floor(log(6, (value * 2 ^ iteration) / (x)))
end

function sum_array(arr, start_i, end_i)::Float64
    sum::Float64 = 0.0
    for i in start_i:end_i
        sum += arr[i]
    end
    sum
end

function collatz(x)
    if x % 2 == 0
        return x / 2
    end
    return 3 * x + 1
end

function collatz_conj(x)
    n_set = [0]
    even_count = 0
    odd_count = 0
    iteration = 0
    while x != 1
        if x % 2 == 0
            even_count += 1
            n_set[end] += 1
        else
            append!(n_set, 0)
            odd_count += 1
        end
        x = collatz(x)
        iteration += 1
    end

    (reverse!(n_set), iteration, even_count, odd_count)
end

function inverse_collatz_conj(n; Type=Float64)
    f = length(n)
    inner_sum::Float64 = 0
    for i in 2:f
        inner_sum += (2 ^ sum_array(n, i, f - 1)) * 3 ^ (i - 2)
    end
    end_ratio::Type = ((2 ^ n[end]) / (3 ^ (f - 1)))
    return round(end_ratio * (2 ^ sum_array(n, 1, f - 1) - inner_sum))
end

g(x) = (x-1)/3
z(x, m) = m * g(x)

function apply(f, x, m, n=1)
    for _ in 1:n
        x = f(x, m)
    end
    return x
end

function collect_collatz(size)
    xvals = zeros(Float64, size)
    yvals = zeros(Float64, size)

    for idx in 1:size
        i = 2 * idx - 1
        n = collatz_conj(i)[1]
        f = length(n)
        pow = sum_array(n, 1, f - 1)

        xvals[idx] = idx
        yvals[idx] = pow
    end
    (xvals, yvals)
end


text:
const ωe = ROT_VELOCITY

"""
    six_dof_ecef_quaternion_fixed_mass(
        x, mass, inertia, forces, moments, h=[0.0, 0.0, 0.0], k=0.0, ellipsoid=WGS84
        )

Six degrees of freedom dynamic system using quaternions for attitude
representation and assuming fixed mass.

Ellipsoidal Earth Model is used and the ECEF reference frame is considered
inertial.

The effects of the angular momentum produced by spinning rotors is taken into
account with the optional argument `h`.

# Arguments
- `x::13-element Array{Number,1}`: state vector.
    u, v, w: inertial linear velocity expressed in body axis. (m/s)
    p, q, r: inertial rotatinal velocity expressed in body axis. (rad/s)
    q0, q1, q2, q3: attitude given by quaternions.
    xe, ye, ze: position wrt the inertial system origin expressed in Earth Axis. (m)
- `mass::Number`: total mass of the aircraft (kg)
- `inertia::3×3 Array{Number,2}`: inertia tensor (kg·m²)
- `forces::3-element Array{Number,1}`: total forces expressed in body axis. (N)
- `moments::3-element Array{Number,1}`: total moments expressed in body axis.(N·m)
- `h::3-element Array{Number,1}`: Additional angular momentum contributions such
  as those coming from spinning rotors (kg·m²/s).
- `k::Number`: orthonormality error factor.
- `ellipsoid::Ellipsoid`: ellipsoid model to be used.

# Returns
- `x_dot`: state vector derivative according to the equation of motion,
    inertial properties and applied forces and moments.

# Notes
- See [1] (page 41, formula 1.4-23) or [2] (page 372, formulas 10-11 to 10-14)
  for more information on quaternions <-> Euler angles conversions.
- Orthonormality error factor is related to a numerical stability artifact used
  in angular kinematic equations. Let λ = k * (1 - q0² - q1² - q2² - q3²) be
  the orthonormality error. The term k·λ·q is added to the angular kinematic
  equations in order to reduce the numerical integration error. According to
  reference [2] k·Δt ≤ 1. See [2] (page 372) for more information on
  orthonormality error factor.
- Implementation based on [1]. However, [2] can also be read.

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 45, formula 1.5-1)

- [2] Zipfel, P. H. (2007). Modeling and simulation of aerospace vehicle
 dynamics. American Institute of Aeronautics and Astronautics.
 (page 396, figure 10.6)
"""
function six_dof_ecef_quaternion_fixed_mass(
    x, mass, inertia, forces, moments, h=[0.0, 0.0, 0.0], k=0.0, ellipsoid=WGS84
    )

    vb = x[1:3]    # u, v, w
    ωb = x[4:6]    # p, q, r
    quat = x[7:10]   # q0, q1, q2, q3
    pos  = x[11:13]  # px, py, pz (ecef)

    Fb = forces
    J  = inertia
    Tb = moments

    Ωe = [0.0   ωe    0.0;
          -ωe   0.0   0.0;
          0.0   0.0   0.0]

    p, q, r = ωb

    ΩB = [0.0   -r    q;
           r    0.0  -p;
          -q     p   0.0]

    Ωq = [0.0    p    q    r;
          -p    0.0  -r    q;
          -q     r   0.0  -p;
          -r    -q    p   0.0]

    lat, lon, height = ecef2llh(pos...; ellipsoid=ellipsoid)
    B = rot_matrix_body2ecef(lat, lon, quat...)

    # Linear kinematic equations
    pos_dot = Ωe * pos + B' * vb
    # Linear momentum equations
    # Note that gravity (B g(p) - B Ωe²) p is included in forces (gravity + prop + aero)
    vb_dot = -(ΩB + B*Ωe) * vb + Fb / mass
    # Angular momentum equations
    Jinv = inv(J)
    # ωb_dot = (-Jinv * ΩB * J) * ωb + Jinv * Tb  # no gyroscopic effects
    ωb_dot = Jinv * (-ΩB * (J * ωb + h) + Tb)
    # Angular Kinematic equations
    # Normalization λ*t < 1 (See Zipfel Chapter 4.3.3.4 p.126)
    # Zipfel, P. H. (2007). Modeling and simulation of aerospace vehicle dynamics.
    # American Institute of Aeronautics and Astronautics.
    λ = k * (1.0 - quat' * quat)
    quat_dot = (-0.5 * Ωq * quat) + λ * quat

    [vb_dot..., ωb_dot..., quat..., pos_dot...]
end


text:
function Grid=SphereGrid(OrientFace,Param)
Grid.nBar=[ 0  1   0   1
           -1  0  -1   0];
  Grid.Dim=3;
  Grid.Type='Quad';
  Grid.Rad=Param.RadEarth;
  Grid.Form='Sphere';
  lat1=0;
  lat2=lat1+pi/1000;
  lon1=0;
  lon2=lon1+pi/1000;

  NumNodes=4;
  Nodes(1:NumNodes)=Node([0 0 0],0);
  NodeNumber=1;
  Nodes(NodeNumber)=Node(sphere2cart(lat1,lon1,Grid.Rad),NodeNumber);
  NodeNumber=2;
  Nodes(NodeNumber)=Node(sphere2cart(lat2,lon1,Grid.Rad),NodeNumber);
  NodeNumber=3;
  Nodes(NodeNumber)=Node(sphere2cart(lat2,lon2,Grid.Rad),NodeNumber);
  NodeNumber=4;
  Nodes(NodeNumber)=Node(sphere2cart(lat1,lon2,Grid.Rad),NodeNumber);
  Grid.Nodes=Nodes;

  NumEdges=4;
  NumEdgesI=4;
  Edges(1:NumEdges)=Edge([1 2],Grid,0,0,'');
  EdgeNumber=1;
  Edges(EdgeNumber)=Edge([1 2],Grid,EdgeNumber,EdgeNumber,'');
  EdgeNumber=2;
  Edges(EdgeNumber)=Edge([2 3],Grid,EdgeNumber,EdgeNumber,'');
  EdgeNumber=3;
  Edges(EdgeNumber)=Edge([3 4],Grid,EdgeNumber,EdgeNumber,'');
  EdgeNumber=4;
  Edges(EdgeNumber)=Edge([4 1],Grid,EdgeNumber,EdgeNumber,'');
  Grid.Edges=Edges;

  NumFaces=1;
  FaceNumber=1;
  [Faces(FaceNumber),Grid]=Face([1 2 3 4],Grid,FaceNumber,'',OrientFace);
  Grid.Faces=Faces;

  Grid.NumNodes=size(Grid.Nodes,2);
  Grid.NumEdges=size(Grid.Edges,2);
  Grid.NumEdgesI=size(Grid.Edges,2);
  Grid.NumEdgesB=0;
  Grid.NumFaces=size(Grid.Faces,2);
  Grid.Dim=3;
  Grid=Orientation(Grid);
  Grid=Renumbering(Grid);
  Grid=FacesInNodes(Grid);

end


text:
"""
    struct CylindricalImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}

Impurity density model which assumes a linear gradient in impurity density in each spatial dimension
of a cylindrical coordinate system.

## Fields
* `offsets::NTuple{3,T}`: impurity density values at the origin of each dimension.
* `gradients::NTuple{3,T}`: linear slopes in `r` and `z` direction.

## Definition in Configuration File

A `CylindricalImpurityDensity` is defined in the configuration file through the field `impurity_density`
(of a `passive` or `surrounding`) with `name: cylindrical` and optional fields `r` and `z`
that can each contain `init` for initial values at 0 and `gradient` for gradients in that dimension.

An example definition of a cylindrical impurity density looks like this:
```yaml
impurity_density:
  name: cylindrical
  r:  # impurity profile with linear gradient in r
    init: 1.0e10     # 1/m³
    gradient: 1.0e11 # 1/m⁴
```
"""
struct CylindricalImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
    offsets::NTuple{3, T}
    gradients::NTuple{3, T}
end

function ImpurityDensity(T::DataType, t::Val{:cylindrical}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, input_units::NamedTuple)
    offsets, gradients = zeros(T,3), zeros(T,3)
    density_unit = input_units.length^(-3)
    density_gradient_unit = input_units.length^(-4)
    if prod(map(k -> k in ["r","z"], collect(keys(dict)))) @warn "Only r and z are supported in the cylindrical impurity density model.\nChange the impurity density model in the config file or remove all other entries." end
    if haskey(dict, "r")
        if haskey(dict["r"], "init")     offsets[1]   = _parse_value(T, dict["r"]["init"], density_unit) end
        if haskey(dict["r"], "gradient") gradients[1] = _parse_value(T, dict["r"]["gradient"], density_gradient_unit) end
    end
    if haskey(dict, "z")
        if haskey(dict["z"], "init")     offsets[3]   = _parse_value(T, dict["z"]["init"], density_unit) end
        if haskey(dict["z"], "gradient") gradients[3] = _parse_value(T, dict["z"]["gradient"], density_gradient_unit) end
    end
    CylindricalImpurityDensity{T}( NTuple{3, T}(offsets), NTuple{3, T}(gradients) )
end

function get_impurity_density(lcdm::CylindricalImpurityDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    pt::CylindricalPoint{T} = CylindricalPoint(pt)
    ρ::T = 0
    for i in eachindex(lcdm.offsets)
        ρ += (lcdm.offsets[i] + pt[i] * lcdm.gradients[i])
    end
    return ρ
end


text:
using CBLS
using ConstraintModels
using Dictionaries
using JuMP
using LocalSearchSolvers
using MathOptInterface
using Test

@testset "ConstraintModels.jl" begin
    include("instances.jl")
    include("raw_solver.jl")
    include("MOI_wrapper.jl")
    include("JuMP.jl")
end


text:
using QuDirac
using Base.Test

default_inner(UndefinedInner())

include("constest.jl")
include("abstractinnertests.jl")
include("generaltests.jl")

default_inner(KroneckerDelta())

include("constest.jl")
include("orthotests.jl")
include("generaltests.jl")
include("laddertests.jl")


text:
import Distributions
import Statistics

abstract type TheoreticalFittedUncertainScalar <: TheoreticalDistributionScalarValue end

Broadcast.broadcastable(uv::TheoreticalFittedUncertainScalar) = Ref(uv.distribution)

"""
UncertainScalarTheoreticalFit

An empirical value represented by a distribution estimated from actual data.

## Fields
- **`distribution`** The distribution describing the value.
- **`values`**: The values from which `distribution` is estimated.
"""
struct UncertainScalarTheoreticalFit{D <: Distribution, T} <: TheoreticalFittedUncertainScalar
    distribution::FittedDistribution{D} # S may be Continuous or Discrete
    values::AbstractVector{T}
end

"""
    ConstrainedUncertainScalarTheoreticalFit

An empirical value represented by a distribution estimated from actual data.

## Fields
- **`distribution`** The truncated version of the distribution describing the
    value.
- **`values`**: The values from which the original distribution was estimated.
"""
struct ConstrainedUncertainScalarTheoreticalFit{D <: Distribution, T} <: TheoreticalFittedUncertainScalar
    distribution::FittedDistribution{D} # S may be Continuous or Discrete
    values::AbstractVector{T}
end

""" Truncate a fitted distribution. """
Distributions.truncated(fd::FittedDistribution, lower, upper) =
    Distributions.truncated(fd.distribution, lower, upper)

Base.rand(fd::UncertainScalarTheoreticalFit) = rand(fd.distribution.distribution)
Base.rand(fd::UncertainScalarTheoreticalFit, n::Int) = rand(fd.distribution.distribution, n)

# For the fitted distributions, we need to access the fitted distribution's distribution
Distributions.pdf(fd::UncertainScalarTheoreticalFit, x) = pdf(fd.distribution.distribution, x)
StatsBase.mode(uv::UncertainScalarTheoreticalFit) = mode(uv.distribution.distribution)
Statistics.mean(uv::UncertainScalarTheoreticalFit) = mean(uv.distribution.distribution)
Statistics.median(uv::UncertainScalarTheoreticalFit) = median(uv.distribution.distribution)
Statistics.quantile(uv::UncertainScalarTheoreticalFit, q) = quantile(uv.distribution.distribution, q)
Statistics.std(uv::UncertainScalarTheoreticalFit) = std(uv.distribution.distribution)
Statistics.var(uv::UncertainScalarTheoreticalFit) = var(uv.distribution.distribution)

export
TheoreticalFittedUncertainScalar,
UncertainScalarTheoreticalFit,
ConstrainedUncertainScalarTheoreticalFit


text:
module StaticPermutations

export AbstractPermutation, Permutation, NoPermutation

export
    identity_permutation,
    isidentity,
    append,
    prepend

import Base: ==, *, /, \

include("types.jl")
include("operations.jl")
include("arrays.jl")

end


text:
module DeBruijn

using Reexport
@reexport using ..LambdaCalculus

import Base: length, show
import ..LambdaCalculus: boundvartype, freevartype, reify

export App,
    Abs,
    Index,
    Term,
    Var

"""Lambda terms using [De Bruijn indexing](De Bruijn indexing), built using only the following rule:

    <Term> := <Number>            (variable)
            | λ <Term>            (abstraction)
            | (<Term> <Term>)     (application)

Each De Bruijn index is a natural number that represents an occurrence of a variable in a ``λ``-term,
and denotes the number of binders that are in scope between that occurrence and its corresponding
binder.

!!! warning "Important note"

    Since we are in Julia, indices start at ``1``.
"""
abstract type Term <: AbstractTerm end

struct Var <: Term
    index::Index

    Var(index) = index ≥ 1 ? new(index) : error("index must be at least 1")
end

struct Abs <: Term
    body::Term
end

struct App <: Term
    car::Term
    cdr::Term
end

show(io::IO, t::Abs) = print(io, "(λ", ".", t.body, ")")
show(io::IO, t::App) = print(io, "(", t.car, " ", t.cdr, ")")
show(io::IO, t::Var) = print(io, "⟨", t.index, "⟩")

length(t::Var) = t.index + 1
length(t::App) = length(t.car) + length(t.cdr) + 2
length(t::Abs) = length(t.body) + 2

reify(v::Var) = :(Var($(v.index)))
reify(t::Abs) = :(Abs($(reify(t.body))))
reify(t::App) = :(App($(reify(t.car)), $(reify(t.cdr))))

freevartype(::Type{<:Term}) = Var
boundvartype(::Type{<:Term}) = Var

include("syntactic.jl")
include("meta.jl")
include("sampling.jl")

end # module DeBruijn


text:
# General api
"""
    formula(trm::TableRegressionModel)
    formula(model::MixedModel)

Unify formula api.
"""
formula(model) = error("formula is not defined for $(typeof(model)).")

"""
    nestedmodels(trm::TableRegressionModel{<: LinearModel}; null::Bool = false, <keyword arguments>)
    nestedmodels(trm::TableRegressionModel{<: GeneralizedLinearModel}; null::Bool = false, <keyword arguments>)
    nestedmodels(model::LinearMixedModel; null::Bool = false, <keyword arguments>)

    nestedmodels(::Type{LinearModel}, formula, data; null::Bool = true, <keyword arguments>)
    nestedmodels(::Type{GeneralizedLinearModel}, formula, data, distr::UnivariateDistribution, link::Link = canonicallink(d); null::Bool = true, <keyword arguments>)
    nestedmodels(::Type{LinearMixedModel}, f::FormulaTerm, tbl; null::Bool = true, wts = [], contrasts = Dict{Symbol, Any}(), verbose::Bool = false, REML::Bool = false)

Generate nested models from a model or formula and data. \n
The null model will be a model with at least one factor (including intercept) if the link function does not allow factors to be 0 (factors in denominators). \n
* `InverseLink` for `Gamma`
* `InverseSquareLink` for `InverseGaussian`
* `LinearModel` fitted with `CholeskyPivoted` when `dropcollinear = true`
Otherwise, it will be an empty model.
"""
nestedmodels(model) = error("nestedmodels is not defined for $(typeof(model)).")

# implement drop1/add1 in R?

const FixDispDist = Union{Bernoulli, Binomial, Poisson}
"""
    canonicalgoodnessoffit(::FixDispDist) = LRT
    canonicalgoodnessoffit(::UnivariateDistribution) = FTest

    const FixDispDist = Union{Bernoulli, Binomial, Poisson}

Return LRT if the distribution has a fixed dispersion.
"""
canonicalgoodnessoffit(::FixDispDist) = LRT
canonicalgoodnessoffit(::UnivariateDistribution) = FTest

"""
    anova(<models>...; test::Type{<: GoodnessOfFit})

Analysis of variance. \n
Return `AnovaResult{M, test, N}`. See `AnovaResult` for details.

* `models`: model objects
    1. `TableRegressionModel{<: LinearModel}` fit by `GLM.lm`
    2. `TableRegressionModel{<: GeneralizedLinearModel}` fit by `GLM.glm`
    3. `LinearMixedModel` fit by `MixedAnova.lme` or `fit(LinearMixedModel, ...)`
    If mutiple models are provided, they should be nested and the last one is the most saturated.
* `test`: test statistics for goodness of fit. Available tests are `LikelihoodRatioTest` (`LRT`) and `FTest`. \n
    If no test argument is provided, the function will automatically determine based on the model type:
    1. `TableRegressionModel{<: LinearModel}`: `FTest`.
    2. `TableRegressionModel{<: GeneralizedLinearModel}`: based on distribution function, see `canonicalgoodnessoffit`.
    3. `LinearMixedModel`: `FTest` for one model fit; `LRT` for nested models.

When multiple models are provided:
* `check`: allows to check if models are nested. Defalut value is true. Some checkers are not implemented now.
* `isnested`: true when models are checked as nested (manually or automatically). Defalut value is false.

For fitting new models and conducting anova at the same time, \n
see `anova_lm` for `LinearModel`, `anova_glm` for `GeneralizedLinearModel`, `anova_lme` for `LinearMixedModel`.
"""
anova(model) = error("anova is not defined for $(typeof(model)).")

"""
    anova(::Type{FTest}, <model>; kwargs...)
    anova(::Type{FTest}, <models>...; kwargs...)

Analysis of Variance by F-test. \n
Return `AnovaResult{M, FTest, N}`. See `AnovaResult` for details.

* `type` specifies type of anova:
    1. One `LinearModel` or `GeneralizedLinearModel`: 1, 2, 3 are valid
    2. One `LinearMixedModel`: 1, 3 are valid.
    3. Others: only 1 is valid.
* `adjust_sigma` determines if adjusting to REML when `LinearMixedModel` is fit by maximum likelihood.  \n
    The result will be slightly deviated from that of model fit by REML.
"""
anova(::Type{FTest}, model) = error("anova by F-test is not defined for $(typeof(model)).")

"""
    anova(::Type{LRT}, <model>; kwargs...)
    anova(::Type{LRT}, <models>...; kwargs...)

Analysis of Variance by likelihood-ratio test. \n
Return `AnovaResult{M, LRT, N}`. See `AnovaResult` for details.
"""
anova(::Type{LRT}, model) = error("anova by likelihood-ratio test is not defined for $(typeof(model)).")

# across different kind of models
function _lrt_nested(models::NTuple{N, RegressionModel}, df, dev, σ²; nestedwarn::Bool = true) where N
    nestedwarn && @warn "Could not check whether models are nested: results may not be meaningful"
    Δdf = _diff(df)
    Δdev = _diffn(dev)
    lrstat = Δdev ./ σ²
    pval = map(zip(Δdf, lrstat)) do (dof, lr)
        lr > 0 ? ccdf(Chisq(dof), lr) : NaN
    end
    AnovaResult{LRT}(models, 1, df, dev, (NaN, lrstat...), (NaN, pval...), NamedTuple())
end

"""
    isnullable(::CholeskyPivoted)
    isnullable(::Cholesky
    isnullable(::InverseLink)
    isnullable(::InverseSquareLink)
    isnullable(::Link)
    isnullable(::LinearModel)
    isnullable(model::GeneralizedLinearModel)
    isnullable(::LinearMixedModel)

Return `true` if empty model can be fitted.
"""
isnullable(m) = error("isnullable is not defined for $(typeof(m)).")

"""
    nobs(aov::AnovaResult{<: Tuple})
    nobs(aov::AnovaResult)

Apply `nobs` to all models in `aov.model`
"""
nobs(aov::AnovaResult{<: Tuple}) = nobs.(aov.model)
nobs(aov::AnovaResult) = nobs(aov.model)

"""
    dof(aov::AnovaResult)

Degree of freedom of models or factors.
"""
dof(aov::AnovaResult) = aov.dof

"""
    dof_residual(aov::AnovaResult{<: Tuple})
    dof_residual(aov::AnovaResult)
    dof_residual(aov::AnovaResult{<: MixedModel, FTest})

Degree of freedom of residuals.
Default is applying `dof_residual` to models in `aov.model`.
For `MixedModels` applying `FTest`, it is calculated by between-within method. See `calcdof` for details.
"""
dof_residual(aov::AnovaResult{<: Tuple}) = dof_residual.(aov.model)
dof_residual(aov::AnovaResult) = dof_residual(aov.model)

"""
    deviance(aov::AnovaResult)

Return the stored devaince. The value repressents different statistics for different models and tests.
1. `LinearModel`: Sum of Squares.
2. `GeneralizedLinearModel`: `deviance(model)`
3. `LinearMixedModel`: `NaN` when applying `FTest`; `-2loglikelihood(model) == deviance(model)` when applying `LRT`.
When `LinearModel` is compared to `LinearMixedModel`, the deviance is alternatively `-2loglikelihood(model)`.
"""
deviance(aov::AnovaResult) = aov.deviance

"""
    teststat(aov::AnovaResult)

Values of test statiscics of `anova`.
"""
teststat(aov::AnovaResult) = aov.teststat

"""
    teststat(aov::AnovaResult)

P-values of test statiscics of `anova`.
"""
pval(aov::AnovaResult) = aov.pval

"""
    anova_test(::AnovaResult)

Test statiscics of `anova`.
"""
anova_test(::AnovaResult{M, T}) where {M, T <: GoodnessOfFit} = T

"""
    anova_type(aov::AnovaResult)

Type of `anova`.
"""
anova_type(aov::AnovaResult) = aov.type

# Calculate dof from assign
function StatsBase.dof(v::Vector{Int})
    dofv = zeros(Int, v[end])
    prev = 1
    ind = 1
    n = length(v)
    while ind <= n
        v[ind] == prev || (prev = v[ind])
        dofv[prev] += 1
        ind += 1
    end
    Int.(dofv)
end


text:
# task file
const _GW_TASK_FILE_NAME = "task.toml"
_task_file(taskdir::String) = joinpath(taskdir, _GW_TASK_FILE_NAME)
_task_file(gwt::GWTask) = _task_file(task_dir(gwt))

# tid
const _GW_TASK_TID_KEY = "tid"

# exptime
const _GW_TASK_EXPTIME_KEY = "exptime"
const _GW_DFLT_TASK_VTIME = 120.0

# status
const _GW_TASK_STATUS_KEY = "runstate"
const _GW_TASK_PENDING_STATUS = "pending"
const _GW_TASK_RUNNING_STATUS = "running"
const _GW_TASK_DONE_STATUS = "done"
const _GW_TASK_ERROR_STATUS = "error"
const _GW_TASK_UNKNOWN_STATUS = "unknown"

_task_toml(gwt::GWTask) = get!(gwt, _GW_TASK_FILE_NAME) do
    Dict{String, Any}(
        _GW_TASK_TID_KEY => task_id(gwt)
    )
end

function _write_task_toml!(gwt::GWTask;
        extime = _GW_TASK_TID_KEY,
        status = _GW_TASK_PENDING_STATUS,
        desc = ""
    )

    toml = _task_toml(gwt)
    toml[_GW_TASK_EXPTIME_KEY] = extime
    toml[_GW_TASK_STATUS_KEY] = status
    toml[_GW_TASK_DESC_KEY] = desc

    tfile = _task_file(gwt)
    _write_toml(tfile, toml)

    return toml
end

function _read_task_toml!(gwt::GWTask)
    tfile = _task_file(gwt)
    toml = _read_toml(tfile)
    merge!(_task_toml(gwt), toml)
    return toml
end

function _get_task_status(gwt::GWTask)
    toml = _task_toml(gwt)
    get!(toml, _GW_TASK_STATUS_KEY, _GW_TASK_UNKNOWN_STATUS)
end

function _up_task_status!(gwt::GWTask, status)
    _read_task_toml!(gwt)
    toml = _task_toml(gwt)
    toml[_GW_TASK_STATUS_KEY] = status
    _write_toml(_task_file(gwt), toml)
end

_is_pending_task(gwt::GWTask) = (_get_task_status(gwt) == _GW_TASK_PENDING_STATUS)

function _get_task_extime(gwt::GWTask)
    toml = _task_toml(gwt)
    get!(toml, _GW_TASK_EXPTIME_KEY, -1.0)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
# * Atomic System of Equations

"""
    AtomicEquations(atom, equations, integrals)

Structure representing the (e.g. Hartree–Fock) `equations` for `atom`,
along with all `integrals` that are shared between the `equations`.
"""
mutable struct AtomicEquations{T, B<:AbstractQuasiMatrix,
                               O<:AbstractOrbital,
                               A<:Atom{T,B,O},
                               Equations<:AbstractVector{<:AtomicOrbitalEquation}}
    atom::A
    equations::Equations
    integrals::Vector{OrbitalIntegral}
    observables::Dict{Symbol,<:Observable}
end

# Iteration interface
Base.length(hfeqs::AtomicEquations) = length(hfeqs.equations)
Base.iterate(iter::AtomicEquations, args...) = iterate(iter.equations, args...)

"""
    update!(equations::AtomicEquations[, atom::Atom])

Recompute all integrals using the current values for the radial
orbitals (optionally specifying which `atom` from which the orbitals
are taken).
"""
function SCF.update!(equations::AtomicEquations, args...; kwargs...)
    Threads.@threads for integral in equations.integrals
        SCF.update!(integral, args...; kwargs...)
    end
    # Updating the operators is apparently not thread-safe at the
    # moment. Why?
    for eq in equations.equations
        update!(SCF.hamiltonian(eq), args...; kwargs...)
    end
end

"""
    energy_matrix!(H, hfeqs::AtomicEquations[, which=:energy])

Compute the energy matrix by computing the energy observable and
storing it in `H`. Requires that `hfeqs` has the `:energy` and
`:kinetic_energy` [`Observable`](@ref)s registered (this is the
default).
"""
function SCF.energy_matrix!(H::HM, hfeqs::AtomicEquations,
                            which::Symbol=:total_energy) where {HM<:AbstractMatrix}
    observable = hfeqs.observables[which]
    observe!(H, observable)
    H
end

# * Orbital symmetries
"""
    find_symmetries(orbitals)

Group all orbitals according to their symmetries, e.g. ℓ for
`Orbital`s. This is used to determine which off-diagonal Lagrange
multipliers are necessary to maintain orthogonality.
"""
find_symmetries(orbitals::Vector{O}) where {O<:AbstractOrbital} =
    merge!(vcat, [Dict(symmetry(orb) => O[orb])
                  for orb in orbitals]...)

SCF.symmetries(atom::Atom) =
    map(values(find_symmetries(atom.orbitals))) do orbitals
        orbital_index.(Ref(atom), orbitals)
    end

# * Setup orbital equations

function get_operator(::FieldFreeOneBodyHamiltonian, atom::Atom,
                      orbital::aO, source_orbital::bO; kwargs...) where {aO,bO}
    orbital == source_orbital ||
        throw(ArgumentError("FieldFreeOneBodyHamiltonian not pertaining to $(orbital)"))
    AtomicOneBodyHamiltonian(atom, orbital)
end

for (i,HT) in enumerate([:KineticEnergyHamiltonian, :PotentialEnergyHamiltonian])
    @eval begin
        get_operator(::$HT, atom::Atom, orbital::aO, ::bO; kwargs...) where {aO,bO} =
            AtomicOneBodyHamiltonian(one_body_hamiltonian(Tuple, atom, orbital)[$i],
                                     orbital)
    end
end

function get_operator(op::CoulombPotentialMultipole, atom::Atom,
                      orbital::aO, ::bO; kwargs...) where {aO,bO}
    a,b = op.a[1],op.b[1]
    HFPotential(:exchange, op.o.k, a, a, atom, op.o.g; kwargs...)
end

get_operator(top::IdentityOperator, atom::Atom,
             ::aO, source_orbital::bO; kwargs...) where {aO,bO} =
                 SourceTerm(top, source_orbital, view(atom, source_orbital))
SCF.update!(::IdentityOperator; kwargs...) = nothing
SCF.update!(::IdentityOperator, ::Atom; kwargs...) = nothing

function get_operator(op::RadialOperator, atom::Atom, orbital::aO, source_orbital::bO; kwargs...) where {aO, bO}
    if orbital == source_orbital
        # In this case, the operator is diagonal in orbital space,
        # i.e. it maps an orbital onto itself.
        op
    else
        # In this case, the operator maps from source_orbital to
        # orbital.
        SourceTerm(op, source_orbital, view(atom, source_orbital))
    end
end

function get_operator(M::AbstractMatrix, atom::Atom, a::aO, b::bO; kwargs...) where {aO, bO}
    R = radial_basis(atom)
    get_operator(R*M*R', atom, a, b; kwargs...)
end

# RadialOperators (which are built from matrices), are independent of
# the atom.
SCF.update!(::RadialOperator; kwargs...) = nothing
SCF.update!(::RadialOperator, ::Atom; kwargs...) = nothing

# function get_operator(top::Projector, ::Atom, orbital::aO, source_orbital::bO) where {aO,bO}
#     orbital == source_orbital || return 0
#     SourceTerm(IdentityOperator{1}(), source_orbital,
#                top.ϕs[findfirst(isequal(source_orbital), top.orbitals)])
# end

get_operator(op::QO, atom::Atom, ::aO, ::bO; kwargs...) where {QO,aO,bO} =
    throw(ArgumentError("Unsupported operator $op"))

"""
    generate_atomic_orbital_equations(atom::Atom, eqs::MCEquationSystem,
                                      integrals, integral_map)

For each variationally derived orbital equation in `eqs`, generate the
corresponding [`AtomicOrbitalEquation`](@ref).
"""
function generate_atomic_orbital_equations(atom::Atom{T,B,O}, eqs::MCEquationSystem,
                                           integrals::Vector,
                                           integral_map::Dict{Any,Int},
                                           symmetries::Dict;
                                           verbosity=0,
                                           kwargs...) where {T,B,O}
    map(eqs.equations) do equation
        orbital = equation.orbital
        terms = Vector{OrbitalHamiltonianTerm{O,O,T}}()

        for (integral,equation_terms) in equation.terms
            if integral > 0
                operator = get_integral(integrals, integral_map, eqs.integrals[integral])
                # We first add all terms with operators diagonal in
                # orbital space.
                pushterms!(terms, operator, filter(t -> t.source_orbital == orbital, equation_terms),
                           integrals, integral_map, eqs.integrals)
                # We then add those terms are that are off-diagonal in
                # orbital space.
                for t in filter(t -> t.source_orbital ≠ orbital, equation_terms)
                    pushterms!(terms, SourceTerm(operator, t.source_orbital, view(atom, t.source_orbital)), [t],
                               integrals, integral_map, eqs.integrals)
                end
            else
                for t in equation_terms
                    operator = get_operator(t.operator, atom, orbital, t.source_orbital; kwargs...)
                    iszero(operator) && continue

                    pushterms!(terms, operator, [t],
                               integrals, integral_map, eqs.integrals)
                end
            end
        end

        # Find all other orbitals of the same symmetry as the current
        # one. These will be used to create a projector, that projects
        # out their components.
        #
        # TODO: Think of what the non-orthogonalities due to
        # `overlaps` imply for the Lagrange multipliers/projectors.
        symmetry_orbitals = filter(!isequal(orbital), symmetries[symmetry(orbital)])
        verbosity > 1 && println("Symmetry: ", symmetry_orbitals)

        AtomicOrbitalEquation(atom, equation, orbital, terms, symmetry_orbitals)
    end
end

atomic_hamiltonian(::Atom{T,B,O,TC,CV,P}) where {T,B,O,TC,CV,P<:AbstractPotential} =
    FieldFreeOneBodyHamiltonian() + CoulombInteraction()

"""
    diff(atom; H=atomic_hamiltonian(atom), overlaps=[], selector=outsidecoremodel, verbosity=0)

Differentiate the energy expression of the Hamiltonian `H` associated
with the `atom`'s configurations(s) with respect to the atomic
orbitals to derive the Hartree–Fock equations for the orbitals.

By default, the Hamiltonian
`H=FieldFreeOneBodyHamiltonian()+CoulombInteraction()`.

Non-orthogonality between orbitals can be specified by providing
`OrbitalOverlap`s between these pairs. Only those electrons not
modelled by `atom.potential` of each configuration are considered for
generating the energy expression, this can be changed by choosing
another value for `selector`.
"""
function Base.diff(atom::Atom,
                   energy_expression::EnergyExpression;
                   H::QuantumOperator=atomic_hamiltonian(atom),
                   overlaps::Vector{<:OrbitalOverlap}=OrbitalOverlap[],
                   selector::Function=cfg -> outsidecoremodel(cfg, atom.potential),
                   configurations = selector.(atom.configurations),
                   orbitals = unique_orbitals(configurations),
                   symmetries = find_symmetries(orbitals),
                   observables::Dict{Symbol,Tuple{<:QuantumOperator,Bool}} =
                   Dict{Symbol,Tuple{<:QuantumOperator,Bool}}(
                       :total_energy => (H,false),
                       :double_counted_energy => (H,true),
                       :kinetic_energy => (KineticEnergyHamiltonian(),false),
                   ),
                   modify_eoms!::Function = eqs -> nothing,
                   modify_integrals!::Function = (eqs,integrals,integral_map) -> nothing,
                   modify_equations!::Function = hfeqs -> nothing,
                   verbosity=0)
    eqs = diff(energy_expression, conj.(orbitals))
    modify_eoms!(eqs)

    if verbosity > 0
        println("Energy expression:")
        display(energy_expression)
        println()

        println("Hartree–Fock equations:")
        display(eqs)
        println()

        println("Symmetries:")
        display(symmetries)
        println()

        if verbosity > 1
            println("Common integrals:")
            display(eqs.integrals)
            println()
        end
    end

    integrals = Vector{OrbitalIntegral}()
    integral_map = Dict{Any,Int}()
    poisson_cache = Dict{Int,CoulombIntegrals.PoissonCache}()
    append_common_integrals!(integrals, integral_map,
                             atom, eqs.integrals,
                             poisson_cache=poisson_cache)
    modify_integrals!(eqs, integrals, integral_map)

    hfeqs = generate_atomic_orbital_equations(atom, eqs,
                                              integrals, integral_map,
                                              symmetries; verbosity=verbosity,
                                              poisson_cache=poisson_cache)
    modify_equations!(hfeqs)

    observables = map(collect(pairs(observables))) do (k,(operator,double_counted))
        verbosity > 2 && println("Observable: $k ($operator)")
        k => Observable(operator, atom, overlaps,
                        integrals, integral_map,
                        symmetries, selector;
                        double_counted=double_counted,
                        poisson_cache=poisson_cache,
                        verbosity=verbosity)
    end |> Dict{Symbol,Observable}

    AtomicEquations(atom, hfeqs, integrals, observables)
end

Base.Matrix(H::QuantumOperator, atom::Atom;
            overlaps::Vector{<:OrbitalOverlap}=OrbitalOverlap[],
            selector::Function=cfg -> outsidecoremodel(cfg, atom.potential),
            configurations = selector.(atom.configurations),
            kwargs...) =
    Matrix(H, configurations, overlaps)

function Base.diff(atom::Atom; H::QuantumOperator=atomic_hamiltonian(atom),
                   kwargs...)
    energy_expression = Matrix(H, atom; kwargs...)
    diff(atom, energy_expression; H=H, kwargs...)
end

# * Overlap matrix
function SCF.overlap_matrix(atom::Atom)
    R = radial_basis(atom)
    R'R
end


text:
module Julia3D

include("../julia_src/j3dapp.jl")

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.shortestpath" begin
    graph = plaingraph(10, [
        (1, 2), (2, 3), (1, 4), (4, 5), (3, 7),
        (7, 4), (7, 8), (8, 9), (9, 10)
    ])
    @test distance(graph, 1, 10) == 5
    @test distance(graph, 1, 1) === nothing
    @test distance(graph, 1, 6) === nothing
    @test length(reachablenodes(graph, 1)) == 8
    @test isreachable(graph, 4, 10)
    @test !isreachable(graph, 6, 1)
    @test eccentricity(graph, 8) == 3
    @test diameter(graph) == 5
    @test count(i->i==Inf, distancematrix(graph)) == 28
    @test shortestpathnodes(graph, 1, 10) == [1, 4, 7, 8, 9, 10]
    @test shortestpathnodes(graph, 5, 8) == [5, 4, 7, 8]
    @test length(longestshortestpathnodes(graph)) == 6

    dgraph = plaindigraph(10, [
        (1, 4), (2, 4), (3, 7), (4, 5), (4, 6),
        (4, 7), (6, 9), (7, 8), (7, 9), (7, 10)
    ])
    @test distance(dgraph, 1, 10) == 3
    @test reversedistance(dgraph, 1, 10) === nothing
    @test length(reachablenodes(dgraph, 4)) == 6
    @test !isreachable(dgraph, 3, 6)
    @test eccentricity(dgraph, 2) == 3
    @test diameter(dgraph) == 3
    @test count(i->i==3.0, distancematrix(dgraph)) == 6
    @test shortestpathnodes(dgraph, 3, 10) == [3, 7, 10]
    @test reverseshortestpathnodes(dgraph, 10, 3) == [10, 7, 3]
end # graph.shortestpath


text:
  # A RockSample lower bound estimation procedure for solvers that use
  # weighted particle belief representation (defined in types.jl)
  # TODO: think how to make this more generic
  #
  # Strategy: Compute a representative state by setting the state of
  # each rock to the one that occurs more frequently in the particle set.
  # Then compute the best sequence of actions for the resulting
  # state. Apply this sequence of actions to each particle and average
  # to get a lower bound.
  #
  # Possible improvement: If a rock is sampled while replaying the action
  # sequence, use dynamic programming to look forward in the action
  # sequence to determine if it would be a better idea to first sense the
  # rock instead. (sensing eliminates the bad rocks in the particle set)

mutable struct RockSampleParticleLB
    weight_sum_of_state::Vector{Float64}
    best_action::RockSampleAction

    function RockSampleParticleLB(pomdp::RockSample)
        this = new()
        this.weight_sum_of_state = Array{Float64}(pomdp.n_states)
        fill!(this.weight_sum_of_state, -Inf)
        this.best_action = RockSampleAction(-1)

        return this
    end
end

function init_bound(lb::RockSampleParticleLB,
                    pomdp::RockSample,
                    config::DESPOTConfig)
    # nothing to do for now
end

function lower_bound(lb::RockSampleParticleLB,
                     pomdp::RockSample,
                     particles::Vector{DESPOTParticle{RockSampleState}},
                     ub_actions::Vector{RockSampleAction},
                     config::DESPOTConfig)

    weight_sum::Float64 = 0.0
    state_seen::Dict{Int64,Int64} = Dict{Int64,Int64}()
    s_index::Int64 = -1

    # Since for this problem the cell that the rover is in is deterministic, picking pretty much
    # any particle state is ok
    if length(particles) > 0
        if isterminal(pomdp, particles[1].state)
            return 0.0 # lower bound value and best action
        end
    end

    # The expected value of sampling a rock, over all particles
    expected_sampling_value::Array{Float64} = fill(0.0, pomdp.n_rocks)
    seen_ptr::Int64 = 0

    # Compute the expected sampling value of each rock. Instead of factoring
    # the weight of each particle, we first record the weight of each state.
    # This is so that the inner loop that updates the expected value of each
    # rock runs once per state seen, instead of once per particle seen. If
    # there are lots of common states between particles, this gives a
    # significant speedup to the search because the lower bound is the
    # bottleneck.

    for p in particles
        if lb.weight_sum_of_state[state_index(pomdp,p.state)+1] == -Inf #Array
            lb.weight_sum_of_state[state_index(pomdp,p.state)+1] = p.weight
            state_seen[seen_ptr] = state_index(pomdp,p.state)
            seen_ptr += 1
        else
            lb.weight_sum_of_state[state_index(pomdp,p.state)+1] += p.weight;
        end
    end

    for i in 0:(seen_ptr-1)
        s_index = state_seen[i]
        weight_sum += lb.weight_sum_of_state[s_index+1]
        for j in 0:pomdp.n_rocks-1
            expected_sampling_value[j+1] += lb.weight_sum_of_state[s_index+1] * (rock_status(j, s_index) ? 10.0 : -10.0)
        end
    end

    # Reset for next use
    fill!(lb.weight_sum_of_state, -Inf)

    most_likely_rock_set = 0
    for i in 0:pomdp.n_rocks-1
        expected_sampling_value[i+1] /= weight_sum
        # Threshold the average to good or bad
        if expected_sampling_value[i+1] > -config.tiny
            most_likely_rock_set |= (1 << i)
        end
        if abs(expected_sampling_value[i+1]) < config.tiny
            expected_sampling_value[i+1] = 0.0
        end
    end

    # Since for this problem the cell that the rover is in is deterministic, picking pretty much
    # any particle state is ok
    s = RockSampleState(
        make_state_index(pomdp, cell_of(pomdp, particles[1].state), most_likely_rock_set))

    # Sequence of actions taken in the optimal policy
    optimal_policy = Vector{RockSampleAction}()
    ret::Float64 = 0.0
    reward::Float64 = 0.0
    prev_cell_coord = [0,0] # initial value - should cause error if not properly assigned
    next_s::RockSampleState = RockSampleState()
    r::Float64 = 0.0
    trans_distribution = create_transition_distribution(pomdp)
    rng = DESPOTRandomNumber(0) # dummy RNG

    while true
        a = ub_actions[state_index(pomdp,s)+1]
        trans_distribution.state = s
        trans_distribution.action = a
        next_s = POMDPs.rand(rng, trans_distribution, next_s)
        if isterminal(pomdp, next_s)
            prev_cell_coord[1] = pomdp.cell_to_coords[cell_of(pomdp, s)+1][1]
            prev_cell_coord[2] = pomdp.cell_to_coords[cell_of(pomdp, s)+1][2]
            ret = 10.0
            break
        end
        push!(optimal_policy,a)
        if length(optimal_policy) == config.search_depth
            prev_cell_coord[1] = pomdp.cell_to_coords[cell_of(pomdp, next_s)+1][1]
            prev_cell_coord[2] = pomdp.cell_to_coords[cell_of(pomdp, next_s)+1][2]
            ret = 0.0
            break
        end
        s = next_s #TODO: deepcopy?
    end

    lb.best_action = (length(optimal_policy) == 0) ? RockSampleAction(2) : optimal_policy[1]

    # Execute the sequence backwards to allow using the DP trick
    for i = length(optimal_policy):-1:1
        act = optimal_policy[i]
        ret *= pomdp.discount
        if action_index(pomdp,act) == 4
            rock = pomdp.rock_at_cell[cell_num(pomdp, prev_cell_coord[1], prev_cell_coord[2])+1]
            if rock != -1
                ret = expected_sampling_value[rock+1] + ret # expected sampling value is an array
            end
            continue
        end

        # Move in the opposite direction since we're going backwards
        if action_index(pomdp, act) == 0
            prev_cell_coord[1] += 1
        elseif action_index(pomdp, act) == 1
            prev_cell_coord[1] -= 1
        elseif action_index(pomdp, act) == 2
            prev_cell_coord[2] -= 1
        elseif action_index(pomdp, act) == 3
            prev_cell_coord[2] += 1
        else
            @assert(false)
        end
    end
    return ret
end