text:
import Base: *, ∈, isempty

export InverseLinearMap,
       an_element,
       constraints_list

"""
    InverseLinearMap{N, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}} <: AbstractAffineMap{N, S}

Given a linear transformation ``M``, this type represents the linear
transformation ``M⁻¹⋅X`` of a set ``X`` without actually computing ``M⁻¹``.

### Fields

- `M` -- invertible matrix
- `X` -- set

### Notes

This type is parametric in the elements of the inverse linear map, `NM`, which is
independent of the numeric type of the wrapped set (`N`). Typically `NM = N`,
but there may be exceptions, e.g., if `NM` is an interval that holds numbers
of type `N`, where `N` is a floating point number type such as `Float64`.

### Examples

For the examples we create a ``3×3`` matrix and a unit three-dimensional square.

```jldoctest ilp_constructor
julia> A = [1 2 3; 2 3 1; 3 1 2];

julia> X = BallInf([0, 0, 0], 1);

julia> ilm = InverseLinearMap(A, X)
InverseLinearMap{Int64, BallInf{Int64, Vector{Int64}}, Int64, Matrix{Int64}}([1 2 3; 2 3 1; 3 1 2], BallInf{Int64, Vector{Int64}}([0, 0, 0], 1))
```

Applying an inverse linear map to a `InverseLinearMap` object combines the two maps into
a single `InverseLinearMap` instance.

```jldoctest ilp_constructor
julia> B = transpose(A); ilm2 = InverseLinearMap(B, ilm)
InverseLinearMap{Int64, BallInf{Int64, Vector{Int64}}, Int64, Matrix{Int64}}([14 11 11; 11 14 11; 11 11 14], BallInf{Int64, Vector{Int64}}([0, 0, 0], 1))

julia> ilm2.M == A*B
true
```

The application of an `InverseLinearMap` to a `ZeroSet` or an `EmptySet` is
simplified automatically.

```jldoctest ilp_constructor
julia> InverseLinearMap(A, ZeroSet{Int}(3))
ZeroSet{Int64}(3)
```
"""
struct InverseLinearMap{N, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}} <: AbstractAffineMap{N, S}
    M::MAT
    X::S

    # default constructor with dimension match check
    function InverseLinearMap(M::MAT, X::S; check_invertibility::Bool=false) where {N, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}}
        @assert dim(X) == size(M, 1) "a linear map of size $(size(M)) cannot " *
                                     "be applied to a set of dimension $(dim(X))"
        if check_invertibility
            @assert isinvertible(M) "the linear map is not invertible"
        end
        return new{N, S, NM, MAT}(M, X)
    end
end

# convenience constructor from a UniformScaling
function InverseLinearMap(M::UniformScaling{N}, X::LazySet; check_invertibility::Bool=false) where {N}
    return InverseLinearMap(M.λ, X, check_invertibility=check_invertibility)
end

# convenience constructor from a scalar
function InverseLinearMap(α::N, X::LazySet; check_invertibility::Bool=false) where {N<:Real}
    if check_invertibility
        @assert !iszero(α) "the linear map is not invertible"
    end

    if α == one(N)
        return X
    end

    D = Diagonal(fill(α, dim(X)))
    return InverseLinearMap(D, X, check_invertibility=false)
end

# combine two linear maps into a single linear map
function InverseLinearMap(M::AbstractMatrix, ilm::InverseLinearMap)
    return InverseLinearMap(ilm.M * M, ilm.X)
end

# ZeroSet is "almost absorbing" for InverseLinearMap (only the dimension changes)
function InverseLinearMap(M::AbstractMatrix, Z::ZeroSet)
    @assert dim(Z) == size(M, 2) "a linear map of size $(size(M)) cannot " *
            "be applied to a set of dimension $(dim(Z))"
    return Z
end

# EmptySet is absorbing for LinearMap
function InverseLinearMap(M::AbstractMatrix, ∅::EmptySet)
    @assert dim(∅) == size(M, 2) "a linear map of size $(size(M)) cannot " *
            "be applied to a set of dimension $(dim(∅))"
    return ∅
end

# --- AbstractAffineMap interface functions ---

function matrix(ilm::InverseLinearMap)
    return inv(ilm.M)
end

function vector(ilm::InverseLinearMap{N}) where {N}
    return spzeros(N, dim(ilm))
end

function set(ilm::InverseLinearMap)
    return ilm.X
end

"""
    dim(ilm::InverseLinearMap)

Return the dimension of an inverse linear map.

### Input

- `ilm` -- inverse linear map

### Output

The ambient dimension of the inverse linear map.
"""
function dim(ilm::InverseLinearMap)
    return size(ilm.M, 1)
end

"""
    σ(d::AbstractVector, ilm::InverseLinearMap)

Return the support vector of the inverse linear map.

### Input

- `d`   -- direction
- `ilm` -- inverse linear map

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``L = M^{-1}⋅X``, where ``M`` is a matrix and ``X`` is a set, since
(M^T)^{-1}=(M^{-1})^T, it follows that ``σ(d, L) = M^{-1}⋅σ((M^T)^{-1} d, X)``
for any direction ``d``.
"""
function σ(d::AbstractVector, ilm::InverseLinearMap)
    y = transpose(ilm.M) \ d
    return ilm.M \ σ(y, ilm.X)
end

"""
    ρ(d::AbstractVector, ilm::InverseLinearMap)

Return the support function of the inverse linear map.

### Input

- `d`      -- direction
- `ilm`    -- inverse linear map

### Output

The support function in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``L = M^{-1}⋅X``, where ``M`` is a matrix and ``X`` is a set, it follows
that ``ρ(d, L) = ρ((M^T)^{-1} d, X)`` for any direction ``d``.
"""
function ρ(d::AbstractVector, ilm::InverseLinearMap)
    y = transpose(ilm.M) \ d
    return ρ(y, ilm.X)
end

"""
    ∈(x::AbstractVector, ilm::InverseLinearMap)

Check whether a given point is contained in the inverse linear map of a set.

### Input

- `x`   -- point/vector
- `ilm` -- inverse linear map of a set

### Output

`true` iff ``x ∈ ilm``.

### Algorithm

This implementation does not explicitly invert the matrix since it uses the
property ``x ∈ M^{-1}⋅X`` iff ``M⋅x ∈ X``.

### Examples

```jldoctest
julia> ilm = LinearMap([0.5 0.0; 0.0 -0.5], BallInf([0., 0.], 1.));

julia> [1.0, 1.0] ∈ ilm
false

julia> [0.1, 0.1] ∈ ilm
true
```
"""
function ∈(x::AbstractVector, ilm::InverseLinearMap)
    y = ilm.M * x
    return y ∈ ilm.X
end

"""
    an_element(ilm::InverseLinearMap)

Return some element of an inverse linear map.

### Input

- `ilm` -- inverse linear map

### Output

An element in the inverse linear map.
It relies on the `an_element` function of the wrapped set.
"""
function an_element(lm::InverseLinearMap)
    return lm.M \ an_element(lm.X)
end

"""
    vertices_list(ilm::InverseLinearMap; prune::Bool=true)

Return the list of vertices of a (polyhedral) inverse linear map.

### Input

- `ilm`   -- inverse linear map
- `prune` -- (optional, default: `true`) if `true` removes redundant vertices

### Output

A list of vertices.

### Algorithm

We assume that the underlying set `X` is polyhedral.
Then the result is just the inverse linear map applied to the vertices of `X`.
"""
function vertices_list(ilm::InverseLinearMap; prune::Bool=true)
    # collect vertices list of wrapped set
    vlist_X = vertices_list(ilm.X)

    # create resulting vertices list
    vlist = Vector{eltype(vlist_X)}(undef, length(vlist_X))
    @inbounds for (i, vi) in enumerate(vlist_X)
        vlist[i] = ilm.M \ vi
    end

    return prune ? convex_hull(vlist) : vlist
end

"""
    constraints_list(ilm::InverseLinearMap)

Return the list of constraints of a (polyhedral) inverse linear map.

### Input

- `ilm` -- inverse linear map

### Output

The list of constraints of the inverse linear map.

### Notes

We assume that the underlying set `X` is polyhedral, i.e., offers a method
`constraints_list(X)`.

### Algorithm

We fall back to a concrete set representation and apply `linear_map_inverse`.
"""
function constraints_list(ilm::InverseLinearMap)
    return constraints_list(linear_map_inverse(ilm.M, ilm.X))
end

"""
    linear_map(M::AbstractMatrix, ilm::InverseLinearMap)

Return the linear map of a lazy inverse linear map.

### Input

- `M`   -- matrix
- `ilm` -- inverse linear map

### Output

The set representing the linear map of the lazy inverse linear map of a set.

### Notes

This implementation is inefficient because it will compute the
concrete inverse of ``M``, which is what `InverseLinearMap` is supposed to avoid.
"""
function linear_map(M::AbstractMatrix, ilm::InverseLinearMap)
    return linear_map(M * inv(ilm.M), ilm.X)
end

function concretize(ilm::InverseLinearMap)
    return linear_map(inv(ilm.M), concretize(ilm.X))
end


text:
const HubbardStateType = Vector{Int}

struct HubbardHilbertSpace{LatticeType<:AbstractSiteNetwork,IndexType<:AbstractUniqueVector{HubbardStateType}} <: HilbertSpace{HubbardStateType}
    lattice::LatticeType
    indexer::IndexType
end

HubbardHilbertSpace(lattice) = HubbardHilbertSpace(lattice, UniqueVector{HubbardStateType}())

# -- 0
# u- 1
# -d 2
# ud 3

# only the first two bits are allowed to be set
is_valid_hubbard_site_state(site_state::Integer) = (site_state & ~(zero(site_state) ⊻ 3)) == 0

function get_σz(::HubbardHilbertSpace, site_state::Integer)
    @assert is_valid_hubbard_site_state(site_state)
    # maps [0, 1, 2, 3] to [0, 1, -1, 0]
    return (site_state & 1) - ((site_state & 2) >> 1)
end

#"""
#Returns charge on site
#"""
function get_charge(::HubbardHilbertSpace, site_state::Integer)
    @assert is_valid_hubbard_site_state(site_state)
    # maps [0, 1, 2, 3] to [0, 1, 1, 2]
    return (site_state & 1) + ((site_state & 2) >> 1)
end

function get_total_charge(hs::HubbardHilbertSpace, state::Vector)
    rv = 0
    for x in state
        rv += get_charge(hs, x)
    end
    return rv
end

get_total_charge(hs::HubbardHilbertSpace, stateidx::Int) = get_total_charge(hs, hs.indexer[stateidx])

# In our second quantization, we use the convention that all up creation
# operators come before (or after) the dn creation operators.

# Returns the phase tracker for both spin up and spin dn in the
# respective bits of the return value, under the condition that a
# particle is hopping between sites x1 and x2.
function phase_tracker(state::Vector, x1::Integer, x2::Integer)
    # re-order such that x1 <= x2
    if x1 > x2
        x1, x2 = x2, x1
    end
    @assert x1 <= x2

    rv = 0
    for i in x1+1:x2-1
        rv ⊻= state[i]
    end
    return rv
end

function apply_total_spin_operator(f, hs::HubbardHilbertSpace, s_j::Integer)
    state = hs.indexer[s_j]
    diagonal = 0.0

    for x in 1:length(hs.lattice)
        for x_r in 1:length(hs.lattice)
            # 0.5 * (S^+_i S^-_j + S^-_i S^+_j)
            if x == x_r
                if state[x] == 1 || state[x] == 2
                    diagonal += 0.5
                end
            elseif (state[x] == 1 && state[x_r] == 2) || (state[x] == 2 && state[x_r] == 1)
                other = copy(state)
                other[x], other[x_r] = other[x_r], other[x]
                s_i = findfirst!(isequal(other), hs.indexer)

                # figure out whether we need to pick up a phase for the fermions
                pt = phase_tracker(state, x, x_r)
                # the following line maps [0, 1, 2, 3] to [1, -1, -1, 1]
                mult = 1 - ((pt & 2) ⊻ ((pt & 1) << 1))
                @assert mult == 1 || mult == -1

                # the minus sign comes from working through the site hops carefully.
                f(s_i, -0.5 * mult)
            end
            # S^z_i S^z_j
            diagonal += 0.25 * get_σz(hs, state[x]) * get_σz(hs, state[x_r])
        end
    end

    f(s_j, diagonal)
    nothing
end

function apply_total_pseudospin_operator(f, hs::HubbardHilbertSpace, s_j::Integer)
    @assert isbipartite(hs.lattice)

    state = hs.indexer[s_j]
    diagonal = 0.0

    for x in 1:length(hs.lattice)
        sublattice_index_x = sublattice_index(hs.lattice, x)
        for x_r in 1:length(hs.lattice)
            # 0.5 * (S^+_i S^-_j + S^-_i S^+_j) under duality
            if (x == x_r)
                if state[x] == 0 || state[x] == 3
                    diagonal += 0.5
                end
            elseif (state[x] == 0 && state[x_r] == 3) || (state[x] == 3 && state[x_r] == 0)
                other = copy(state)
                other[x], other[x_r] = other[x_r], other[x]
                s_i = findfirst!(isequal(other), hs.indexer)

                # figure out whether we need to pick up a phase for the fermions
                pt = phase_tracker(state, x, x_r)
                # the following line maps [0, 1, 2, 3] to [1, -1, -1, 1]
                mult = 1 - ((pt & 2) ⊻ ((pt & 1) << 1))
                @assert mult == 1 || mult == -1

                sublattice_parity = sublattice_index_x ⊻ sublattice_index(hs.lattice, x_r)
                @assert sublattice_parity & 1 == sublattice_parity
                sublattice_factor = 1 - 2 * sublattice_parity
                f(s_i, 0.5 * sublattice_factor * mult)
            end
            # S^z_i S^z_j under duality
            diagonal += 0.25 * (get_charge(hs, state[x]) - 1) * (get_charge(hs, state[x_r]) - 1)
        end
    end

    f(s_j, diagonal)
    nothing
end

struct HubbardParameters
    t::Float64
    U::Float64
    V::Float64
    W::Float64
    J_xy::Float64
    J_z::Float64
    μ::Float64
    ϵ_total_spin::Float64
    ϵ_total_pseudospin::Float64
    t2::Float64

    function HubbardParameters(;
                               t::Real=0.0,
                               U::Real=0.0,
                               V::Real=0.0,
                               W::Real=0.0,
                               J::Real=0.0,
                               J_xy::Real=0.0,
                               J_z::Real=0.0,
                               μ::Real=0.0,
                               ϵ_total_spin::Real=0.0,
                               ϵ_total_pseudospin::Real=0.0,
                               t2::Real=0.0)
        if J != 0
            J_xy == J_z == 0 || throw(ArgumentError("If J is provided, J_xy and J_z must not be."))
            J_xy = J_z = J
        end

        new(t, U, V, W, J_xy, J_z, μ, ϵ_total_spin, ϵ_total_pseudospin, t2)
    end
end

hubbard_hamiltonian(; kwargs...) = hubbard_hamiltonian(HubbardParameters(; kwargs...))

function hubbard_neighbor_terms(f, hs, s_j, x::Int, x_r::Int, η::Rational{Int}, t_up, t_dn, J_xy=0, J_z=0, V=0, W=0)
    state = hs.indexer[s_j]
    diagonal = 0.0

    # figure out whether we need to pick up a phase
    pt = phase_tracker(state, x, x_r)
    x_r_up_phase = 1 - ((pt & 1) << 1)
    x_r_dn_phase = 1 - (pt & 2)
    e_iθ = exp_2πiη(η)
    @assert x_r_up_phase == 1 || x_r_up_phase == -1
    @assert x_r_dn_phase == 1 || x_r_dn_phase == -1

    # up spin hopping
    b_up = 1 # const
    if state[x_r] & b_up != 0
        # backward
        if state[x] & b_up == 0
            other = copy(state)
            other[x_r] -= b_up
            other[x] += b_up
            s_i = findfirst!(isequal(other), hs.indexer)
            f(s_i, -t_up * x_r_up_phase * e_iθ)
        end
    end
    if state[x] & b_up != 0
        # forward
        if state[x_r] & b_up == 0
            other = copy(state)
            other[x] -= b_up
            other[x_r] += b_up
            s_i = findfirst!(isequal(other), hs.indexer)
            f(s_i, -t_up * x_r_up_phase * conj(e_iθ))
        end
    end

    # dn spin hopping
    b_dn = 2 # const
    if state[x_r] & b_dn != 0
        # backward
        if state[x] & b_dn == 0
            other = copy(state)
            other[x_r] -= b_dn
            other[x] += b_dn
            s_i = findfirst!(isequal(other), hs.indexer)
            f(s_i, -t_dn * x_r_dn_phase * e_iθ)
        end
    end
    if state[x] & b_dn != 0
        # forward
        if state[x_r] & b_dn == 0
            other = copy(state)
            other[x] -= b_dn
            other[x_r] += b_dn
            s_i = findfirst!(isequal(other), hs.indexer)
            f(s_i, -t_dn * x_r_dn_phase * conj(e_iθ))
        end
    end

    # exchange J
    if J_xy != 0
        # 0.5 * (S^+_i S^-_j + S^-_i S^+_j)
        if ((state[x] == 1 && state[x_r] == 2)
            || (state[x] == 2 && state[x_r] == 1))
            other = copy(state)
            other[x], other[x_r] = other[x_r], other[x]
            s_i = findfirst!(isequal(other), hs.indexer)
            # the minus sign comes from working through the site hops carefully.
            f(s_i, -0.5 * J_xy * x_r_up_phase * x_r_dn_phase)
        end
    end
    if J_z != 0
        # S^z_i S^z_j
        diagonal += 0.25 * J_z * get_σz(hs, state[x]) * get_σz(hs, state[x_r])
    end

    # Neighbor repulsion "V"
    if V != 0
        c = get_charge(hs, state[x]) * get_charge(hs, state[x_r])
        diagonal += V * c
    end

    # Neighbor doublon repulsion "W"
    if W != 0
        if (state[x] == 3 && state[x_r] == 3)
            diagonal += W
        end
    end

    f(s_j, diagonal)
    nothing
end

function hubbard_hamiltonian(p::HubbardParameters)
    return function apply_hamiltonian(f, hs::HubbardHilbertSpace, s_j::Integer)
        state = hs.indexer[s_j]
        diagonal = 0.0

        neighborsη(hs.lattice) do x::Int, x_r::Int, η::Rational{Int}
            hubbard_neighbor_terms(f, hs, s_j, x, x_r, η, p.t, p.t, p.J_xy, p.J_z, p.V, p.W)
        end

        if p.t2 != 0
            neighborsη(hs.lattice, Val{2}) do x::Int, x_r::Int, η::Rational{Int}
                hubbard_neighbor_terms(f, hs, s_j, x, x_r, η, p.t2, p.t2)
            end
        end

        # Hubbard U
        if p.U != 0
            doubly_occupied_sites = 0
            for x in 1:length(hs.lattice)
                doubly_occupied_sites += ifelse(state[x] == 3, 1, 0)
            end
            diagonal += p.U * doubly_occupied_sites
        end

        # chemical potential μ
        if p.μ != 0
            total_charge = 0
            for x in 1:length(hs.lattice)
                total_charge += get_charge(hs, state[x])
            end
            diagonal += -p.μ * total_charge
        end

        # total spin term
        if p.ϵ_total_spin != 0
            apply_total_spin_operator(edapply(f, p.ϵ_total_spin), hs, s_j)
        end

        # total pseudospin term
        if p.ϵ_total_pseudospin != 0
            apply_total_pseudospin_operator(edapply(f, p.ϵ_total_pseudospin), hs, s_j)
        end

        f(s_j, diagonal)
        nothing
    end
end

function seed_state!(hs::HubbardHilbertSpace; N_up::Int, N_dn::Int)
    if !(0 <= N_up <= length(hs.lattice))
        throw(ArgumentError("Invalid N_up provided for size $(length(hs.lattice)) lattice: $(N_up)"))
    end
    if !(0 <= N_dn <= length(hs.lattice))
        throw(ArgumentError("Invalid N_dn provided for size $(length(hs.lattice)) lattice: $(N_dn)"))
    end
    state = zeros(Int, length(hs.lattice))
    for i in 1:N_up
        state[i] ⊻= 1
    end
    for i in 1:N_dn
        state[i] ⊻= 2
    end
    findfirst!(isequal(state), hs.indexer)
    return hs
end

# Deprecated 2018-10-16
@deprecate seed_state!(hs::HubbardHilbertSpace, N_up::Integer, N_dn::Integer) seed_state!(hs, N_up=N_up, N_dn=N_dn)

function permutation_parity(perm::Vector{Int})
    parity = 0
    nparticles = length(perm)
    inspected = falses(nparticles)
    for i in 1:nparticles
        inspected[i] && continue
        j = i
        while true
            j = perm[j]
            parity ⊻= 1
            @assert !inspected[j]
            inspected[j] = true
            j == i && break
        end
        parity ⊻= 1
    end
    return parity
end

function translateη(hs::HubbardHilbertSpace, ltrc::LatticeTranslationCache, j::Integer)
    state = hs.indexer[j]
    sz = length(state)
    newstate = zero(state)
    phase = 0//1

    parity = 0
    for bit in 1:2
        # Construct the sequence of creation operators
        cdagger = Int[]
        for i in 1:sz
            if state[i] & bit != 0
                new_site_index, η = translateη(ltrc, i)
                @assert 0 < new_site_index <= sz
                newstate[new_site_index] ⊻= bit
                push!(cdagger, new_site_index)
                phase -= η
            end
        end

        # Determine the parity of the permutation that orders them
        perm = sortperm(cdagger, alg=TimSort)
        parity ⊻= permutation_parity(perm)
    end

    @assert parity == 0 || parity == 1
    if parity != 0
        phase += 1//2
    end

    return findfirst!(isequal(newstate), hs.indexer), phase
end

function site_spinflip(site_state)
    let x = (site_state + 1) & 2
        return site_state ⊻ (x | (x >> 1))
    end
end

function spinflipη(hs::HubbardHilbertSpace, j::Integer)
    state = hs.indexer[j]
    i = findfirst!(isequal(map(site_spinflip, state)), hs.indexer)

    # We are implementing the transformation c_↑ ↦ c_↓ and vice versa.  Since
    # our second-quantization convention assumes all spin-up operators come
    # before (or, equivalently, after), spin-down operators, we need to account
    # for this permutation of operators.  The permutation is odd if and only if
    # there is an odd number of up spins and an odd number of down spins.  The
    # following line picks up a sign when this is the case.
    η = (reduce(xor, state) == 3) ? 1//2 : 0//1

    return i, η
end

function particleholeη(hs::HubbardHilbertSpace, j::Integer)
    state = hs.indexer[j]
    i = findfirst!(isequal(map(x -> x ⊻ 3, state)), hs.indexer)
    return i, 0//1
end

ext: jl


text:
# These Symata functions (symbols) translate between Symata and Julia data and expressions

### JVar

@sjdoc JVar """
    JVar(x)

return the Julia value of the Symbol that `x` evaluates to.

For example, if `a = 1` in Julia and `b = a` in Symata, then `JVar(b)` evaluates to `1`.
"""

@sjseealso_group(Jxpr, JVar)
apprules(mx::Mxpr{:JVar}) = Core.eval(Main,symname(mx[1]))

### SetJ

@mkapprule SetJ nargs => 2  """
    SetJ(x,val)

sets the Julia symbol `x` to `val`.

Variables and functions in Symata are separate from those in Julia, ie, their table of bindings to symbols are separate.
"""

exportj(lhs::SymString) = exportj(lhs,doeval(lhs))
## FIXME: find out what we want on the next line and do it properly.
exportj(lhs::SymString, rhs) = Main.eval(Expr(:(=),symname(Symbol(lhs)),rhs))
exportj(lhs::SymString, rhs::Symbol) = exportj(lhs,QuoteNode(rhs))
@doap SetJ(lhs::SymString,rhs) = exportj(lhs,rhs)

### ExportJ

@mkapprule ExportJ """
    ExportJ(x,y,...)

sets the Julia symbols `x`,`y`, etc. to the Symata-value of `x`,`y`, etc.
`ExportJ(x)` is equivalent to `SetJ(x,x)`.
"""

@doap ExportJ(vars::SymString...) = foreach(exportj, vars)

### Jxpr

@sjdoc Jxpr """
    Jxpr

allows embedding Julia expressions. A Jxpr is entered like this `J( expr )`.

`expr` is interpreted as a Julia expression and it is wrapped expression with head `Jxpr`, which is then evaluated when
`Jxpr` is evaluated. You never see the head `Jxpr`. For example,
 `m = J( collect(1:10) )`  creates a Julia array and binds it to the Symata symbol `m`.
"""

@sjexamp( Jxpr,
         "This creates a Julia Array{Int,1} and \"binds\" it to the Symata symbol m.",
         ("m = J( collect(1:3) )",
          "3-element Array{Int64,1}:\n 1\n 2\n 3"))

@sjexamp( Jxpr,
         "Call a Julia function",
         ("tf = J( time )",""),
         ("tf()","1.424287593897437e9"))

# quote, i.e. :( expr ) is parsed as a Julia expression and is wrapped as  < -- we will remove this
# Mxpr with head Jxpr. It is evaluated here.
@mkapprule Jxpr :nodefault => true

@doap function Jxpr(args...)
    local res
    for x in args
        res = eval(x)
    end
    res
end

@sjdoc J """
    J(expr1,expr2,...)

Evaluate Julia expressions `expr1`, `expr2`, ... in Julia and return the final result.
"""

@sjseealso_group(J, Compile, SymataCall)

#### Unpack

@sjdoc Unpack """
    Unpack(a)

unpack a Julia typed `AbstractArray` into an Symata `List` expression.

Only `1`-d is supported. If `a` is a Julia `Dict`, then a `List` of `Lists` of
key,value pairs is returned.
"""

@sjexamp( Unpack,
         "This creates a List of three random Float64's.",
         ("Unpack( J(rand(3)) )", "[0.5548766917324894,0.034964001133465095,0.9122052258982192]"))

@mkapprule Unpack :nargs => 1

@doap Unpack(a::AbstractArray) = unpacktoList(a)

"""
    unpacktoList(obj::AbstractArray{T,1})

Converts a Julia array to a Symata list of type `ListT`. Only one-dimensional
arrays are supported.
"""
function unpacktoList(obj)
    mx = mxpr(:List, do_unpack(obj))
    setfixed(mx)
    setcanon(mx)
    return mx
end
do_unpack(obj) = copyto!(newargs(length(obj)),obj)

## We copy this and change it a bit to tolistoflists in mxpr_utils
function do_unpack(dict::Dict)
    args = newargs(length(dict))
    i = 0
    for (k,v) in dict
        i += 1
        args[i] = mxpr(:List,k,v)
    end
    return args
end

#### Pack

@sjdoc Pack """
    Pack(expr)

pack the arguments of the Symata expression `expr` into a Julia array,
whose element type is the minimum required to holds the arguments of `expr`.
"""

@sjexamp( Pack,
         "This returns a Julia array of element type Int [1,2,3].",
         ("ClearAll(f)",""),
         ("Pack(f(1,2,3))","3-element Array{Int64,1}: [1,2,3]"))

@sjseealso_group(Pack,Unpack)

"""
    typejoin_array(a::Array)

return the typejoin of all elements in `a`.
I saw this in Base somewhere. Probably more efficient there.
"""
function typejoin_array(a)
    length(a) == 0 && return Any
    T = typeof(a[1])
@inbounds  for i in 2:length(a)
        T = typejoin(T, typeof(a[i]))
    end
    T
end

# 1-d unpack
function apprules(mx::Mxpr{:Pack})
    a = margs(margs(mx)[1])
    T = typejoin_array(a)
    do_pack(T,a)
end

do_pack(T,sjobj) = copyto!(Array{T}(undef, length(sjobj)), sjobj)

#### Translate Symata to Julia

##
# Wrap Expr to prevent Symata from evaluating it.

abstract type AbstractMtoE end

mutable struct MtoECompile <: AbstractMtoE
end

mutable struct MtoEPlain <: AbstractMtoE
end

mxpr_to_expr(x, aux::AbstractMtoE) = x

# Don't really want to do this with everything !
function mxpr_to_expr(s::Symbol, aux::AbstractMtoE)
    s
#    Symbol(lowercase(string(s)))
end

function mxpr_to_expr(s::Symbol,aux::MtoEPlain)
    Symbol(lowercase(string(s)))  # only lowercase some
end

const MTOJSYM_COMPILE = Dict(
    :Times => :mmul,
    :Plus => :mplus,
    :Power => :mpow,
    :Abs => :mabs,
    :Exp => :Exp  # Is this necessary ?
)

mtojsym_compile(s::Symbol) =  get(MTOJSYM_COMPILE, s, mtojsym(s))

const MTOJSYM_PLAIN = Dict(
                             :Times => :(*),
                             :Plus => :(+),
                             :Power => :(^),
                             :Abs => :abs,
                             :E => :e
                             )

function mtojsym_plain(s::Symbol)
    Symbol(lowercase(string(get(MTOJSYM_PLAIN, s, mtojsym(s)))))
end

function mxpr_to_expr(mx::Mxpr,aux::MtoECompile)
    head = mtojsym_compile(mhead(mx))
    mxpr_to_expr_body(head,mx,aux)
end

function mxpr_to_expr(mx::Mxpr{:Power}, aux::MtoECompile)
    if mx[1] == :E
        #        return Expr(:call, :Exp, margs(mx)[2:end]...)
        return Expr(:call, :Exp, [mxpr_to_expr(x,aux) for x in margs(mx)[2:end]]...)
    else
#        return Expr(:call, :mpow, margs(mx)[1:end]...)
        return Expr(:call, :mpow, [mxpr_to_expr(x,aux) for x in margs(mx)]...)
    end
end

# convert List to Julia array
function mxpr_to_expr(mx::Mxpr{:List},aux::MtoECompile)
    Expr(:vect, [mxpr_to_expr(x,aux) for x in margs(mx)]...)
end

function mxpr_to_expr(mx::Mxpr,aux::MtoEPlain)
    head = mtojsym_plain(mhead(mx))
    mxpr_to_expr_body(head,mx,aux)
end

function mxpr_to_expr_body(head,mx::Mxpr,aux::AbstractMtoE)
    isempty(mx) && return :( $(head)() )
    a = margs(mx)
    a1 = mxpr_to_expr(a[1],aux)
    ex = :( $(head)($(a1)) )
    length(mx) == 1 && return ex
    for i in 2:length(a)
        push!(ex.args, mxpr_to_expr(a[i],aux))
    end
    ex
end

#### Compile

@mkapprule Compile """
    f = Compile(expr)

convert `expr` to a compiled function.

    f = Compile( [x,y,...], expr)

create an anonymous Julia function `(x,y,...) -> expr` from Symata expression `expr` and bind
the result to `f`.

```
symata> f = Compile(x^2 + y^2)
symata> f(3,4)
        25
```

Compile works by creating an anonymous Julia function from the Symata expression `expr`. This anonymous
function is just-in-time compiled for the arguments that it is called with.
If no arguments are specified, then the function signature is the sorted list of free symbols found in `expr`. Free symbols are all symbols
that are not bound in the Julia. More precisely, in the Symata module scope. This
excludes `e`, for instance. Alternatively, the arguments may be specified, for example `Compile( [x,y], expr)`.

By default `expr` is not evaluated. This works for literal expressions, such as `x^2`. If the expression must be evaluated, wrap it
in `Evaluate`. For example

```
symata 1> expr = x^2 + y^3
Out(1) = x^2 + y^3

symata 2> f = Compile(Evaluate(expr))
Out(2) = (::#1) (generic function with 1 method)

symata 3> f(2,3)
Out(3) = 31
```
"""
@doap function Compile(a::Mxpr{:List}, body)
    aux = MtoECompile()
    jexpr = Expr(:function, Expr(:tuple, [mxpr_to_expr(x, aux) for x in margs(a)]...) , mxpr_to_expr(body, aux))
    Core.eval(Main, jexpr)
end

# This is not reliable. It is not worth the complexity of fixing it, I think.
# We require an explicit list of variables, instead.
# We could look at Attributes, Protected, etc. and just use plain variables.
# @doap function Compile(body)
#     aux = MtoECompile()
#     jbody = mxpr_to_expr(body,aux)
#     Core.eval(Main, Expr(:function, Expr(:tuple, freesyms(jbody)...), jbody))
# end

## NOTE: freesyms is unrelated to getfreesyms and setfreesyms
"""
    freesyms(ex::Expr)

return a list of Symata-unbound symbols at any depth in expression `ex`, by
walking through the expression. This may not be useful anymore.
"""
function freesyms(x)
    syms = Dict()
    freesyms!(x, syms)
    sort(collect(keys(syms)))
end

function freesyms!(ex::Expr, syms::Dict)
    a = ( ex.head == :call ?  view(ex.args,2:length(ex.args)) : ex.args )
    foreach(x -> freesyms!(x,syms), a)
end

# We have been using isdefined for a long time,
# but this is incorrect. We want isbound, which
# checks for Symata binding.
# Neither is correct. For instance, E is not Symata bound.
# E is bound in Julia Main. (we did it). So, the fact
# that this was worked for several examples, is fortuitous.
function freesyms!(s::Symbol, syms)
#    if ! isdefined(Symata,s)
    if ! isbound(s)
        syms[s] = 1
    end
end

freesyms!(x, syms) = nothing

####  ToJuliaExpression

@mkapprule ToJuliaExpression

@doap function ToJuliaExpression(x)
    mxpr(:JuliaExpression, mxpr_to_expr(x, MtoEPlain()))
end

"""
    @symExpr(expr)

Symata-evaluates `expr`, translates the result to a Julia expression, and inserts it into the surrounding code. This works
just like any Julia macro for inserting code, except that the code is generated by evaluating a Symata expression.

The following evaluates Symata code, translates the result to Julia and uses it as the body of a Julia function.
```
julia> f1(z) = @symExpr Together(PolyLog(-1,z) + (1-z))
```
"""
macro symExpr(expr)
    mx = symataevaluate(expr, EvaluateJuliaSyntaxSimple())
    mx = doeval(mx)
    local ex
    if isa(mx,Mxpr{:JuliaExpression})
       ex = mx[1]
    else
       ex = mxpr_to_expr(mx, MtoECompile())
    end
    :(($(esc(ex))))  ## We need to escape this so that the code is evaluated in the correct module
end

### ToJuliaString

@sjdoc ToJuliaString """
    ToJuliaString(expr)

translate `expr` to Julia, returning a string.
If the option `NoSymata => False` is given, then translated code requires `using Symata`.
Currently, only arithemetic, trigonometric functions, and a few other things are translated.
"""
@mkapprule ToJuliaString :nodefault => true  :options => Dict( :NoSymata => true )

protect(:NoSymata)

function do_ToJuliaString(mx, x; NoSymata=true)
    aux = NoSymata ? MtoEPlain() : MtoECompile()
    string(mxpr_to_expr(x, aux))
end

@mkapprule SymataCall  """
    SymataCall(var, body)
    SymataCall(varlist, body)

return a Julia function that substitues `var` or variable in `varlist` into `body` and
submits the result to the Symata evaluation sequence.

The variables are localized before writing the Julia function.
"""
@doap SymataCall(var::Symbol, ex::Mxpr) = wrap_symata(doeval(ex),var)
@doap SymataCall(vars::Mxpr{:List}, ex::Mxpr) = wrap_symata(doeval(ex),vars...)

"""
    f = wrap_symata(expr::Mxpr,var)

return a Julia function that Symata-evaluates an expression when called.

The function corresponds to `f(var) = expr`.
This wraps `expr` in a single-argument Julia function.

First creates a lexical scope, that is a new gensym symbol `nvar` is generated and is substituted
everywhere in `expr` for `var`. Calling `f(x)`  Symata-binds `x` to `nvar` and `expr` is
Symata-evaluated and the result is returned.

`wrap_symata` is used by `NIntegrate` and `Plot`.
"""
function wrap_symata(var0,expr0::Mxpr)   # Why are arguments reversed here ? Is this ever called ?
    (var,expr) = localize_variable(var0,expr0)
#    deepunsetfixed(expr)
    function (x)
        setsymval(var,x)
        # deepunsetfixed(expr)
        # this seems wasteful. When called from NIntegrate, it is not needed. Don't know why.
        # because we needed HoldAll
        doeval(expr)
    end
end

function wrap_symata(expr0::Mxpr,vars0...)
    vars = Array{Symbol}(undef, length(vars0))
    expr = expr0
    for i in 1:length(vars0)
        (vars[i],expr) = localize_variable(vars0[i],expr)
    end
    deepunsetfixed(expr)
    function (args...)
        for i in 1:length(args)
            setsymval(vars[i],args[i])
        end
        deepunsetfixed(expr)
        doeval(expr)
    end
end

function wrap_symata(expr0::SJSym,var0)
    var0 == expr0 && return (x) -> x
    sym = get_localized_symbol(var0)
    return (sym,expr0)   # this does not appear to be useful
end

#### CodeNative

# This does not yet work. Symata evaluates the symbol `f` to
# the function object, so we lose its name. We probably need
# a robust way to associate the symbol with the function.
# ... Julia numbers anonymous functions and names others.
# we can get this information somehow
# ok, it is string(func::Function). Then we can keep a
# hash table of Symata symbols bound to functions.
# ... a PITA.
# Better to define a Julia type symfunction or something
# that stores the function and the name.
# we can make it callable and it passes the call to its function
# can also do this at the Symata level.

# FIXME. implement functions according to the description above.

# @mkapprule CodeNative
# @doap function CodeNative(f::Symbol, arg)
#     thetypes = (margs(args)...)
#     code_native(eval(f), thetypes)
# end


text:
# ---
# title: 1269. Number of Ways to Stay in the Same Place After Some Steps
# id: problem1269
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-ways-to-stay-in-the-same-place-after-some-steps/description/>
# hidden: true
# ---
#
# You have a pointer at index `0` in an array of size `arrLen`. At each step,
# you can move 1 position to the left, 1 position to the right in the array or
# stay in the same place  (The pointer should not be placed outside the array at
# any time).
#
# Given two integers `steps` and `arrLen`, return the number of ways such that
# your pointer still at index `0` after **exactly**` steps` steps.
#
# Since the answer may be too large, return it **modulo**  `10^9 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: steps = 3, arrLen = 2
#     Output: 4
#     Explanation: There are 4 differents ways to stay at index 0 after 3 steps.
#     Right, Left, Stay
#     Stay, Right, Left
#     Right, Stay, Left
#     Stay, Stay, Stay
#
#
# **Example 2:**
#
#
#
#     Input: steps = 2, arrLen = 4
#     Output: 2
#     Explanation: There are 2 differents ways to stay at index 0 after 2 steps
#     Right, Left
#     Stay, Stay
#
#
# **Example 3:**
#
#
#
#     Input: steps = 4, arrLen = 2
#     Output: 8
#
#
#
#
# **Constraints:**
#
#   * `1 <= steps <= 500`
#   * `1 <= arrLen <= 10^6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export rne_park_forward, rne_park_backward, rne_park

function Adj(G, g)
    R = G[1:3, 1:3]
    p = G[1:3, 4]
    [        R  zeros(eltype(G),3,3)
     skew(p)*R                     R] * g
end

function Adjdual(G, g)
    R = G[1:3, 1:3]
    p = G[1:3, 4]
    [        R zeros(eltype(G),3,3)
     skew(p)*R                   R]' * g
end

function adj(g, h)
    wg = g[1:3, 1]
    vg = g[4:6, 1]
    [skew(wg) zeros(eltype(g),3,3)
     skew(vg)             skew(wg)] * h
end

function adjdual(g, h)
    wg = g[1:3, 1]
    vg = g[4:6, 1]
    [skew(wg) zeros(eltype(g),3,3)
     skew(vg)             skew(wg)]' * h
end

function rne_park_forward(q, dq, ddq, Ti_inv, S, g, ifunc=identity)

    typ = eltype(q)
    dof = length(q)

    V = [Array{typ}(undef, 6) for _ in 1:dof+1]
    dV = [Array{typ}(undef, 6) for _ in 1:dof+1]

    V[1] = zeros(typ, 6)
    dV[1] = [zeros(typ, 3); -g]

    # Forward
    for i in 1:dof
        V[i+1] = ifunc(Adj(Ti_inv[i], V[i])) + ifunc(S[i]*dq[i])
        V[i+1] = ifunc(V[i+1])

        dV[i+1] = ifunc(S[i]*ddq[i]) + ifunc(Adj(Ti_inv[i], dV[i])) +
                ifunc(adj(ifunc(Adj(Ti_inv[i], V[i])), ifunc(S[i]*dq[i])))
        dV[i+1] = ifunc(dV[i+1])
    end

    V, dV
end

function rne_park_backward(Tdh_inv, S, L, l, m, V, dV, ifunc=identity)

    typ = eltype(m)
    dof = length(m)

    # extend Tdh_inv so that Tdh_inv[dof+1] return identity
    Tdh_inv = [Tdh_inv..., Matrix{Int}(I, 4, 4)]

    F = [Array{typ}(undef, 6) for _ in 1:dof+1]

    F[dof+1] = zeros(typ, 6)

    tau = zeros(typ, dof)

    # Backward
    for i in reverse(1:dof)
        Llm = [       L[i]        skew(l[i])
               -skew(l[i]) m[i].*Matrix{Int}(I, 3, 3)]

        F[i] = Adjdual(Tdh_inv[i+1], F[i+1]) + Llm*dV[i+1] - adjdual(V[i+1],  Llm*V[i+1])
        F[i] = ifunc(F[i])

        tau[i] = ifunc( dot(S[i], F[i]) )
    end

    tau
end

function rne_park(q, dq, ddq, Ti_inv, S, g, L, l, m, ifunc=identity)
    V, dV = rne_park_forward(q, dq, ddq, Ti_inv, S, g, ifunc)
    tau = rne_park_backward(Ti_inv, S, L, l, m, V, dV, ifunc)
    tau
end


text:
using JSON

export split_in_clades, dates_to_tree!, closest_tree_leaf, cluster_lbi_maxima!

"""
Split strains in clades using a timetree. Two leaves `n1` and `n2` belong to the same clade if and only if `lca(n1,n2)` lives after `datethr`.
"""
function split_in_clades(t::Tree{LBIData}, strainnames::Array{String,1}, datethr)
	clades = Dict()
	for s in strainnames
		found = false
		for id in keys(clades)
			if lca(t.lleaves[s], t.lleaves[id]).data.date >= datethr
				push!(clades[id], s)
				found = true
				break
			end
		end
		if !found
			clades[s] = [s]
		end
	end
	return clades
end

function split_in_clades(t::Tree{LBIData}, strainnames)
	println("Splitting in clades using Kmeans")
end

"""
Stores the dates contained in parsed json file `json` into `t`. The json file should be the output of `augur refine --timetree`
"""
function dates_to_tree!(t::Tree{LBIData}, json::Dict)
	for (n,dat) in json["nodes"]
		if haskey(t.lnodes, n)
			t.lnodes[n].data.date = Date(dat["date"])
		end
	end
end
dates_to_tree!(t::Tree{LBIData}, jsonfile::String) = dates_to_tree!(t, JSON.Parser.parsefile(jsonfile))

"""
Honestly, I see nothing super easy but going through all the leaves.

If I wanted to actually make this efficient, it would need to be some variation of Dijkstra's algorithm: go in a circle around `r` until you encounter a leaf. Actually it's quite easy on a tree : since you can't loop, you don't have to keep a registry of already visited nodes.
"""
function closest_tree_leaf(r::TreeNode, t::Tree)
	m = Inf
	out = ""
	for n in values(t.leaves)
		tmp = TreeTools.node_divtime(r, n)
		if tmp < m
			m = tmp
			out = n.label
		end
	end
	return (out, m)
end
closest_tree_leaf(r::String, t::Tree) = closest_tree_leaf(t.lnodes[r], t)

function closest_tree_leaf(r::TreeNode, t::Tree, candidates::Array{<:AbstractString,1})
	m = Inf
	out = ""
	for s in candidates
		n = t.lnodes[s]
		tmp = TreeTools.node_divtime(r, n)
		if tmp < m
			m = tmp
			out = n.label
		end
	end
	return (out, m)
end
closest_tree_leaf(r::String, t::Tree, candidates::Array{<:AbstractString,1}) = closest_tree_leaf(t.lnodes[r], t, candidates)

"""
	cluster_lbi_maxima!(t::Tree, fp)

Cluster sequences of each datebin in `fp.datebin` using local maximas of the LBI for the tree of that datebin.
"""
function cluster_lbi_maxima!(t::Tree, fp)
	LLM = local_lbi_maximas!(t, fp);
	clusters = Dict()
	lbis = Dict()
	for (db, llm) in LLM
		centers = [x.label for x in llm]
		clusters[db] = Dict(k=>[] for k in centers)
		lbis[db] = Dict(k=>v for (k,v) in llm)
		pop = fp.datebin[db]
		# For each strain in the current pop, find which LBI maximum it is closest to
		for s in pop
			closest_center = closest_tree_leaf(t.lleaves[s[:strain]], t, centers)[1]
			push!(clusters[db][closest_center], s[:strain])
		end
	end
	for (db,cl) in clusters
		for (k,c) in cl
			if isempty(c)
				delete!(clusters[db], k)
				delete!(lbis[db], k)
			end
		end
	end
	return clusters, lbis
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
include("bitindex.jl")
include("bitpar-compiler.jl")

@inline function reversebits(x::T, ::BitsPerSymbol{2}) where T <: Base.BitUnsigned
     mask = 0x33333333333333333333333333333333 % T
     x = ((x >> 2) & mask) | ((x & mask) << 2)
     return reversebits(x, BitsPerSymbol{4}())
end

@inline function reversebits(x::T, ::BitsPerSymbol{4}) where T <: Base.BitUnsigned
     mask = 0x0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F % T
     x = ((x >> 4) & mask) | ((x & mask) << 4)
     return bswap(x)
end

reversebits(x::T, ::BitsPerSymbol{8}) where T <: Base.BitUnsigned = bswap(x)

@inline function complement_bitpar(x::Unsigned, ::T) where {T<:NucleicAcidAlphabet{2}}
    return ~x
end

@inline function complement_bitpar(x::Unsigned, ::T) where {T<:NucleicAcidAlphabet{4}}
    return (
        ((x & repeatpattern(typeof(x), 0x11)) << 3) | ((x & repeatpattern(typeof(x), 0x88)) >> 3) |
        ((x & repeatpattern(typeof(x), 0x22)) << 1) | ((x & repeatpattern(typeof(x), 0x44)) >> 1)
    )
end

@inline function gc_bitcount(x::Unsigned, ::BitsPerSymbol{2})
    msk = repeatpattern(typeof(x), 0x55)
    c = x & msk
    g = (x >> 1) & msk
    return count_ones(c ⊻ g)
end

@inline function gc_bitcount(x::Unsigned, ::BitsPerSymbol{4})
    a =  x & repeatpattern(typeof(x), 0x11)
    c = (x & repeatpattern(typeof(x), 0x22)) >> 1
    g = (x & repeatpattern(typeof(x), 0x44)) >> 2
    t = (x & repeatpattern(typeof(x), 0x88)) >> 3
    return count_ones((c | g) & ~(a | t))
end

@inline function mismatch_bitcount(a::UInt64, b::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    return count_nonzero_nibbles(a ⊻ b)
end

@inline function mismatch_bitcount(a::UInt64, b::UInt64, ::T) where {T<:NucleicAcidAlphabet{2}}
    return count_nonzero_bitpairs(a ⊻ b)
end

@inline function match_bitcount(a::UInt64, b::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    return count_0000_nibbles(a ⊻ b)
end

@inline function match_bitcount(a::UInt64, b::UInt64, ::T) where {T<:NucleicAcidAlphabet{2}}
    return count_00_bitpairs(a ⊻ b)
end

@inline function ambiguous_bitcount(x::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    return count_nonzero_nibbles(enumerate_nibbles(x) & 0xEEEEEEEEEEEEEEEE)
end

@inline function ambiguous_bitcount(a::UInt64, b::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    return count_nonzero_nibbles((enumerate_nibbles(a) | enumerate_nibbles(b)) & 0xEEEEEEEEEEEEEEEE)
end

@inline function gap_bitcount(x::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    return count_zero_nibbles(x)
end

@inline function gap_bitcount(a::UInt64, b::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    # Count the gaps in a, count the gaps in b, subtract the number of shared gaps.
    return count_0000_nibbles(a) + count_0000_nibbles(b) - count_0000_nibbles(a | b)
end

@inline function certain_bitcount(x::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    x = enumerate_nibbles(x)
    x = x ⊻ 0x1111111111111111
    return count_0000_nibbles(x)
end

@inline function certain_bitcount(a::UInt64, b::UInt64, ::T) where {T<:NucleicAcidAlphabet{4}}
    x = enumerate_nibbles(a) ⊻ 0x1111111111111111
    y = enumerate_nibbles(b) ⊻ 0x1111111111111111
    return count_0000_nibbles(x | y)
end

@inline count_a(x::Unsigned) = count_00_bitpairs(x)
@inline count_c(x::Unsigned) = count_01_bitpairs(x)
@inline count_g(x::Unsigned) = count_10_bitpairs(x)
@inline count_t(x::Unsigned) = count_11_bitpairs(x)


text:
using KernelAbstractions: @kernel, @index
using KernelAbstractions.Extras.LoopInfo: @unroll

using Oceananigans.Architectures: arch_array, architecture
using Oceananigans.Operators: Δzᶜᶜᶜ

const SingleColumnGrid = AbstractGrid{<:AbstractFloat, <:Flat, <:Flat, <:Bounded}

"""
    regrid!(a, b)

Regrid field `b` onto the grid of field `a`.

!!! warning "Functionality limitation"
    Currently `regrid!` only regrids in the vertical ``z`` direction and works only on
    fields that have data only in ``z`` direction.

Example
=======

Generate a tracer field on a vertically stretched grid and regrid it on a regular grid.

```jldoctest
using Oceananigans

Nz, Lz = 2, 1.0
topology = (Flat, Flat, Bounded)

input_grid = RectilinearGrid(size=Nz, z = [0, Lz/3, Lz], topology=topology, halo=1)
input_field = CenterField(input_grid)
input_field[1, 1, 1:Nz] = [2, 3]

output_grid = RectilinearGrid(size=Nz, z=(0, Lz), topology=topology, halo=1)
output_field = CenterField(output_grid)

regrid!(output_field, input_field)

output_field[1, 1, :]

# output
4-element OffsetArray(::Vector{Float64}, 0:3) with eltype Float64 with indices 0:3:
 0.0
 2.333333333333334
 3.0
 0.0
```
"""
regrid!(a, b) = regrid!(a, a.grid, b.grid, b)

function we_can_regrid(a, target_grid, source_grid, b)
    # Only 1D regridding in the vertical is supported, so check that
    #   1. source and target grid are in the same "class" and
    #   2. source and target Field have same horizontal size
    typeof(source_grid).name.wrapper === typeof(target_grid).name.wrapper &&
        size(a)[1:2] === size(b)[1:2] && return true

    return false
end

function regrid!(a, target_grid, source_grid, b)
    if we_can_regrid(a, target_grid, source_grid, b)
        arch = architecture(a)
        source_z_faces = znodes(Face, source_grid)

        event = launch!(arch, target_grid, :xy, _regrid!, a, b, target_grid, source_grid, source_z_faces)
        wait(device(arch), event)

        return a
    else
        msg = """Regridding
                 $(summary(b)) on $(summary(source_grid))
                 to $(summary(a)) on $(summary(target_grid))
                 is not supported."""

        return throw(ArgumentError(msg))
    end
end

#####
##### Regridding for single column grids
#####

@kernel function _regrid!(target_field, source_field, target_grid, source_grid, source_z_faces)
    i, j = @index(Global, NTuple)

    Nx_target, Ny_target, Nz_target = size(target_grid)
    Nx_source, Ny_source, Nz_source = size(source_grid)
    i_src = ifelse(Nx_target == Nx_source, i, 1)
    j_src = ifelse(Ny_target == Ny_source, j, 1)

    @unroll for k = 1:target_grid.Nz
        @inbounds target_field[i, j, k] = 0

        z₋ = znode(Center(), Center(), Face(), i, j, k,   target_grid)
        z₊ = znode(Center(), Center(), Face(), i, j, k+1, target_grid)

        # Integrate source field from z₋ to z₊
        k₋_src = searchsortedfirst(source_z_faces, z₋)
        k₊_src = searchsortedfirst(source_z_faces, z₊) - 1

        # Add contribution from all full cells in the integration range
        @unroll for k_src = k₋_src:k₊_src
            @inbounds target_field[i, j, k] += source_field[i_src, j_src, k_src] * Δzᶜᶜᶜ(i_src, j_src, k_src, source_grid)
        end

        zk₋_src = znode(Center(), Center(), Face(), i_src, j_src, k₋_src, source_grid)
        zk₊_src = znode(Center(), Center(), Face(), i_src, j_src, k₊_src+1, source_grid)

        # Add contribution to integral from fractional bottom part,
        # if that region is a part of the grid.
        @inbounds target_field[i, j, k] += source_field[i_src, j_src, k₋_src - 1] * (zk₋_src - z₋)

        # Add contribution to integral from fractional top part
        @inbounds target_field[i, j, k] += source_field[i_src, j_src, k₊_src] * (z₊ - zk₊_src)

        @inbounds target_field[i, j, k] /= Δzᶜᶜᶜ(i, j, k, target_grid)
    end
end


text:
###################################################
## output.jl
##      Output data and solutions to ASCII readable format
## Authors: Arthur Delarue, Sébastien Martin, 2018
###################################################

"""
	Convert number of seconds since midnight to DateTime format
"""
using Dates

function totime(x)
	# newx = round(x, 3)  #shuai
	newx  = round(x, digits=3)
	hours = div(newx, 3600)
	minutes = div(newx % 3600, 60)
	seconds = floor(newx % 60)
	thousands = Int(round(1000 * (newx - floor(newx))))
	return Dates.Time(hours, minutes, seconds, thousands)
end

"""
    Output full routing solution to dataframe
"""
function output(data::SchoolBusData)
    !data.withFinalBuses && error("Final routes not computed")
    bus_id = Int[]
    start_loc = Int[]
	end_loc = Int[]
    start_time = Dates.Time[]
	end_time = Dates.Time[]
	action = String[]
	for (j, bus) in enumerate(data.buses)
		for (i, school) in enumerate(bus.schools)
			route = bus.routes[i]
			src = [data.stops[school][data.routes[school][route].stops[end]].uniqueId + 100000]
			dst = [data.schools[school].id + 200000]
			currenttime = data.schools[school].starttime - data.schools[school].dwelltime
			dst_t = [totime(currenttime)]
			currenttime -= traveltime(data,
			                          data.stops[school][data.routes[school][route].stops[end]],
			                          data.schools[school])
			src_t = [totime(currenttime)]
			act = ["travel"]
			for k in length(data.routes[school][route].stops):-1:2
				stop = data.routes[school][route].stops[k]
				# add waiting time at stop
				pushfirst!(act, "pickup")
				pushfirst!(dst_t, src_t[1])
				pushfirst!(dst, src[1])
				pushfirst!(src, src[1])
				currenttime -= stopTime(data, data.stops[school][stop])
				pushfirst!(src_t, totime(currenttime))
				# add travel time from previous stop
				pushfirst!(act, "travel")
				pushfirst!(dst_t, src_t[1])
				pushfirst!(dst, src[1])
				prev = data.routes[school][route].stops[k-1]
				pushfirst!(src, data.stops[school][prev].uniqueId + 100000)
				currenttime -= traveltime(data, data.stops[school][prev], data.stops[school][stop])
				pushfirst!(src_t, totime(currenttime))
			end
			stop = data.routes[school][route].stops[1]
			pushfirst!(act, "pickup")
			pushfirst!(dst_t, src_t[1])
			pushfirst!(dst, src[1])
			pushfirst!(src, src[1])
			currenttime -= stopTime(data, data.stops[school][stop])
			pushfirst!(src_t, totime(currenttime))
			if i == 1
				pushfirst!(act, "deadhead")
				pushfirst!(dst_t, src_t[1])
				pushfirst!(dst, src[1])
				pushfirst!(src, 900000 + bus.yard)
				currenttime -= traveltime(data, data.yards[bus.yard], data.stops[school][stop])
				pushfirst!(src_t, totime(currenttime))
			else
				pushfirst!(act, "deadhead")
				pushfirst!(dst_t, src_t[1])
				pushfirst!(dst, src[1])
				pushfirst!(src, 200000 + data.schools[bus.schools[i-1]].id)
				currenttime -= traveltime(data, data.schools[bus.schools[i-1]],
				                          data.stops[school][stop])
				pushfirst!(src_t, totime(currenttime))
				pushfirst!(act, "wait")
				pushfirst!(dst_t, src_t[1])
				pushfirst!(dst, src[1])
				pushfirst!(src, src[1])
				currenttime = data.schools[bus.schools[i-1]].starttime
				pushfirst!(src_t, totime(currenttime))
				pushfirst!(act, "dropoff")
				pushfirst!(dst_t, src_t[1])
				pushfirst!(dst, src[1])
				pushfirst!(src, src[1])
				currenttime -= data.schools[bus.schools[i-1]].dwelltime
				pushfirst!(src_t, totime(currenttime))
			end
			# add this school to the itinerary
			append!(start_loc, src)
			append!(end_loc, dst)
			append!(start_time, src_t)
			append!(end_time, dst_t)
			append!(bus_id, j * ones(Int, length(src)))
			append!(action, act)
		end
		# add final school and return to the yard
		append!(action, ["dropoff", "deadhead"])
		append!(bus_id, [j, j])
		push!(start_loc, end_loc[end])
		push!(start_time, end_time[end])
		push!(end_loc, end_loc[end])
		currenttime = data.schools[bus.schools[end]].starttime
		push!(end_time, totime(currenttime))
		push!(start_loc, end_loc[end])
		push!(start_time, end_time[end])
		push!(end_loc, bus.yard + 900000)
		currenttime += traveltime(data, data.schools[bus.schools[end]], data.yards[bus.yard])
		push!(end_time, totime(currenttime))
	end
	df = DataFrame()
	df[:bus_id] = bus_id
	df[:start_loc] = start_loc
	df[:end_loc] = end_loc
	df[:action] = action
	df[:start_time] = start_time
	df[:end_time] = end_time
	return df
end

data1 = output(data1)
g


text:
ENV["OMP_NUM_THREADS"] = 4
using Pardiso
using Elsa, BenchmarkTools, Arpack, LinearAlgebra, LinearMaps, KrylovKit, ArnoldiMethod,
      IterativeSolvers, SuiteSparse

function solver()
    ps = MKLPardisoSolver()
    set_matrixtype!(ps, Pardiso.COMPLEX_HERM_INDEF)
    # set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    pardisoinit(ps)
    return ps
end
function lmap(h::AbstractMatrix{Tv}) where {Tv}
    ps = solver()
    A_pardiso = get_matrix(ps, h, :C);
    cols, rows = size(h)
    # b = zeros(eltype(h), cols)
    b = Vector(h[:,3])
    set_phase!(ps, Pardiso.ANALYSIS)
    # set_perm!(ps, randperm(n))
    pardiso(ps, A_pardiso, b)
    set_phase!(ps, Pardiso.NUM_FACT)
    pardiso(ps, A_pardiso, b)
    set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    LinearMap{Tv}((x, y) -> pardiso(ps, x, A_pardiso, y), cols, rows,
        ismutating = true, ishermitian = true)
end

function lmap2(h::AbstractMatrix{Tv}) where {Tv}
    cols, rows = size(h)
    # F = ldlt(h)
    F = lu(h)
    LinearMap{Tv}((x, y) -> ldiv!(x, F, y), cols, rows,
        ismutating = true, ishermitian = true)
end

sortreal(list, n = min(length(list), 10); rev = false) = sort(sort(list, by = abs, rev = rev)[1:n], by = real)

# sys = System(LatticePresets.honeycomb(), Model(Hopping(1))) |> grow(region = Region(:circle, 300))
sys = System(LatticePresets.honeycomb(), Model(Onsite(.2), Hopping(1, range = 1/√3))) |> grow(supercell = 30)
h = hamiltonian(sys, k = (0.21321,0.234213));
@time l = lmap(h);
@time l2 = lmap2(h);
lit = LOBPCGIterator(l, I, true, rand(eltype(h), size(h,1), 10));
lit2 = LOBPCGIterator(l2, I, true, rand(eltype(h), size(h,1), 10));

lob(lit) = sortreal(1 ./ lobpcg!(lit).λ);
arp(h) = sortreal(eigs(h, nev = 10, sigma = .01im)[1]);
ame(l) = sortreal(1 ./ partialschur(l, nev = 10)[1].eigenvalues);
kry(l) = sortreal(1 ./ eigsolve(x -> l * x, complex.(rand(size(l,1))), 10)[1]);

@btime lob($lit)
@btime arp($h)
@btime ame($l)
@btime kry($l)

@btime lob($lit2)
@btime arp($h)
@btime ame($l2)
@btime kry($l2)


text:
export searchMatrix

"""
[240] Medium [Search a 2D Matrix II](https://leetcode.com/problems/search-a-2d-matrix-ii/description/)

Write an efficient algorithm that searches for a value in an *m* x *n* matrix. This matrix has the following properties:

+ Integers in each row are sorted in ascending from left to right.
+ Integers in each column are sorted in ascending from top to bottom.

## Example:

Consider the following matrix:

```
[
    [1,   4,  7, 11, 15],
    [2,   5,  8, 12, 19],
    [3,   6,  9, 16, 22],
    [10, 13, 14, 17, 24],
    [18, 21, 23, 26, 30]
]
```

Given target = `5`, return `true`.

Given target = `20`, return `false`.
"""
function searchMatrix(matrix::Matrix{Int}, target::Int)::Bool
    row, col = size(matrix)
    i, j = row, 1

    while (i > 0 && j <= col)
        last = matrix[i, j]

        if (last == target)
            return true
        elseif last > target
            i -= 1
        else
            j += 1
        end
    end

    false
end


text:
# clustering julia

# load libs
using VegaLite, VegaDatasets, DataFrames, Statistics, JSON, CSV, Distances, Clustering

# loading the data
download("https://raw.githubusercontent.com/ageron/handson-ml/master/datasets/housing/housing.csv","newhouses.csv");
@time houses = CSV.read("newhouses.csv", DataFrame);

#
cd("G:\\My Drive\\18. Github Repo\\DataScience")
describe(houses);
names(houses)

# vegalite
cali_shape = JSON.parsefile("data/california-counties.json")
VV = VegaDatasets.VegaJSONDataset(cali_shape,"data/california-counties.json")

# plotting using the median house price
@vlplot(width=500, height=300) +
@vlplot(
    mark={:geoshape, fill=:black, stroke=:white},
    data={values=VV, format={type=:topojson, feature=:cb_2015_california_county_20m}},
    projection={type=:albersUsa}) +
@vlplot(:circle, data=houses, projection={type=:albersUsa}, longitude="longitude:q", latitude="latitude:q", size={value=12}, color="median_house_value:q")

#
bucketprice = Int.(div.(houses[!,:median_house_value],50000));
insertcols!(houses,3, :cprice => bucketprice);

@vlplot(width=500, height=300) +
@vlplot(
    mark={:geoshape, fill=:black, stroke=:white},
    data={values=VV, format={type=:topojson, feature=:cb_2015_california_county_20m}},
    projection={type=:albersUsa}) +
@vlplot(:circle, data=houses, projection={type=:albersUsa}, longitude="longitude:q", latitude="latitude:q", size={value=12}, color="cprice:n")

# clustering
X = Matrix{Float64}(houses[!, [:latitude, :longitude]]);
@time C = kmeans(X', 10);
insertcols!(houses, 3, :cluster10 => C.assignments);

@vlplot(width=500, height=300) +
@vlplot(
    mark={:geoshape, fill=:black, stroke=:white},
    data={values=VV, format={type=:topojson, feature=:cb_2015_california_county_20m}},
    projection={type=:albersUsa}) +
@vlplot(:circle, data=houses, projection={type=:albersUsa}, longitude="longitude:q", latitude="latitude:q", size={value=12}, color="cluster10:n")

# using the distances
D = pairwise(Euclidean(), X', X', dims=2);
@time K = kmedoids(D, 10);
insertcols!(houses, 3, :medoid10 => K.assignments);

@vlplot(width=500, height=300) +
@vlplot(
    mark={:geoshape, fill=:black, stroke=:white},
    data={values=VV, format={type=:topojson, feature=:cb_2015_california_county_20m}},
    projection={type=:albersUsa}) +
@vlplot(:circle, data=houses, projection={type=:albersUsa}, longitude="longitude:q", latitude="latitude:q", size={value=12}, color="medoid10:n")

# hierachiecal clustering
@time K = hclust(D);
@time L = cutree(K, k=10);
insertcols!(houses, 3, :hclust => L);

@vlplot(width=500, height=300) +
@vlplot(
    mark={:geoshape, fill=:black, stroke=:white},
    data={values=VV, format={type=:topojson, feature=:cb_2015_california_county_20m}},
    projection={type=:albersUsa}) +
@vlplot(:circle, data=houses, projection={type=:albersUsa}, longitude="longitude:q", latitude="latitude:q", size={value=12}, color="hclust:n")

max_stars_repo_path: 05. Clustering.jl

max_issues_repo_path: 05. Clustering.jl

max_forks_repo_path: 05. Clustering.jl


text:
## info on the current run for the extraction
## analysis of the .csv produced

using DataFrames, Query
using CSV, Glob, Dates
using Printf

function perf_status(file)
    if isfile(file)
        perf= CSV.File(file, delim=";")
        nstarcluster= sum(perf.nmax)

        x = @from i in perf begin
            @where i.cycle == 1
            @select i
            @collect DataFrame
        end
        nstartotal= sum(x.nstar)
    else
        nstarcluster= 0
        nstartotal= 0
    end
    return(nstarcluster, nstartotal)
end

function specialstr(str, CODE)
    d= Dict( "PURPLE" => "\033[95m" ,
             "CYAN" => "\033[96m" ,
             "DARKCYAN" => "\033[36m" ,
             "BLUE" => "\033[94m" ,
             "GREEN" => "\033[92m" ,
             "YELLOW" => "\033[93m" ,
             "RED" => "\033[91m" ,
             "BOLD" => "\033[1m" ,
             "UNDERLINE" => "\033[4m" ,
             "END" => "\033[0m")

    if haskey(d, CODE)
        res= @sprintf("%s%s%s", d[CODE], str, d["END"])
        return(res)
    else
        return(str)
    end
end

function bold(s) specialstr(s,"BOLD") end
function yellow(s) specialstr(s,"YELLOW") end

function main(csvlist)
    println("## Extraction status")
    println("## Based on files *.done.csv and *.perf.csv")
    println("## $(now())")
    println("##")

    nvotTotal= 0 ; ncycleTotal= 0 ; nstarclusterTotal= 0 ; nstarTotal= 0
    for file in csvlist
        done= CSV.File(file, delim=";")

        nvot= length(done.votname)
        ncycle= sum(done.cycle)
        fileperf= split(file,".")[1]*".perf.csv"
        nstarcluster, nstartotal= perf_status(fileperf)

        println("### $file")
        println("### votables: $nvot")
        println("### cycles: $ncycle")
        println("### stars extracted: $nstarcluster")
        println("### total star analyzed: $nstartotal")
        println("###")

        nvotTotal += nvot
        ncycleTotal += ncycle
        nstarclusterTotal += nstarcluster
        nstarTotal += nstartotal
    end
    println("##")
    println(yellow("## SUMMARY:"))
    println(yellow("## votables: $nvotTotal"))
    println(yellow("## cycles: $ncycleTotal"))
    println(yellow("## stars extracted: $nstarclusterTotal"))
    println(yellow("## stars analyzed: $nstarTotal"))

end

list= glob("*.done.csv")
main(list)

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
export
    AbstractState,
    StateVector,
    StateMatrix,

    𝛒,
    purity

abstract type AbstractState end

################
# state vector #
################

"""
    StateVector{T <: Number} <: AbstractState

Vector representation for pure quantum state.
There are various constructures to construct different pure quantum states.
"""
mutable struct StateVector{T<:Number}<:AbstractState
    v::Vector{T}
    dim::Int64
end

function Base.show(io::IO, state::StateVector{T}) where {T}
    println(io, "StateVector{$T}(dim=$(state.dim), vec=[")
    Base.print_matrix(IOContext(io, :limit=>true), state.v)
    print(io, "\n])")
end

function Base.:+(s1::StateVector, s2::StateVector)
    s1, s2 = (s1.dim ≥ s2.dim) ? (s1, s2) : (s2, s1)
    new_state = copy(s1)
    new_state.v[1:length(s2.v)] .+= s2.v

    return new_state
end

function Base.:*(c::Number, s::StateVector)
    new_state = copy(s)
    new_state.v .*= c

    return new_state
end

Base.:*(s::StateVector, c::Number) = c * s

"""
    Base.vec(state::StateVector{<:Number})

To get the vector of a pure quantum state.

# Examples
```julia-repl
julia> state = FockState(1);

julia> vec(state)
$DIM-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im
     ⋮
 0.0 + 0.0im
```
"""
Base.vec(state::StateVector{<:Number}) = state.v

"""
    𝛒(state::StateVector{<:Number})

To get the density matrix of a pure quantum state.

# Examples
```julia-repl
julia> state = FockState(1);

julia> 𝛒(state)
$DIM×$DIM Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
    ⋮                             ⋱
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
```
"""
𝛒(state::StateVector{<:Number}) = state.v * state.v'

"""
    purity(state::StateVector{<:Number})

Calculate purity for quantum state in vector representation.

# Examples
```jldoctest
julia> state = VacuumState();

julia> purity(state)
1.0
```
"""
function purity(state::StateVector{<:Number})
    𝛒 = state.v * state.v'
    𝛒 /= tr(𝛒)

    return real(tr(𝛒^2))
end

"""
    Base.copy(state::StateVector{<:Number})

Return a new instance of a `StateVector`

# Examples
```jldoctest
julia> state = VacuumState();

julia> new_state = copy(state);
```
"""
function Base.copy(state::StateVector{T}) where {T<:Number}
    return StateVector{T}(copy(state.v), state.dim)
end

################
# state matrix #
################

"""
    StateMatrix{T <: Number} <: AbstractState

Density Matrix representation for pure and mixed quantum state.
There are various constructures to construct different pure and mixed quantum states.
"""
mutable struct StateMatrix{T<:Number}<:AbstractState
    𝛒::Matrix{T}
    dim::Int64
end

function Base.show(io::IO, state::StateMatrix{T}) where {T}
    println(io, "StateMatrix{$T}(dim=$(state.dim), 𝛒=[")
    Base.print_matrix(IOContext(io, :limit=>true), state.𝛒)
    print(io, "\n])")
end

function Base.:+(s1::StateMatrix, s2::StateMatrix)
    s1, s2 = (s1.dim ≥ s2.dim) ? (s1, s2) : (s2, s1)
    new_state = copy(s1)
    new_state.𝛒[1:size(s2.𝛒, 1), 1:size(s2.𝛒, 2)] .+= s2.𝛒

    return new_state
end

function Base.:*(c::Number, s::StateMatrix)
    new_state = copy(s)
    new_state.𝛒 .*= c

    return new_state
end

Base.:*(s::StateMatrix, c::Number) = c * s

"""
    StateMatrix(state::StateVector{<:Number})

Convert a `StateVector` to a `StateMatrix`.

# Examples
```jldoctest
julia> state = VacuumState();

julia> StateMatrix(state);
```
"""
function StateMatrix(state::StateVector{T}) where {T<:Number}
    𝛒 = state.v * state.v'

    return StateMatrix{T}(𝛒, state.dim)
end

"""
    𝛒(state::StateMatrix{<:Number})

To get the density matrix of a pure quantum state.

# Examples
```julia-repl
julia> state = FockState(1);

julia> 𝛒(state)
$DIM×$DIM Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
    ⋮                             ⋱
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
```
"""
𝛒(state::StateMatrix{<:Number}) = state.𝛒

"""
    purity(state::StateMatrix{<:Number})

Calculate purity for quantum state in density matrix representation.

# Examples
```jldoctest
julia> state = VacuumState(rep=StateMatrix);

julia> purity(state)
1.0
```
"""
function purity(state::StateMatrix{<:Number})
    𝛒 = state.𝛒
    𝛒 /= tr(𝛒)

    return real(tr(𝛒^2))
end

"""
    Base.copy(state::StateMatrix{<:Number})

Return a new instance of a `StateMatrix`

# Examples
```jldoctest
julia> state = VacuumState(rep=StateMatrix);

julia> new_state = copy(state);
```
"""
function Base.copy(state::StateMatrix{T}) where {T<:Number}
    return StateMatrix{T}(copy(state.𝛒), state.dim)
end


text:
"""
    CA_matches(sequence::String, CAs::Vector{Rule}, CA_weights::Dict{String,Int64}, occurrence_weighting::Bool)

Counts the number of CA's matched by a sequence (only support for simple rules).

# Arguments
- `sequence::String`: sequence to count matches.
- `CAs::Vector{Rule}`: list of all CA's.
- `CA_weights::Dict{String,Int64}`: weights to use for CA counts.
- `occurrence_weighting::Bool`: whether to use occurrence weighting during counting.
"""
function CA_matches(sequence::String, CAs::Vector{Rule}, CA_weights::Dict{String,Int64}, occurrence_weighting::Bool ; protein::Bool=false)

    letter_transformations = Dict{Char,Vector{Char}}('A'=>['A'],
                                                     'T'=>['T'],
                                                     'C'=>['C'],
                                                     'G'=>['G'],
                                                     'U'=>['U'],
                                                     'R'=>['A','G'],
                                                     'Y'=>['C','T'],
                                                     'S'=>['G','C'],
                                                     'W'=>['A','T'],
                                                     'K'=>['G','T'],
                                                     'M'=>['A','C'],
                                                     'B'=>['C','G','T'],
                                                     'D'=>['A','G','T'],
                                                     'H'=>['A','C','T'],
                                                     'V'=>['A','C','G'],
                                                     'N'=>['A','T','C','G'],
                                                     '-'=>['-'])

    used_idxs = Array{Int,1}()

    # Initialize CAOS score
    total_score = 0

    # Iterate through each rule
    for rule in CAs

        # Initialize variables
        is_match = true
        score = 1

        # Check if the sequence follows the rule
        for (iter,idx) in enumerate(rule.idxs)

            if idx in used_idxs
                is_match = false
                break
            end

            if haskey(letter_transformations, sequence[idx]) && !protein
                seq_letters = letter_transformations[sequence[idx]]
            else
                seq_letters = sequence[idx]
            end

            if !(rule.char_attr[iter] in seq_letters)
                is_match = false
                break
            end
        end

        # Keep track of score
        if is_match

            push!(used_idxs, (rule.idxs)[1])

            # Add value of pure/private rule
            if rule.is_pure
                if length(rule.char_attr) == 1
                    score *= CA_weights["sPu"]
                else
                    score *= CA_weights["cPu"]
                end
            else
                if length(rule.char_attr) == 1
                    score *= CA_weights["sPr"]
                    if occurrence_weighting
                        score *= rule.occurances/rule.num_group
                    end
                else
                    score *= CA_weights["cPr"]
                    if occurrence_weighting
                        score *= rule.occurances/rule.num_group
                    end
                end
            end

            # Add to total score
            total_score += score
        end
    end

    # Return CAOS score
    return Int(round(total_score))
end

"""
    classify_sequence(sequence::String, tree::Node, CA_weights::Dict{String,Int64},
    all_CA_weights::Dict{Int64,Dict{String,Int64}}, occurrence_weighting::Bool,
    depth::Int64, tiebreaker::Vector{Dict{String,Int64}} ; blast_results=["Fake Label"], combo_classification::Bool=false, protein::Bool=false)

Classifies an input sequence given a phylogentic tree.

# Arguments
- `sequence::String`: sequence to count matches.
- `tree::Node`: the tree represented as a Node.
- `CA_weights::Dict{String,Int64}`: weights to use for CA counts.
- `all_CA_weights::Dict{Int64,Dict{String,Int64}}`: all sets of weights to use for CA counts.
- `occurrence_weighting::Bool`: whether to use occurrence weighting during counting.
- `depth::Int64`: current depth of the tree.
- `tiebreaker::Vector{Dict{String,Int64}}`: tiebreaking procedures to use.
- `blast_results=["Fake Label"]`: list of blast results.
- `combo_classification::Bool=false`: whether to use both Blast and CAOS for classification.
"""
function classify_sequence(sequence::String, tree::Node, CA_weights::Dict{String,Int64}, all_CA_weights::Dict{Int64,Dict{String,Int64}}, occurrence_weighting::Bool, depth::Int64, tiebreaker::Vector{Dict{String,Int64}} ; blast_results=["Fake Label"], combo_classification::Bool=false, protein::Bool=false)

    if haskey(all_CA_weights, depth)
        CA_weights = all_CA_weights[depth]
    end

    # If at a terminal node, return the taxa label
    if length(tree.children) == 0
        return tree.taxa_label

    # Get the CA score for each child
    else

        if false && combo_classification && (length(blast_results) >= 10)
            for child in tree.children
                descendents = get_descendents(child)
                all_blast = true
                for blast_result in blast_results[1:10]
                    if !(blast_result in descendents)
                        all_blast = false
                        break
                    end
                end
                if all_bast
                    return classify_sequence(sequence,child,CA_weights,all_CA_weights,occurrence_weighting,depth+1,tiebreaker, blast_results=blast_results, combo_classification=combo_classification, protein=protein)
                end
            end
        end

        child_CA_score = Array{Int,1}()
        for child in tree.children
            push!(child_CA_score, CA_matches(sequence, child.CAs, CA_weights, occurrence_weighting, protein=protein))
        end

        max_child_idx = argmax(child_CA_score)

        # Select the child with highest CA score and descend in that direction
        if length(max_child_idx) == 1
            return classify_sequence(sequence,tree.children[max_child_idx[1]],CA_weights,all_CA_weights,occurrence_weighting,depth+1,tiebreaker, blast_results=blast_results, combo_classification=combo_classification, protein=protein)

        # If there is a tie, return either use the tiebreaker or return the tree
        elseif tiebreaker[1] == Dict{String,Int64}()
            if combo_classification && (length(blast_results) >= 5)
                for child in tree.children
                    descendents = get_descendents(child)
                    top_blasts = true
                    for blast_result in blast_results[1:5]
                        if !(blast_result in descendents)
                            top_blasts = false
                            break
                        end
                    end
                    if top_blasts
                        return classify_sequence(sequence,child,CA_weights,all_CA_weights,occurrence_weighting,depth+1,tiebreaker, blast_results=blast_results, combo_classification=combo_classification, protein=protein)
                    end
                end
            end
            return tree
        else
            return classify_sequence(sequence,tree,tiebreaker[1],Dict{Int,Dict{String,Int64}}(),occurrence_weighting,depth,tiebreaker[2:end], blast_results=blast_results, combo_classification=combo_classification, protein=protein)
        end
    end
end


text:
struct HitPoint
    point::AbstractVec
    normal::AbstractVec
    u::AbstractFloat
    v::AbstractFloat
    front::Bool
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
nnodes(t) = _nnodes(NodeType(t), t)
_nnodes(::LeafNode, _) = 1
_nnodes(::InnerNode, t) = 1 + mapreduce(nnodes, +, _iter(children(t)); init=0)

nleafs(t) = _nleafs(NodeType(t), t)
_nleafs(::LeafNode, _) = 1
_nleafs(::InnerNode, t) = mapreduce(nleafs, +, _iter(children(t)); init=0) + Int(nchildren(t) == 0)

treeheight(t) = _treeheight(NodeType(t), t)
_treeheight(::LeafNode, _) = 0
_treeheight(::InnerNode, t) = mapreduce(treeheight, max, _iter(children(t)); init=-1) + 1


text:
# To ensure that tests can run on travis we have to do a little
# hackadoodle here. The tests require a license file. We include
# a license file that is only valid for one day (the day when
# change is submitted).
# If there is no valid license file, we default to that file.

using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const FALLBACK_URL = "mosek://solve.mosek.com:30080"

using MosekTools
const optimizer = Mosek.Optimizer()
MOI.set(optimizer, MOI.RawOptimizerAttribute("fallback"), FALLBACK_URL)
MOI.set(optimizer, MOI.Silent(), true)

function MosekOptimizerWithFallback()
    optimizer = Mosek.Optimizer()
    MOI.set(optimizer, MOI.RawOptimizerAttribute("fallback"), FALLBACK_URL)
    MOI.set(optimizer, MOI.Silent(), true)
    return optimizer
end

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "Mosek"
end

@testset "Parameters" begin
    optimizer = MosekOptimizerWithFallback()
    MOI.set(optimizer, MOI.RawOptimizerAttribute("fallback"), FALLBACK_URL)
    @testset "Double Parameter" begin
        MOI.set(optimizer, MOI.RawOptimizerAttribute("INTPNT_CO_TOL_DFEAS"), 1e-7)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-7
        MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1e-8)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-8
        @testset "with integer value" begin
            MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1)
            @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1
        end
    end
    @testset "Integer Parameter" begin
        MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_INTPNT_MAX_ITERATIONS"), 100)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 100
        MOI.set(optimizer, MOI.RawOptimizerAttribute("INTPNT_MAX_ITERATIONS"), 200)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 200
        @testset "with enum value" begin
            MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_OPTIMIZER"), MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
            @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_OPTIMIZER")) == convert(Int32, MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
        end
    end
    @testset "String Parameter" begin
        MOI.set(optimizer, MOI.RawOptimizerAttribute("PARAM_WRITE_FILE_NAME"), "foo.txt")
        # Needs https://github.com/JuliaOpt/Mosek.jl/pull/174
        #@test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "foo.txt"
        MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_SPAR_PARAM_WRITE_FILE_NAME"), "bar.txt")
        #@test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "bar.txt"
    end
    @testset "TimeLimitSec" begin
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
        MOI.set(optimizer, MOI.TimeLimitSec(), 1.0)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_OPTIMIZER_MAX_TIME")) == 1.0
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === 1.0
        MOI.set(optimizer, MOI.TimeLimitSec(), nothing)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
    end
end

@testset "supports_incremental_interface" begin
    @test MOI.supports_incremental_interface(optimizer)
end

const config = MOIT.Config(
    Float64, atol=1e-3, rtol=1e-3,
    exclude=Any[MOI.ConstraintName, MOI.VariableBasisStatus, MOI.ConstraintBasisStatus], # result in errors for now
)

@testset "Direct optimizer tests" begin
    MOIT.runtests(optimizer, config,
        exclude=[
            # FIXME
            # Expression: MOI.add_constraint(model, x, set2)
            #   Expected: MathOptInterface.LowerBoundAlreadySet{MathOptInterface.EqualTo{Float64}, MathOptInterface.GreaterThan{Float64}}(MathOptInterface.VariableIndex(1))
            #     Thrown: ErrorException("Cannot put multiple bound sets of the same type on a variable")
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
            # FIXME ArgumentError: MosekTools.Optimizer does not support getting the attribute MathOptInterface.VariablePrimalStart().
            "test_model_VariablePrimalStart",
            # FIXME
            "test_model_duplicate_VariableName",
            # FIXME `MOI.ListOfConstraintAttributesSet` incorrect
            "test_model_ListOfConstraintAttributesSet",
            #  Expression: MOI.set(model, MOI.ConstraintName(), c, "c1")
            #    Expected: MathOptInterface.UnsupportedAttribute{MathOptInterface.ConstraintName}(MathOptInterface.ConstraintName(), "`ConstraintName`s are not supported for `VariableIndex` constraints.")
            #  No exception thrown
            "test_model_VariableIndex_ConstraintName",
            # FIXME Needs https://github.com/jump-dev/MathOptInterface.jl/pull/1787
            "test_variable_solve_ZeroOne_with_0_upper_bound",
        ],
    )
end

@testset "Bridged and cached" begin
    #model = MOIB.full_bridge_optimizer(Mosek.Optimizer(), Float64)
    model = MOIB.full_bridge_optimizer(MosekOptimizerWithFallback(), Float64)
    MOI.set(model, MOI.Silent(), true)

    # linear and basic tests
    MOIT.runtests(model, config,
        exclude=[
            "test_basic_VectorAffineFunction_PositiveSemidefiniteConeSquare", # AssertionError: (m.x_sd[ref2id(vi)]).matrix == -1 src/variable.jl:173
            "test_basic_VectorOfVariables_PositiveSemidefiniteConeSquare",
            "test_basic_VectorAffineFunction_LogDetConeTriangle",
            "test_basic_VectorAffineFunction_NormNuclearCone",
            "test_basic_VectorOfVariables_NormNuclearCone",
            "test_basic_VectorAffineFunction_NormSpectralCone",
            "test_basic_VectorOfVariables_NormSpectralCone",
            "test_basic_VectorAffineFunction_RootDetConeTriangle",
            "test_basic_VectorOfVariables_RootDetConeTriangle",
            "test_basic_VectorAffineFunction_PositiveSemidefiniteConeTriangle", # TODO: implement get ConstraintSet for SAF
            "test_basic_VectorOfVariables_PositiveSemidefiniteConeTriangle",
            "test_conic_LogDetConeTriangle_VectorOfVariables",
            "test_variable_solve_ZeroOne_with_0_upper_bound",
            "test_variable_solve_ZeroOne_with_upper_bound",
            "test_model_ListOfConstraintAttributesSet", # list not properly set
            "BoundAlreadySet", # TODO throw error if bound already set
            "test_model_duplicate_VariableName",
            "test_model_VariablePrimalStart", # able to set but not to get VariablePrimalStart
            "test_objective_set_via_modify",
            # FIXME Mosek.MosekError(1307, "Variable '' (9) is a member of cone '' (0).")
            "test_basic_VectorQuadraticFunction_LogDetConeTriangle",
            "test_basic_VectorOfVariables_LogDetConeTriangle",
            # FIXME Needs https://github.com/jump-dev/MathOptInterface.jl/pull/1787
            "test_constraint_ZeroOne_bounds",
            "test_constraint_ZeroOne_bounds_2",
            "test_constraint_ZeroOne_bounds_3",
            "test_variable_solve_ZeroOne_with_0_upper_bound",
            "test_variable_solve_ZeroOne_with_upper_bound",
        ],
    )

    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)

    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                #Mosek.Optimizer(),
                MosekOptimizerWithFallback(),
        ),
        Float64,
    )
    MOI.set(model, MOI.Silent(), true)

    MOI.Test.runtests(model, config,
        exclude=[
            # FIXME Mosek.MosekError(1307, "Variable '' (1) is a member of cone '' (0).") src/msk_functions.jl:477
            "test_conic_LogDetConeTriangle_VectorOfVariables",
            # FIXME Needs https://github.com/jump-dev/MathOptInterface.jl/pull/1787
            "test_constraint_ZeroOne_bounds",
            "test_variable_solve_ZeroOne_with_0_upper_bound",
            "test_variable_solve_ZeroOne_with_upper_bound",
        ],
    )
end

@testset "Matrix name" begin
    MOI.empty!(optimizer)
    x, cx = MOI.add_constrained_variables(optimizer, MOI.PositiveSemidefiniteConeTriangle(3))
    err = MOI.UnsupportedAttribute{MOI.VariableName}
    @test_throws err MOI.set(optimizer, MOI.VariableName(), x[1], "a")
    MOI.empty!(optimizer)
    model = MOI.Utilities.CachingOptimizer(MOI.Utilities.Model{Float64}(), optimizer)
    x, cx = MOI.add_constrained_variables(model, MOI.PositiveSemidefiniteConeTriangle(3))
    MOI.set(model, MOI.VariableName(), x[1], "a")
    MOI.Utilities.attach_optimizer(model) # Should drop errors silently in the copy
    @test "a" == MOI.get(model, MOI.VariableName(), x[1])
end

# See https://github.com/jump-dev/MosekTools.jl/issues/95
@testset "Mapping enums" begin
    MOI.empty!(optimizer)
    # Force variable bridging to test attribute substitution
    bridged = MOI.Bridges.Variable.Zeros{Float64}(optimizer)
    cache = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    model = MOI.Utilities.CachingOptimizer(cache, bridged)
    attr = MOI.RawOptimizerAttribute("MSK_IPAR_INTPNT_SOLVE_FORM")
    value = MosekTools.MSK_SOLVE_DUAL
    MOI.add_constrained_variables(model, MOI.Zeros(1))
    MOI.Utilities.attach_optimizer(model)
    # The function should not be used so we can use anything here as first argument
    @test MOI.Utilities.substitute_variables(x -> x^2, value) == value
    MOI.set(model, attr, value)
    # Currently Mosek returns `value.value` but we don't want the tests to fail if this gets fixed in the future
    # so we also allow `value`.
    @test MOI.get(model, attr) == value.value || MOI.get(model, attr) == value
    @test MOI.get(optimizer, attr) == value.value || MOI.get(model, attr) == value
end


text:
using Documenter
using Literate

const literate_dir = joinpath(@__DIR__, "literate")
const generated_dir = joinpath(@__DIR__, "src", "generated")

@info "Loading Catlab.jl"
using Catlab

@info "Building Literate.jl docs"
for (root, dirs, files) in walkdir(literate_dir)
  out_dir = joinpath(generated_dir, relpath(root, literate_dir))
  for file in files
    if last(splitext(file)) == ".jl"
      Literate.markdown(joinpath(root, file), out_dir;
        documenter=true, credit=false)
      Literate.notebook(joinpath(root, file), out_dir;
        execute=true, documenter=true, credit=false)
    end
  end
end

@info "Building Documenter.jl docs"
makedocs(
  modules     = [Catlab],
  format      = Documenter.HTML(),
  sitename    = "Catlab.jl",
  doctest     = false,
  pages       = Any[
    "Catlab.jl" => "index.md",
    "Tutorials" => Any[
      "generated/wiring_diagrams/wiring_diagram_basics.md",
      "generated/wiring_diagrams/diagrams_and_expressions.md",
      "generated/graphics/graphviz_wiring_diagrams.md",
      "generated/graphics/tikz_wiring_diagrams.md",
      "generated/programs/algebraic_nets.md",
    ],
    "APIs"      => Any[
      "apis/index.md",
      "apis/catlab.md",
      "apis/doctrines.md",
      "apis/wiring_diagrams.md",
      "apis/graphics.md",
      "apis/programs.md",
    ]
  ]
)

@info "Deploying docs"
deploydocs(
  target = "build",
  repo   = "github.com/epatters/Catlab.jl.git",
  branch = "gh-pages",
  deps   = nothing,
  make   = nothing
)


text:
# +Internal
#
# Print first line of code
#
# For instance:
# #+BEGIN_SRC julia :eval never
# function foo()
# #+END_SRC
#
function org_string_code(di::Documented_Item)::String
    s=raw_string_code(di)
    if !isempty(s)
        s="#+BEGIN_SRC julia :eval never :exports code\n$(s)\n#+END_SRC\n"
    end
    return s
end

# +Internal
#
# Print first line of code AND body
#
# For instance:
# #+BEGIN_SRC julia :eval never
# function foo()
#   # foo body
# end
# #+END_SRC
#
function org_string_code_with_body(di::Documented_Item)::String
    s=raw_string_code_with_body(di)
    if !isempty(s)
        s="#+BEGIN_SRC julia :eval never :exports code\n$(s)\n#+END_SRC\n"
    end
    return s
end


text:
using LinearAlgebra
using ArnoldiMethod: lock!

function example()
    n = 10
    Q = qr(rand(n, n)).Q * Matrix(I, n, n)
    R = diagm(0 => 1:10) .+ triu(randn(n, n) .* 0.01, 1)
    R[2,3] = 1.0
    R[3,2] = -1.0

    # Some Schur decomp
    A = Q * R * Q'

    lock!(R, Q, 1, [8,9,10])
    truncate!()

    @show norm(A * Q - Q * R)

    return R

end


text:
struct Nested <: InferenceType end
const NESTED = Nested()

nested_results(traces, log_weights) = NonparametricSamplingResults{Nested}(NESTED, log_weights, Array{Any, 1}(), traces)

function sample(r :: SamplingResults{Nested}, k, n :: Int)
    v = r[k]
    log_Z = logsumexp(r.log_weights)
    weights = StatsBase.Weights(exp.(r.log_weights .- log_Z))
    if n == 1
        StatsBase.sample(v, weights)
    else
        StatsBase.sample(v, weights, n)
    end
end

@doc raw"""
    function nested(f, replace_fn; params = (), num_points :: Int64 = 1)

Generic implementation of nested sampling (Skilling, Nested sampling for general Bayesian computation, Bayesian Analysis, 2006).
The number of sampling iterations is a function of `num_points` aka ``N``, and the empirical entropy of the sampling distribution, given at the
``n``-th iteration by ``H_n \approx \sum_k \hat p_k \log \hat p_k^{-1}``, where ``\hat p_k = L_k w_k / Z_k``, ``L_k`` is the likelihood
value, ``w_k`` is the difference in prior mass, and ``Z_k`` is the current estimate of the partition function. The number of sampling iterations
is equal to ``\min_n \{n > 0: n > NH_n\}``.

Args:
+ `f`: stochastic function. Must have signature `f(t :: Trace, params...)`
+ `replace_fn`: function that returns a tuple `(new_trace :: Trace, new_log_likelihood :: Float64)`. The input signature of this function
    must be `replace_fn(f :: F, log_likelihood :: Float64, params...) where F <: Function`. It must guarantee that `new_log_likelihood > log_likelihood`.
+ `params`: any parameters to pass to `f`
+ `num_points`: the number of likelihood points to keep track of
"""
function nested(f, replace_fn; params = (), num_points :: Int64 = 1)
    traces = [trace() for _ in 1:num_points]
    for t in traces
        f(t, params...)
    end
    loglikelihoods = [loglikelihood(t) for t in traces]
    weights = Float64[]
    pts = Array{Trace, 1}()
    x_last = 1.0
    Z = 0.0  # partition function estimate
    n = 1  # number of iterations

    while true
        min_l_ix = argmin(loglikelihoods)
        log_l_n = loglikelihoods[min_l_ix]

        x_n = exp(-n / num_points)
        w_n = x_last - x_n
        x_last = x_n

        weight_n = log_l_n + log(w_n)
        Z += exp(weight_n)
        push!(weights, weight_n)
        push!(pts, traces[min_l_ix])

        new_t, this_log_l = replace_fn(f, log_l_n, params...)
        traces[min_l_ix] = new_t
        loglikelihoods[min_l_ix] = this_log_l

        # estimate entropy
        normed_weights = weights .- log(Z)
        if n > 1
            entropy = -1.0 * sum(exp.(normed_weights) .* normed_weights)
            breakpoint = Int(num_points * ceil(entropy))
            if n > breakpoint
                break
            end
        end
        n += 1
    end
    # approximation for the remainder of the sampled traces
    append!(pts, traces)
    x_last = exp(-(n - 1) / num_points)
    x_n = exp(-n / num_points)
    w_n = x_last - x_n
    for t in traces
        push!(weights, loglikelihood(t) + log(w_n))
    end
    nested_results(pts, weights)
end

@doc raw"""
    nested(f; params = (), num_points :: Int64 = 1)

Run nested sampling using internal rejection method.
"""
nested(f; params = (), num_points :: Int64 = 1) = nested(f, rejection; params = params, num_points = num_points)

export nested, sample


text:
if haskey(ENV, "GITHUB_ACTIONS")
    # We're being run as part of a Github action. The most likely case is that
    # this is the auto-merge action as part of the General registry.
    # For now, skip the installation.
    @info("Detected a Github action. Skipping installation.")
    exit(0)
end

using Libdl

const DEPS_FILE = joinpath(@__DIR__, "deps.jl")

if isfile(DEPS_FILE)
    rm(DEPS_FILE)
end

function write_depsfile(path)
    if Sys.iswindows()
        # When `path` gets written out to a file, it will escape any
        # backslashes, so we need to doubly escape them. If your path uses
        # forward slashes, this operation won't do anything.
        path = replace(path, "\\" => "\\\\")
    end
    open(DEPS_FILE, "w") do io
        println(io, "const libgurobi = \"$(path)\"")
    end
end

const ALIASES = [
    "gurobi90",
    "gurobi81", "gurobi80",
    "gurobi75", "gurobi70"
]

paths_to_try = copy(ALIASES)

for a in ALIASES
    if haskey(ENV, "GUROBI_HOME")
        if Sys.isunix()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "lib", string("lib", a, ".so")))
        end
        if Sys.iswindows()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "bin", string(a, ".", Libdl.dlext)))
        end
        if Sys.isapple()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "lib", string("lib", a, ".dylib")))
        end
    end
    # gurobi uses .so on OS X for some reason
    if Sys.isapple()
        push!(paths_to_try, string("lib$a.so"))
        push!(paths_to_try, string("lib$a.dylib"))
    end
end

found = false
for l in paths_to_try
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        global found = true
        write_depsfile(l)
        break
    end
end

function diagnose_gurobi_install()
    println("""
    Unable to locate Gurobi installation. Running some common diagnostics.

    Gurobi.jl only supports the following versions:
    """)
    println.(" - ", ALIASES)
    println("""

    Did you download and install one of these versions from gurobi.com?

    """)
    if haskey(ENV, "GUROBI_HOME")
        dir = joinpath(ENV["GUROBI_HOME"], Sys.isunix() ? "lib" : "bin")
        println("""
        Found GUROBI_HOME =  $(ENV["GUROBI_HOME"])

        Does this point to the correct install location?
        - on Windows, this might be `C:\\Program Files\\gurobi810\\win64\\`
        - alternatively, on Windows, this might be `C:/Program Files/gurobi810/win64/`
        - on OSX, this might be `/Library/gurobi810/mac64/`
        - on Unix, this might be `/home/my_user/gurobi810/linux64/`

        Note: this has to be a full path, not a path relative to your current
        directory or your home directory.

        We're going to look for the Gurobi library in this directory:
            $(dir)

        That directory has the following files:
        """)
        try
            for file in readdir(dir)
                println(" - ", joinpath(dir, file))
            end
            println("""

            We were looking for (but could not find) a file named like
            `libgurobiXXX.so`, `libgurobiXXX.dylib`, or `gurobiXXX.dll`. You
            should update your GUROBI_HOME environment variable to point to the
            correct location.
            """)
        catch ex
            if typeof(ex) <: SystemError
                println("""
                Aha! We tried looking in `$(dir)`, but something went wrong. Are
                you sure that your GUROBI_HOME environment variable is correct?
                When combined with the appropriate suffix (e.g., `lib` or
                `bin`, it needs to point to a valid directory.
                """)
            else
                rethrow(ex)
            end
        end
    else
        try
            println("""
            Looking for a version of Gurobi in your path:
            """)
            # Try to call `gurobi_cl`. This should work if Gurobi is on the
            # system path. If it succeeds, it will print out the version.
            run(`gurobi_cl --version`)
            println("""

            We couldn't find the `GUROBI_HOME` environment variable, but we
            found this version of Gurobi on your path. Is this version one of
            the supported versions listed above? If not, you should edit your
            `PATH` to point to the correct version.
            """)
        catch
            println("""

            We could not find a version of Gurobi in your path, and we could
            not find the environment variable `GUROBI_HOME`. You should set
            the `GUROBI_HOME` environment variable to point to the install
            location. For example:
            - on Windows, this might be `C:\\Program Files\\gurobi810\\win64\\`
            - on OSX, this might be `/Library/gurobi810/mac64/`
            - on Unix, this might be `/opt/gurobi810/linux64/`

            Alternatively, you can add the Gurobi install directory to your path.
            """)
        end
    end
end

if !found && !haskey(ENV, "GUROBI_JL_SKIP_LIB_CHECK")
    diagnose_gurobi_install()
    error("""
    Unable to locate Gurobi installation. If the advice above did not help,
    open an issue at https://github.com/JuliaOpt/Gurobi.jl and post the full
    print-out of this diagnostic attempt.
    """)
end