text:
"""
$(TYPEDEF)

The abstract type describing a control system. With state-space dimension `nx`,
control dimension `nu` and sampling time ΔT. By convention, for `ΔT == 0` the
system is continuous.

"""
abstract type ControlSystem{ΔT, nx, nu} end

"""
    $(TYPEDSIGNATURES)

Returns the sampling time ΔT of the system.
"""
samplingtime(::Type{<:ControlSystem{ΔT}}) where {ΔT} = ΔT
samplingtime(cs::ControlSystem) = samplingtime(typeof(cs))

"""
    $(TYPEDSIGNATURES)

Returns true if the system is has a non-zero sampling rate (e.g. discrete or sampled
contiuous system.
"""
issampled(cs::ControlSystem) = !iszero(samplingtime(cs))

"""
    $(TYPEDSIGNATURES)

Returns the number of states of the control system.
"""
n_states(::Type{<:ControlSystem{ΔT, nx}}) where {ΔT, nx} = nx
n_states(cs::ControlSystem) = n_states(typeof(cs))

"""
    $(TYPEDSIGNATURES)

Returns the number of controls of the control system.
"""
n_controls(::Type{<:ControlSystem{ΔT, nx, nu}}) where {ΔT, nx, nu} = nu
n_controls(cs::ControlSystem) = n_controls(typeof(cs))

"-------------------------------- Minimal Interface -------------------------------"

"""
    $(FUNCTIONNAME)(cs::ControlSystem, x::SVector, u::SVector, t::AbstractFloat)

Returns the time derivative of the state `dx` at a given state `x`, control
input `u` and time `t`
"""
function dx end

""""
    $(FUNCTIONNAME)(cs::ControlSystem, xₖ::SVector, uₖ::SVector, k::Int)

Returns the next state (`xₖ₊₁`) for a control system with non-zero `ΔT` when
applying input `uₖ` at state `xₖ` and time step k

We provide a convencience default below.
"""
function next_x end

"""
    $(FUNCTIONNAME)(cs::ControlSystem, x::SVector, u::SVector, t::AbstractFloat)

Returns the continuous time Jacobian linearization of the dynamics defined in
`dx` at a given state `x`, control `u` and time point `t` in terms of a
`LinearSystem` (A, B).  (such that Δẋ ≈ A*Δx + B*Δu)

We provide a convencience default below.
"""
function linearize end

""""
    $(FUNCTIONNAME)(cs::ControlSystem, x::SVector, u::SVector, t::AbstractFloat)

Returns the discrete time (ZoH) Jacobian linearization for a system with
non-zero sampling time ΔT. Input arguments have the same meaning as for
`linearization`.

We provide a convencience default below.
"""
function linearize_discrete end

"""
    $(FUNCTIONNAME)(cs::ControlSystem, x0::SVector, u::SVector, t0::AbstractFloat,
                    ΔT::AbstractFloat, n_intsteps::Float64)

Integrate propagate the state `x0` given the control `u` and the current time
`t0` until `t0+ΔT`. Here, `n_intsteps` is the number of steps that is used between
`t0` and `t0+ΔT`.

We provide a convencience default below.
"""
function integrate end

"--------------------- Optional Interface Extensions ---------------------"

"""
    $(FUNCTIONNAME)(cs::ControlSystem)

Returns the indeces of the x and y coordinate of this sytem as SIndex.
"""
function xyindex end

"""
    $(FUNCTIONNAME)(cs::ControlSystem)

Returns the indices of the states of this system as SIndex.
"""
function xindex end

"--------------------- Convencience Implementations ---------------------"

"""
    $(TYPEDSIGNATURES)

A convencience implementation of `linearize` using `ForwardDiff.jl`. Overload
this with an explicit version to get better performance.
"""
function _linearize_ad(cs::ControlSystem, x::SVector, u::SVector, t::AbstractFloat)
    @warn "You are using the fallback linearization using ForwardDiff. Consider
    implementing a custom `linearize` or `linearize_discrete` for your
    `ControlSystem` type." maxlog=1
    A = ForwardDiff.jacobian(x->dx(cs, x, u, t), x)
    B = ForwardDiff.jacobian(u->dx(cs, x, u, t), u)
    return LinearSystem{0}(A, B)
end
function linearize(cs::ControlSystem, x::SVector, u::SVector, t::AbstractFloat)
    return _linearize_ad(cs, x, u, t)
end

"""
    $(TYPEDSIGNATURES)

A convencience implementation of `linearize_discrete`. For better performance,
this may be overloaded with some explicit (analytic) expressions (that may even
avoid calling `linearize`).
"""
linearize_discrete(cs::ControlSystem,
                   x0::SVector,
                   u0::SVector,
                   t0::AbstractFloat) = discretize(linearize(cs, x0, u0, t0),
                                                   Val{samplingtime(cs)}())

"""
    $(TYPEDSIGNATURES)

A convencience implementation of `integrate` using RungeKutta of order 4 as
convencience default. You may impelement your own integrator here. Consider
using `DifferentialEquations.jl.
"""
function integrate(cs::ControlSystem, x0::SVector, u::SVector, t0::AbstractFloat,
                   ΔT::AbstractFloat, n_intsteps::Int=2)
    Δt = ΔT/n_intsteps
    x = x0
    for t in range(t0, stop=t0+ΔT, length=n_intsteps+1)[1:end-1]
        k1 = Δt * dx(cs, x, u, t);
        k2 = Δt * dx(cs, x + 0.5 * k1, u, t + 0.5 * Δt);
        k3 = Δt * dx(cs, x + 0.5 * k2, u, t + 0.5 * Δt);
        k4 = Δt * dx(cs, x + k3      , u, t + Δt);
        x += (k1 + 2.0 * (k2 + k3) + k4) / 6.0;
    end
    return x
end

"""
    $(TYPEDSIGNATURES)

Integrate to xₖ₊₁ starting from x, applying u.
"""
function next_x(cs::ControlSystem, x::SVector, u::SVector, t::Float64)
    @assert issampled(cs) "next_x requires `ControlSystem` with ΔT > 0."
    return integrate(cs, x, u, t, samplingtime(cs))
end


text:
function create(::Type{HarmonicCurveFitter})
    return jcall(HarmonicCurveFitter, "create", HarmonicCurveFitter, ())
end

function with_max_iterations(obj::HarmonicCurveFitter, arg0::jint)
    return jcall(obj, "withMaxIterations", HarmonicCurveFitter, (jint,), arg0)
end

function with_start_point(obj::HarmonicCurveFitter, arg0::Vector{jdouble})
    return jcall(obj, "withStartPoint", HarmonicCurveFitter, (Vector{jdouble},), arg0)
end

ext: jl


text:
################################################################################
#
#   METHODS FOR CONSTRUCTION INTERACTION MATRICES FOR LATTICES
#
#   STRUCTURE OF THE FILE
#
#   1) INTERACTION MATRICES IN REAL SPACE
#
#   2) INTERACTION MATRICES IN MOMENTUM SPACE
#
#   NOTE: No Majorana fermions so far on the level of matrices, because of some
#   problems with the various gauges that can be tuned. Want: sublattice
#   idenfitication in bipartite systems.
#
################################################################################

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Interaction matrix in REAL space
#
#   Parameters are
#   - lattice: The complete lattice of which the interaction matrix should be constructed
#   - enforce_hermitian (optional): if the matrix should be made hermitian by 0.5*(A + A_dagger)
#
#-----------------------------------------------------------------------------------------------------------------------------
"""
    getInteractionMatrixRealSpace(unitcell::Unitcell  [ ; enforce_hermitian::Bool=false ])
    getInteractionMatrixRealSpace(lattice::Lattice    [ ; enforce_hermitian::Bool=false ])

Constructs the interaction matrix (in *position* space) of a given `Unitcell` or `Lattice` object.
The matrix is a NxN matrix where N is the number of sites in the given object.
Entries (i,j) contain the interactions between sites i and j.

The matrix is of type `Array{Complex,2}`, i.e. it has complex entries to account for eventual complex hopping amplitudes.

Furthermore, the matrix can be constructed hermitian by using `enforce_hermitian=true` in which case
the interactions between i and j are not only added to element (i,j) of the matrix but also to element (j,i) as conjugated.

# Examples

```julia-repl
julia> getInteractionMatrixRealSpace(lattice)
2×2 Array{Complex,2}:
...

julia> getInteractionMatrixRealSpace(unitcell, enforce_hermitian=true)
2×2 Array{Complex,2}:
...

```
"""
function getInteractionMatrixRealSpace(lattice::Lattice; enforce_hermitian::Bool=false)
    # generate a new matrix
    matrix = zeros(Complex, length(lattice.positions),length(lattice.positions))
    # decide if hermitian or not hermitian enforced
    if enforce_hermitian
        # iterate over all connections and add strength to two entries
        for c in lattice.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            # add to the matrix twice
            matrix[index_from, index_to] += 0.5 * strength
            matrix[index_to, index_from] += 0.5 * conj(strength)
        end
    else
        # iterate over all connections and add strength to one entry
        for c in lattice.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            # add to the matrix once
            matrix[index_from, index_to] += strength
        end
    end
    # return the matrix
    return matrix
end
function getInteractionMatrixRealSpace(unitcell::Unitcell; enforce_hermitian::Bool=false)
    # generate a new matrix
    matrix = zeros(Complex, length(unitcell.basis),length(unitcell.basis))
    # decide if hermitian or not hermitian enforced
    if enforce_hermitian
        # iterate over all connections and add strength to two entries
        for c in unitcell.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            # add to the matrix twice
            matrix[index_from, index_to] += 0.5 * strength
            matrix[index_to, index_from] += 0.5 * conj(strength)
        end
    else
        # iterate over all connections and add strength to one entry
        for c in unitcell.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            # add to the matrix once
            matrix[index_from, index_to] += strength
        end
    end
    # return the matrix
    return matrix
end
export getInteractionMatrixRealSpace

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Interaction matrix in RECIPROCAL (K) space
#
#   Parameters are
#   - lattice: The complete lattice of which the interaction matrix should be constructed
#   - k_vector: The reciprocal point k at which the matrix is constructed
#   - enforce_hermitian (optional): if the matrix should be made hermitian by 0.5*(A + A_dagger)
#
#-----------------------------------------------------------------------------------------------------------------------------
"""
    getInteractionMatrixKSpace(unitcell::Unitcell, k_vector::Array{Float64,1} [ ; enforce_hermitian::Bool=false ])
    getInteractionMatrixKSpace(lattice::Lattice,   k_vector::Array{Float64,1} [ ; enforce_hermitian::Bool=false ])

Constructs the interaction matrix (in *momentum* space at point `k_vector`) of a given `Unitcell` or `Lattice` object.
The matrix is a NxN matrix where N is the number of sites in the given object.
Entries (i,j) contain the interactions between sites i and j as well as the phase factor exp(i k*delta).

The matrix is of type `Array{Complex,2}`, i.e. it has complex entries to account for the phase factor.

Furthermore, the matrix can be constructed hermitian by using `enforce_hermitian=true` in which case
the interactions between i and j are not only added to element (i,j) of the matrix but also to element (j,i) as conjugated.

# Examples

```julia-repl
julia> getInteractionMatrixKSpace(lattice, k)
2×2 Array{Complex,2}:
...

julia> getInteractionMatrixKSpace(lattice, [pi/2.0, 0.0], enforce_hermitian=true)
2×2 Array{Complex,2}:
...

```
"""
function getInteractionMatrixKSpace(lattice::Lattice, k_vector::Array{Float64,1}; enforce_hermitian::Bool=false)
    # generate a new matrix
    matrix = zeros(Complex, size(lattice.positions,1),size(lattice.positions,1))
    # decide if hermitian or not hermitian enforced
    if enforce_hermitian
        # iterate over all connections
        for c in lattice.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            wrap        = c[4]
            # get the difference vector
            pos_delta   = lattice.positions[index_to] .- lattice.positions[index_from]
            if size(lattice.lattice_vectors,1) > 0
                for pair in zip(wrap, lattice.lattice_vectors)
                    pos_delta .+= pair[1].*pair[2]
                end
            end
            # add to the matrix twice
            matrix[index_from, index_to] += 0.5 * strength      * exp(-im * sum(pos_delta .* k_vector))
            matrix[index_to, index_from] += 0.5 * conj(strength * exp(-im * sum(pos_delta .* k_vector)))
        end
    else
        # iterate over all connections
        for c in lattice.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            wrap        = c[4]
            # get the difference vector
            pos_delta   = lattice.positions[index_to] .- lattice.positions[index_from]
            if size(lattice.lattice_vectors,1) > 0
                for pair in zip(wrap, lattice.lattice_vectors)
                    pos_delta .+= pair[1].*pair[2]
                end
            end
            # add to the matrix once
            matrix[index_from, index_to] += strength * exp(-im * sum(pos_delta .* k_vector))
        end
    end
    # return the matrix
    return matrix
end
function getInteractionMatrixKSpace(unitcell::Unitcell, k_vector::Array{Float64,1}; enforce_hermitian::Bool=false)
    # generate a new matrix
    matrix = zeros(Complex, size(unitcell.basis,1),size(unitcell.basis,1))
    # decide if hermitian or not hermitian enforced
    if enforce_hermitian
        # iterate over all connections
        for c in unitcell.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            wrap        = c[4]
            # get the difference vector
            pos_delta   = unitcell.basis[index_to] .- unitcell.basis[index_from]
            if size(unitcell.lattice_vectors,1) > 0
                for pair in zip(wrap, unitcell.lattice_vectors)
                    pos_delta .+= pair[1].*pair[2]
                end
            end
            # add to the matrix twice
            matrix[index_from, index_to] += 0.5 * strength      * exp(-im * sum(pos_delta .* k_vector))
            matrix[index_to, index_from] += 0.5 * conj(strength * exp(-im * sum(pos_delta .* k_vector)))
        end
    else
        # iterate over all connections
        for c in unitcell.connections
            # get the indices
            index_from  = Int(c[1])
            index_to    = Int(c[2])
            strength    = c[3]
            wrap        = c[4]
            # get the difference vector
            pos_delta   = unitcell.basis[index_to] .- unitcell.basis[index_from]
            if size(unitcell.lattice_vectors,1) > 0
                for pair in zip(wrap, unitcell.lattice_vectors)
                    pos_delta .+= pair[1].*pair[2]
                end
            end
            # add to the matrix once
            matrix[index_from, index_to] += strength * exp(-im * sum(pos_delta .* k_vector))
        end
    end
    # return the matrix
    return matrix
end
export getInteractionMatrixKSpace


text:
#=================================================================================
  NWI (Nocedal Wright Iterative)) is a trust region sub-problem solver that sol-
  ves the quadratic programming problem subject to the solution being in a
  Euclidean Ball. It is appropriate for all types of Hessians. It solves the
  problem by solving the secular equ  ation by a safeguarded Newton's method.
  The secular equation is simply a the inverse of the difference between the
  step length and the trust region radius. It accepts either the Newton step
  if this is interior and the Hessian is positive definite, or steps to the
  boundary (with some slack) in a clever way.

  NTR accepts all Hessians and is therefore well-suited for Newton's method for
  any vexity. It is expensive for very large problems, and uses the direct
  eigensolution. This could potentially be useful for problems with simple
  structure of the eigenproblem, but I suspect NWI is superior in most cases.
=================================================================================#

struct NWI <: NearlyExactTRSP
end
summary(::NWI) = "Trust Region (Newton, eigen)"

"""
    initial_safeguards(B, h, g, Δ)

Returns a tuple of initial safeguarding values for λ. Newton's method might not
work well without these safeguards when the Hessian is not positive definite.
"""
function initial_safeguards(B, h, g, Δ)
    # equations are on p. 560 of [MORESORENSEN]
    T = eltype(g)
    λS = maximum(-h)

    # they state on the first page that ||⋅|| is the Euclidean norm
    gnorm = norm(g)
    Bnorm = opnorm(B, 1)
    λL = max(T(0), λS, gnorm/Δ - Bnorm)
    λU = gnorm/Δ + Bnorm
    (L=λL, U=λU, S=λS)
end
function safeguard_λ(λ::T, λsg) where T
    # p. 558
    λ = min(max(λ, λsg.L), λsg.U)
    if λ ≤ λsg.S
        λ = max(T(1)/1000*λsg.U, sqrt(λsg.L*λsg.U))
    end
    λ
end

"""
    is_maybe_hard_case(QΛQ, Qt∇f)

Returns a tuple of a Bool, `hardcase` and an integer, `λidx`. `hardcase` is true
if the sub-problem is the "hard case", that is the case where the smallest
eigenvalue is negative. `λidx` is the index of the first eigenvalue not equal
to the smallest eigenvalue. `QΛQ.values` holds the eigenvalues of H sorted low
to high, `Qt∇f` is a vector of the inner products between the eigenvectors and the
gradient.
"""
function is_maybe_hard_case(QΛQ, Qt∇f::AbstractVector{T}) where T
    # If the solution to the trust region sub-problem is on the boundary of the
    # trust region, {w : ||w|| ≤ Δk}, then the solution is usually found by
    # finding a λ ≥ 0 such that ||(B + λI)⁻¹g|| = Δ and x'(B + λI)x > 0, x≠0.
    # However, in the hard case our strategy for iteratively finding such an
    # λ breaks down, because ||p(λ)|| does not go to ∞ as λ→-λⱼ for j such that
    # qₜ'g ≠ 0 (see section 4.3 in [N&W] for more details). Remember, that B
    # does *not* have to be pos. def. in the trust region based methods!

    # Get the eigenvalues
    Λ = QΛQ.values

    # Get number of eigen values
    λnum = length(Λ)

    # The hard case requires a negativ smallest eigenvalue.
    λmin = first(Λ) # eigenvalues are sorted
    if λmin >= T(0)
        return false, 1
    end

    # Assume hard case and verify
    λidx = 1
    hard_case = true
    for (Qt∇f_j, λ_j) in zip(Qt∇f, Λ)
        if abs(λmin - λ_j) > sqrt(eps(T))
            hard_case = true
            break
        else
            if abs(Qt∇f_j) > sqrt(eps(T))
                hard_case = false
                break
            end
        end
        λidx += 1
    end

    hard_case, λidx
end

# Equation 4.38 in N&W (2006)
calc_p!(p, Qt∇f, QΛQ, λ) = calc_p!(p, Qt∇f, QΛQ, λ, 1)

# Equation 4.45 in N&W (2006) since we allow for first_j > 1
function calc_p!(p, Qt∇f, QΛQ, λ::T, first_j) where T
    # Reset search direction to 0
    fill!(p, T(0))

    # Unpack eigenvalues and eigenvectors
    Λ = QΛQ.values
    Q = QΛQ.vectors
    for j = first_j:length(Λ)
        κ = Qt∇f[j] / (Λ[j] + λ)
        @. p = p - κ*Q[:, j]
    end
    p
end

"""
    solve_tr_subproblem!(∇f, H, Δ, s; abstol, maxiter)
Args:
    ∇f: The gradient
    H:  The Hessian
    Δ:  The trust region size, ||s|| <= Δ
    s: Memory allocated for the step size, updated in place
    abstol: The convergence abstol for root finding
    maxiter: The maximum number of root finding iterations

Returns:
    m - The numeric value of the quadratic minimization.
    interior - A boolean indicating whether the solution was interior
    lambda - The chosen regularizing quantity
    hard_case - Whether or not it was a "hard case" as described by N&W (2006)
    solved - Whether or not a solution was reached (as opposed to
      terminating early due to maxiter)
"""
function (ms::NWI)(∇f, H, Δ, p, scheme; abstol=1e-10, maxiter=50)
    T = eltype(p)
    n = length(∇f)
    H = H isa UniformScaling ? Diagonal(copy(∇f).*0 .+ 1) : H
    h = diag(H)

    # Note that currently the eigenvalues are only sorted if H is perfectly
    # symmetric.  (Julia issue #17093)
    if H isa Diagonal
        QΛQ = eigen(H)
    else
        QΛQ = eigen(Symmetric(H))
    end
    Q, Λ = QΛQ.vectors, QΛQ.values
    λmin, λmax = Λ[1], Λ[n]
    # Cache the inner products between the eigenvectors and the gradient.
    Qt∇f = Q' * ∇f

    # These values describe the outcome of the subproblem.  They will be
    # set below and returned at the end.
    solved = true

    # Potentially an unconstrained/interior solution. The smallest eigenvalue is
    # positive, so the Newton step, pN, is fine unless norm(pN, 2) > Δ.
    if λmin >= sqrt(eps(T))
        λ = T(0) # no amount of I is added yet
        p = calc_p!(p, Qt∇f, QΛQ, λ) # calculate the Newton step
        if norm(p, 2) ≤ Δ
            # No shrinkage is necessary: -(H \ ∇f) is the minimizer
            interior = true
            solved = true
            hard_case = false

            m = dot(∇f, p) + dot(p, H * p)/2

            return (p=p, mz=m, interior=interior, λ=λ, hard_case=hard_case, solved=solved, Δ=Δ)
        end
    end

    # Set interior flag
    interior = false

    # The hard case is when the gradient is orthogonal to all
    # eigenvectors associated with the lowest eigenvalue.
    maybe_hard_case, first_j = is_maybe_hard_case(QΛQ, Qt∇f)
    # Solutions smaller than this lower bound on lambda are not allowed:
    # as the shifted Hessian will not be PSD.
    λ_lb = -λmin
    λ = λ_lb

    # Verify that it is actually the hard case situation by calculating the
    # step with λ = λmin (it's currently λ_lb, verify that that is correct).
    if maybe_hard_case
        # The "hard case".
        # λ is taken to be -λmin and we only need to find a multiple of an
        # orthogonal eigenvector that lands the iterate on the boundary.

        # The old p is discarded, and replaced with one that takes into account
        # the first j such that λj ≠ λmin. Formula 4.45 in N&W (2006)
        pλ = calc_p!(p, Qt∇f, QΛQ, λ, first_j)

        # Check if the choice of λ leads to a solution inside the trust region.
        # If it does, then we construct the "hard case solution".
        if norm(pλ, 2) ≤ Δ
            hard_case = true
            solved = true

            tau = sqrt(Δ^2 - norm(pλ, 2)^2)

            @. p = -pλ + tau * Q[:, 1]

            m = dot(∇f, p) + dot(p, H * p)/2

            return (p=p, mz=m, interior=interior, λ=λ, hard_case=hard_case, solved=solved, Δ=Δ)
        end
    end
    # If this is reached, we cannot be in the hard case after all, and we
    # can use Newton's method to find a p such that norm(p, 2) = Δ.
    hard_case = false

    # Algorithim 4.3 of N&W (2006), with s insted of p_l for consistency
    # with Optim.jl

    solved = false
    # We cannot be in the λ = -λ₁ case, as the root is in (-λ₁, ∞) interval.
    λ = λ + sqrt(eps(T))
    isg = initial_safeguards(H, h, ∇f, Δ)
    λ = safeguard_λ(λ, isg)
    for iter in 1:maxiter
        λ_previous = λ
        for i = 1:n
            @inbounds H[i, i] = h[i] + λ
        end

        F = H isa Diagonal ? cholesky(H; check=false) : cholesky(Hermitian(H); check=false)
        if !issuccess(F)
            # We should not be here, but the lower bound might fail for
            # numerical reasons.
            if λ < λ_lb
                λ = T(1)/2 * (λ_previous - λ_lb) + λ_lb
            end
            continue
        end
        R = F.U
        p .= R \ (R' \ -∇f)
        q_l = R' \ p

        p_norm = norm(p, 2)
        λ_update = p_norm^2 * (p_norm - Δ) / (Δ * dot(q_l, q_l))
        λ += λ_update

        # Check that λ is not less than λ_lb, and if so, go
        # half the way to λ_lb. (This should be geometric mean)
        if λ < λ_lb
            λ = T(1)/2 * (λ_previous - λ_lb) + λ_lb
        end
        if abs(λ - λ_previous) ≤ abstol
            solved = true
            break
        end
    end
    for i = 1:n
        @inbounds H[i, i] = h[i]
    end
    m = dot(∇f, p) + dot(p, H * p)/2
    return (p=p, mz=m, interior=interior, λ=λ, hard_case=hard_case, solved=solved, Δ=Δ)
end


text:
module SplittingOperators

export @Lie, @Strang, @TripleJump, @Order6, @Magic, @SuperMagic

"""

    @Lie( push_t, push_v )

    Apply the first order Lie splitting

    push_t and push_v are two function calls with
    `dt` as argument.

"""
macro Lie(push_t, push_v)
    return esc( quote
        $push_t
	$push_v
    end)
end

"""

    @Strang( push_t, push_v )

    Apply the second order Strang splitting

    push_t and push_v are two function calls with
    `dt` as argument.

"""
macro Strang(push_t, push_v)
    return esc(quote
        local full_dt = dt
	dt = 0.5full_dt
        $push_t
	dt = full_dt
        $push_v
	dt = 0.5full_dt
        $push_t
	dt = full_dt
    end)
end

"""

    @Magic( push_t, push_v )

    Apply the second order modified Strang splitting

    with tan(dt/2) instead of dt/2 for push_t and sin(dt)
    instead of dt for push_v

    push_t and push_v are two function calls with
    `dt` as argument.

"""
macro Magic(push_t, push_v)
    return esc(quote
        local full_dt = dt
	  dt = tan(0.5full_dt)
          $push_t
	  dt = sin(full_dt)
          $push_v
	  dt = tan(0.5full_dt)
          $push_t
	  dt = full_dt
    end)
end

"""

    @TripleJump( push_t, push_v )

    Apply the fourth order Triple Jump splitting

    push_t and push_v are two function calls with
    `dt` as argument.

"""
macro TripleJump(push_t, push_v)
    return esc(quote
        local full_dt = dt
        dt =  0.675603595979829full_dt
        $push_t
        dt =  1.351207191959658full_dt
        $push_v
        dt = -0.17560359597982855full_dt
        $push_t
        dt = -1.702414383919315full_dt
        $push_v
        dt = -0.17560359597982855full_dt
        $push_t
        dt =  1.351207191959658full_dt
        $push_v
        dt =  0.675603595979829full_dt
        $push_t
	dt = full_dt
    end)
end

"""

    @Order6( push_t, push_v )

    Apply the sixth order splitting

    push_t and push_v are two function calls with
    `dt` as argument.

"""
macro Order6(push_t, push_v)
    return esc(quote
        local full_dt = dt
        dt =  0.0414649985182624full_dt
        $push_t
        dt =  0.123229775946271full_dt
        $push_v
        dt =  0.198128671918067full_dt
        $push_t
        dt =  0.290553797799558full_dt
        $push_v
        dt = -0.0400061921041533full_dt
        $push_t
        dt = -0.127049212625417full_dt
        $push_v
        dt =  0.0752539843015807full_dt
        $push_t
        dt = -0.246331761062075full_dt
        $push_v
        dt = -0.0115113874206879full_dt
        $push_t
        dt =  0.357208872795928full_dt
        $push_v
        dt =  0.23666992478693111full_dt
        $push_t
        dt =  0.20477705429147008full_dt
        $push_v
        dt =  0.23666992478693111full_dt
        $push_t
        dt =  0.357208872795928full_dt
        $push_v
        dt = -0.0115113874206879full_dt
        $push_t
        dt = -0.246331761062075full_dt
        $push_v
        dt =  0.0752539843015807full_dt
        $push_t
        dt = -0.127049212625417full_dt
        $push_v
        dt = -0.0400061921041533full_dt
        $push_t
        dt =  0.290553797799558full_dt
        $push_v
        dt =  0.198128671918067full_dt
        $push_t
        dt =  0.123229775946271full_dt
        $push_v
        dt =  0.0414649985182624full_dt
        $push_t
	dt = full_dt
    end)
end

end


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavfilter/asrc_abuffer.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_asrc_buffer_add_samples,
    av_asrc_buffer_add_buffer,
    av_asrc_buffer_add_audio_buffer_ref

function av_asrc_buffer_add_samples(abuffersrc,data,linesize,nb_samples::Integer,sample_rate::Integer,sample_fmt::Integer,ch_layout::Int64,planar::Integer,pts::Int64,flags::Integer)
    ccall((:av_asrc_buffer_add_samples,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{Ptr{Uint8}},Ptr{Cint},Cint,Cint,Cint,Int64,Cint,Int64,Cint),abuffersrc,data,linesize,nb_samples,sample_rate,sample_fmt,ch_layout,planar,pts,flags)
end

function av_asrc_buffer_add_buffer(abuffersrc,buf,buf_size::Integer,sample_rate::Integer,sample_fmt::Integer,ch_layout::Int64,planar::Integer,pts::Int64,flags::Integer)
    ccall((:av_asrc_buffer_add_buffer,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{Uint8},Cint,Cint,Cint,Int64,Cint,Int64,Cint),abuffersrc,buf,buf_size,sample_rate,sample_fmt,ch_layout,planar,pts,flags)
end

function av_asrc_buffer_add_audio_buffer_ref(abuffersrc,samplesref,flags::Integer)
    ccall((:av_asrc_buffer_add_audio_buffer_ref,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{AVFilterBufferRef},Cint),abuffersrc,samplesref,flags)
end


text:
include("GroverSearch.jl")
using Yao.Intrinsics
# if basis[abs(locs)] == locs>0?1:0, then flip the sign.
inference_oracle(locs) = control(locs[1:end-1], abs(locs[end])=>Z)

################################################
#                Doing Inference               #
################################################
target_space(oracle) = real(Diagonal(mat(oracle)).diag) .< 0

function inference1(psi::AbstractRegister, evidense::Vector{Int}, num_iter::Int)
    oracle = inference_oracle(evidense)(nqubits(psi))
    ts = target_space(oracle)
    rv1 = Reflect(copy(psi))
    grover = chain(oracle, rv1)
    for i in 1:num_iter
        p_target = norm(statevec(psi)[ts])^2
        println("step $i, overlap = $p_target")
        grover(psi)
    end
    psi
end

function run_inference(psi, evidense)
    # the desired subspace
    subinds = indices_with(nqubits(psi), abs.(evidense), Int.(evidense.>0))

    v_desired = statevec(psi)[subinds+1]
    p = norm(v_desired)^2
    v_desired[:] ./= sqrt(p)

    # search the subspace
    num_iter = num_grover_step(p)
    println("Estimated num_step = ", pi/4/sqrt(p))
    inference1(psi, evidense, num_iter)
    fidelity = (psi.state[subinds+1]'*v_desired) |> abs
    fidelity
end

# the second version

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
# Clusters
# ========

# Clustering.jl : ClusteringResult

"Use `NoClustering()` to avoid the use of clustering where a `Clusters` type is needed."
struct NoClustering <: ClusteringResult end

"Data structure to represent sequence clusters. The sequence data itself is not included."
@auto_hash_equals struct Clusters <: ClusteringResult
    clustersize::Vector{Int}
    clusters::Vector{Int}
    weights::StatsBase.Weights{Float64, Float64, Array{Float64,1}}
end

nelements(cl::Clusters) = length(cl.clusters)

# Clustering.jl : nclusters, counts, assignments

"Get the number of clusters in a `Clusters` object."
Clustering.nclusters(cl::Clusters) = length(cl.clustersize)

"""
Get sample counts of clusters as a `Vector`. Each `k` value is the number of samples
assigned to the k-th cluster.
"""
Clustering.counts(cl::Clusters) = cl.clustersize

"""
Get a vector of assignments, where the `i` value is the index/number of the cluster to
which the i-th sequence is assigned.
"""
Clustering.assignments(cl::Clusters) = cl.clusters

#  convert from ClusteringResult of Clustering.jl to Clusters
# -------------------------------------------------------------------

# needs tests
function Base.convert(::Type{Clusters}, cl::ClusteringResult)
    clustersize = counts(cl)
    clusters = assignments(cl)
    weights = Weights(Float64[1.0/clustersize[k] for k in clusters],
        Float64(length(clustersize)))
    Clusters(clustersize, clusters, weights)
end

@inline Base.convert(::Type{Clusters}, cl::Clusters) = cl # no-op

# weights
# -------

"""
`getweight(c[, i::Int])`

This function returns the weight of the sequence number `i`. getweight should be defined for
any type used for `count!`/`count` in order to use his weigths. If `i` isn't used, this
function returns a vector with the weight of each sequence.
"""
@inline getweight(weight::NoClustering, seq::Int) = 1.0

getweight(cl::Clusters) = cl.weights

getweight(cl::Clusters, seq::Int) = cl.weights[seq]

@inline getweight(cl::Weights, i::Int) = cl[i]

# getweight(cl::ClusteringResult) = getweight(convert(Clusters, cl))
# getweight(cl::ClusteringResult, seq::Int) = getweight(convert(Clusters, cl), seq)


text:
mutable struct Entity
    genes:: Array{Int64,1}
    fitness:: Int64
end

function generateEntity(genesLength)
    return Entity(rand(0:1,genesLength), 0)
end

function calculateFitness!(entity, fitValue)
    entity.fitness = 0
    for i in 1:length(entity.genes)
        entity.fitness += entity.genes[i]
    end
    entity.fitness = abs(entity.fitness - fitValue)
end

function printEntity(entity)
    for i in 1:length(entity.genes)
        print(entity.genes[i])
    end
    print(" $(entity.fitness) \n");
end

function crossover!(entity1,entity2,crossoverPoint)
    for i in 1:crossoverPoint
        x = entity1.genes[i]
        entity1.genes[i] = entity2.genes[i]
        entity2.genes[i] = x
    end
end

function mutate!(entity,mutatePercentage)
    if rand(Float64) < mutatePercentage
        mutationPoint = rand(1:genesLength)
        entity.genes[mutationPoint] = 1 - entity.genes[mutationPoint]
    end
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
struct MergeMapper
    cellmaps::Any
    function MergeMapper()
        south = cell_map_to_south()
        east = cell_map_to_east()
        north = cell_map_to_north()
        west = cell_map_to_west()
        southeast = cell_map_to_south_east()
        northeast = cell_map_to_north_east()
        northwest = cell_map_to_north_west()
        southwest = cell_map_to_south_west()
        new([south, east, north, west, southeast, northeast, northwest, southwest])
    end
end

function Base.getindex(m::MergeMapper, i)
    return m.cellmaps[i]
end

function cell_map_to_south()
    cellmap = CellMap([-1.0, -3.0], [1.0, -1.0])
    return cellmap
end

function cell_map_to_east()
    cellmap = CellMap([1.0, -1.0], [3.0, 1.0])
    return cellmap
end

function cell_map_to_north()
    cellmap = CellMap([-1.0, 1.0], [1.0, 3.0])
    return cellmap
end

function cell_map_to_west()
    cellmap = CellMap([-3.0, -1.0], [-1.0, 1.0])
    return cellmap
end

function cell_map_to_south_east()
    cellmap = CellMap([1.0, -3.0], [3.0, -1.0])
    return cellmap
end

function cell_map_to_north_east()
    cellmap = CellMap([1.0, 1.0], [3.0, 3.0])
    return cellmap
end

function cell_map_to_north_west()
    cellmap = CellMap([-3.0, 1.0], [-1.0, 3.0])
    return cellmap
end

function cell_map_to_south_west()
    cellmap = CellMap([-3.0, -3.0], [-1.0, -1.0])
    return cellmap
end

function map_quadrature(quad, mergemapper, mapid)
    mappedpoints = mergemapper[mapid](quad.points)
    return QuadratureRule(mappedpoints, quad.weights)
end

function update_quadrature!(cellquads::CellQuadratures, s, cellid, quad)
    row = cell_sign_to_row(s)
    idx = cellquads.celltoquad[row, cellid]
    cellquads.quads[idx] = quad
end

function update_quadrature!(interfacequads::InterfaceQuadratures, s, cellid, quad)
    row = cell_sign_to_row(s)
    idx = interfacequads.celltoquad[cellid]
    interfacequads.quads[row, idx] = quad
end

function update_face_quadrature!(facequads::FaceQuadratures, s, faceid, cellid, quad)
    row = cell_sign_to_row(s)
    idx = facequads.facetoquad[row, faceid, cellid]
    @assert idx > 4 "Attempting to modify the uniform cell face quadrature rules"
    facequads.quads[idx] = quad
end

function map_and_update_quadrature!(cellquads, s, cellid, mergemapper, mapid)
    quad = cellquads[s, cellid]
    newquad = map_quadrature(quad, mergemapper, mapid)
    update_quadrature!(cellquads, s, cellid, newquad)
end

function map_and_update_face_quadrature!(facequads, s, cellid, mergemapper, mapid)
    nfaces = number_of_faces_per_cell(facequads)
    for faceid = 1:nfaces
        quad = facequads[s, faceid, cellid]
        newquad = map_quadrature(quad, mergemapper, mapid)
        update_face_quadrature!(facequads, s, faceid, cellid, newquad)
    end
end

struct MergeCutMesh
    mergedwithcell::Any
    cutmesh::Any
    function MergeCutMesh(cutmesh::CutMesh)
        ncells = number_of_cells(cutmesh)
        mergedwithcell = vcat((1:ncells)', (1:ncells)')
        new(mergedwithcell, cutmesh)
    end
end

function number_of_cells(mergecutmesh::MergeCutMesh)
    return number_of_cells(mergecutmesh.cutmesh)
end

function cell_sign(mergecutmesh::MergeCutMesh)
    return cell_sign(mergecutmesh.cutmesh)
end

function cell_sign(mergecutmesh::MergeCutMesh,cellid)
    return cell_sign(mergecutmesh.cutmesh,cellid)
end

function dimension(mergecutmesh::MergeCutMesh)
    return dimension(mergecutmesh.cutmesh)
end

function cell_map(mergecutmesh::MergeCutMesh, s, cellid)
    row = cell_sign_to_row(s)
    mergecellid = mergecutmesh.mergedwithcell[row, cellid]
    return cell_map(mergecutmesh.cutmesh, mergecellid)
end

function cell_map(mergecutmesh::MergeCutMesh, cellid)
    return cell_map(mergecutmesh.cutmesh, cellid)
end

function is_interior_cell(mergecutmesh::MergeCutMesh)
    return is_interior_cell(mergecutmesh.cutmesh)
end

function number_of_nodes(mergecutmesh::MergeCutMesh)
    cellsign = cell_sign(mergecutmesh)
    numnodes = 0
    for (cellid, s) in enumerate(cellsign)
        if s == +1 || s == 0
            nodeids = nodal_connectivity(mergecutmesh, +1, cellid)
            numnodes = max(numnodes, maximum(nodeids))
        end
        if s == -1 || s == 0
            nodeids = nodal_connectivity(mergecutmesh, -1, cellid)
            numnodes = max(numnodes, maximum(nodeids))
        end
    end
    return numnodes
end

function number_of_degrees_of_freedom(mergecutmesh::MergeCutMesh)
    dim = dimension(mergecutmesh)
    numnodes = number_of_nodes(mergecutmesh)
    return dim * numnodes
end

function merge_cells!(mergecutmesh::MergeCutMesh, s, mergeto, mergefrom)
    row = cell_sign_to_row(s)
    ncells = number_of_cells(mergecutmesh)
    @assert 1 <= mergeto <= ncells
    @assert 1 <= mergefrom <= ncells

    mergecutmesh.mergedwithcell[row, mergefrom] = mergeto
end

function Base.show(io::IO, mergecutmesh::MergeCutMesh)
    ncells = number_of_cells(mergecutmesh.cutmesh)
    str = "MergeCutMesh\n\tNum. Cells: $ncells"
    print(io, str)
end

function nodal_connectivity(mergecutmesh::MergeCutMesh, s, cellid)
    row = cell_sign_to_row(s)
    mergecellid = mergecutmesh.mergedwithcell[row, cellid]
    return nodal_connectivity(mergecutmesh.cutmesh, s, mergecellid)
end

function cell_connectivity(mergecutmesh::MergeCutMesh)
    return cell_connectivity(mergecutmesh.cutmesh)
end

function active_cells(mergecutmesh::MergeCutMesh)
    return active_cells(mergecutmesh.cutmesh)
end

function quadrature_areas(cellquads)
    ncells = number_of_cells(cellquads)
    areas = zeros(2, ncells)
    for cellid = 1:ncells
        for s in [+1, -1]
            if has_quadrature(cellquads, s, cellid)
                row = cell_sign_to_row(s)
                areas[row, cellid] = sum(cellquads[s, cellid].weights)
            end
        end
    end
    return areas
end

function is_tiny_cell(cellquads, quadareas; ratio = 0.2)
    ncells = number_of_cells(cellquads)
    istiny = zeros(Bool, 2, ncells)

    tinyarea = ratio * sum(uniform_cell_quadrature(cellquads).weights)

    for cellid = 1:ncells
        for s in [+1, -1]
            if has_quadrature(cellquads, s, cellid)
                row = cell_sign_to_row(s)
                if quadareas[row, cellid] <= tinyarea
                    istiny[row, cellid] = true
                end
            end
        end
    end

    return istiny
end

function merge_cells_in_mesh!(mergecutmesh, cellquads, interfacequads, mergemapper)

    cellconnectivity = cell_connectivity(mergecutmesh)
    activecells = active_cells(mergecutmesh)
    ncells = number_of_cells(mergecutmesh)

    quadareas = quadrature_areas(cellquads)
    istinycell = is_tiny_cell(cellquads, quadareas)

    for cellid = 1:ncells
        for s in [+1, -1]
            row = cell_sign_to_row(s)
            if istinycell[row, cellid]
                nbrcellids = cellconnectivity[:, cellid]
                nbractive =
                    [nbrid == 0 ? false : activecells[row, nbrid] for nbrid in nbrcellids]
                nbrnontiny =
                    [nbrid == 0 ? false : !istinycell[row, nbrid] for nbrid in nbrcellids]
                faceid = findfirst(nbractive .& nbrnontiny)
                !isnothing(faceid) ||
                    error("Could not find appropriate merge direction for cell $cellid, cellsign $s")
                oppositeface = opposite_face(faceid)
                mergecellid = nbrcellids[faceid]

                merge_cells!(mergecutmesh, s, mergecellid, cellid)
                map_and_update_quadrature!(cellquads, s, cellid, mergemapper, oppositeface)
                map_and_update_quadrature!(
                    interfacequads,
                    s,
                    cellid,
                    mergemapper,
                    oppositeface,
                )
            end
        end
    end
end

function merge_cells_in_mesh!(
    mergecutmesh,
    cellquads,
    interfacequads,
    facequads,
    mergemapper,
)

    cellconnectivity = cell_connectivity(mergecutmesh)
    activecells = active_cells(mergecutmesh)
    ncells = number_of_cells(mergecutmesh)

    quadareas = quadrature_areas(cellquads)
    istinycell = is_tiny_cell(cellquads, quadareas)

    for cellid = 1:ncells
        for s in [+1, -1]
            row = cell_sign_to_row(s)
            if istinycell[row, cellid]
                nbrcellids = cellconnectivity[:, cellid]
                nbractive =
                    [nbrid == 0 ? false : activecells[row, nbrid] for nbrid in nbrcellids]
                nbrnontiny =
                    [nbrid == 0 ? false : !istinycell[row, nbrid] for nbrid in nbrcellids]
                faceid = findfirst(nbractive .& nbrnontiny)
                !isnothing(faceid) ||
                    error("Could not find appropriate merge direction for cell $cellid, cellsign $s")
                oppositeface = opposite_face(faceid)
                mergecellid = nbrcellids[faceid]

                merge_cells!(mergecutmesh, s, mergecellid, cellid)
                map_and_update_quadrature!(cellquads, s, cellid, mergemapper, oppositeface)
                map_and_update_quadrature!(
                    interfacequads,
                    s,
                    cellid,
                    mergemapper,
                    oppositeface,
                )
                map_and_update_face_quadrature!(facequads,s,cellid,mergemapper,oppositeface)
            end
        end
    end
end

struct MergedMesh
    mergecutmesh::MergeCutMesh
    nodelabeltonodeid::Any
    numnodes::Any
    function MergedMesh(mergecutmesh::MergeCutMesh)
        activenodelabels = active_node_labels(mergecutmesh)
        maxlabel = maximum(activenodelabels)
        nodelabeltonodeid = zeros(Int, maxlabel)
        for (idx, label) in enumerate(activenodelabels)
            nodelabeltonodeid[label] = idx
        end
        numnodes = length(activenodelabels)
        new(mergecutmesh, nodelabeltonodeid, numnodes)
    end
end

function number_of_cells(mergedmesh::MergedMesh)
    return number_of_cells(mergedmesh.mergecutmesh)
end

function Base.show(io::IO, mergedmesh::MergedMesh)
    ncells = number_of_cells(mergedmesh)
    numnodes = mergedmesh.numnodes
    str = "MergedMesh\n\tNum. Cells: $ncells\n\tNum. Nodes: $numnodes"
    print(io, str)
end

function dimension(mergedmesh::MergedMesh)
    return dimension(mergedmesh.mergecutmesh)
end

function cell_sign(mergedmesh::MergedMesh)
    return cell_sign(mergedmesh.mergecutmesh)
end

function cell_sign(mergedmesh::MergedMesh,cellid)
    return cell_sign(mergedmesh.mergecutmesh,cellid)
end

function is_interior_cell(mergedmesh::MergedMesh)
    return is_interior_cell(mergedmesh.mergecutmesh)
end

function nodal_connectivity(mergedmesh::MergedMesh, s, cellid)
    labels = nodal_connectivity(mergedmesh.mergecutmesh, s, cellid)
    nodeids = mergedmesh.nodelabeltonodeid[labels]
    return nodeids
end

function cell_map(mergedmesh::MergedMesh, s, idx)
    return cell_map(mergedmesh.mergecutmesh, s, idx)
end

function cell_map(mergedmesh::MergedMesh,cellid)
    return cell_map(mergedmesh.mergecutmesh,cellid)
end

function number_of_nodes(mergedmesh::MergedMesh)
    return mergedmesh.numnodes
end

function cell_connectivity(mergedmesh::MergedMesh)
    return cell_connectivity(mergedmesh.mergecutmesh)
end

function number_of_degrees_of_freedom(mergedmesh::MergedMesh)
    dim = dimension(mergedmesh)
    numnodes = number_of_nodes(mergedmesh)
    return dim * numnodes
end

function active_node_labels(mergecutmesh::MergeCutMesh)
    numcells = number_of_cells(mergecutmesh)
    cellsign = cell_sign(mergecutmesh)

    activenodeids = Int[]

    cellids = findall(x -> x == 0 || x == +1, cellsign)
    for cellid in cellids
        append!(activenodeids, nodal_connectivity(mergecutmesh, +1, cellid))
    end

    cellids = findall(x -> x == 0 || x == -1, cellsign)
    for cellid in cellids
        append!(activenodeids, nodal_connectivity(mergecutmesh, -1, cellid))
    end

    sort!(activenodeids)
    unique!(activenodeids)
    return activenodeids
end


text:
const polyCoeffs17 = [2305843009213693950, 2301339409586323472, 2305812222886150007, 2304950122020066599, 432341041511554329, 218546466595982445, 2288774021294666627, 1327872805097898761, 1780609199527208287, 1649816982904334450, 870691130730945933, 845235595464045280, 1023515427770432269, 2127882104260287640, 2304992536952056903, 2305796279969513047, 2305843009213693918, 1]

#ei coefficients for m=2^5...2^30
const skipMat17 = [
    [
        1972101396566503791,
        2051641146438145039,
        533173487516893382,
        545730390508709012,
        1620905962071234677,
        1870596525612362239,
        2112925984824714998,
        868498104066040073,
        1148789919911044037,
        804191351341168970,
        730557211712521605,
        1517348237025551691,
        850529375445951650,
        717602406832799554,
        2195246085431969997,
        366890617173684729,
        1334285358249915285,
    ],
    [
        327007261501768800,
        1668507619896500323,
        570663528499880382,
        849832518716040141,
        167129205132215511,
        2029511566136161224,
        1622479814610858826,
        136068578967795545,
        465570142067313353,
        443857812536418275,
        1789415231111220632,
        722738813040608858,
        1536033290496530295,
        1803331387035840455,
        1337554956082899184,
        1262248399564869336,
        1964812692109901433,
    ],
    [
        1383941110888204807,
        1660111263786065031,
        635055952618544659,
        1275158498626512324,
        897543831095138544,
        1662589174551273379,
        1474986199793476754,
        683298560121080184,
        1364215338737587469,
        2093515129575752146,
        1526800205169804711,
        2260714998735099815,
        2305114977222035959,
        2007484595445621488,
        1660849571985359945,
        413454328436564294,
        1737301931246044859,
    ],
    [
        34041530522269951,
        989904588027666109,
        1424411770440999775,
        2233548010652125537,
        478899446812575694,
        1782594547279649086,
        1366273273270145235,
        771683963354479181,
        1042197874804458977,
        47908093511613389,
        41949539543464032,
        2092607021509534694,
        1553302290867055291,
        736326122801287713,
        1087073336441938432,
        1684732520189930480,
        1523208273561944990,
    ],
    [
        66724810607261494,
        2243218526996946637,
        1995426020901038827,
        1273288642662700447,
        1060296628596275187,
        1033117315917689018,
        1301796675675435740,
        1169744488162372532,
        1921449921433723255,
        1937495610468571249,
        993465857439554119,
        16019142551863611,
        2219970175835711455,
        1006443666620132302,
        924296262882437568,
        2208413115601089136,
        264904656702381098,
    ],
    [
        86455199185446885,
        747612738103107959,
        1407924622912058943,
        198804599081851031,
        25930318115067972,
        412979531876469963,
        715204133791195588,
        1086783985939604076,
        804486103308535628,
        137063081601334719,
        822767962759173950,
        1264697699255051205,
        1772814692491721476,
        2011485527677807227,
        399961272695350873,
        373953638541612193,
        11760493978029076,
    ],
    [
        1708255563439812290,
        774500572940894871,
        99516351798717781,
        231082297833576357,
        2060166097000388081,
        1966934547960002972,
        1891076714188958609,
        37390213127216786,
        1584006554430522049,
        890578844293885408,
        1954994145879248940,
        2062392734010795911,
        1020363785670620886,
        2100645351290474524,
        1589146708828507128,
        1348093719303308944,
        180639173108775193,
    ],
    [
        2202229385283005142,
        1571862202556900423,
        368225414353564040,
        119278548399139866,
        952453769547029233,
        1720919561984421856,
        1430751071426459339,
        1307940815395225146,
        1377228865567960952,
        1990790714390702490,
        677872157032915967,
        390217894376331041,
        738924010463110026,
        224005659670894829,
        2067169282375116419,
        1185837644893947828,
        628166474873784103,
    ],
    [
        909535552350355334,
        1984196621966487204,
        1203870959107599346,
        1777238049233603922,
        717242522487790393,
        345373014150664306,
        285368686010812695,
        1609391685732340886,
        739922468389050788,
        234977750553506588,
        1990140987404995049,
        1604923589910519614,
        2278523533111581954,
        180283447256677998,
        1068630207724457873,
        304098965825306400,
        2062051512419880902,
    ],
    [
        1422118065708901004,
        836642887813982710,
        1888061742065740657,
        274692856042566127,
        1356911043266965940,
        976761670652707464,
        531998614533163061,
        1912546600455521925,
        17045467741578191,
        1654907078127935108,
        2030183713003610570,
        352607571297516120,
        355680148820201396,
        739542852109405450,
        709616794056412549,
        1336256981892236022,
        2069426500410221871,
    ],
    [
        2286984874264166784,
        2166094372062800745,
        879293621066552442,
        2278089229880928468,
        1994187841709719990,
        2093302072134305223,
        361239617152043736,
        2130415391463653737,
        1938772631940986057,
        805733948963353760,
        1424018902747763999,
        773476430100927128,
        2124442120407728160,
        121106965413699034,
        475872283540920661,
        1044131603014983440,
        1923430800318934442,
    ],
    [
        451619876923000687,
        71097883912081557,
        1628549684595849797,
        41885163308977346,
        362429429797508426,
        971288200022972701,
        1418629218562813163,
        752202962754801305,
        2086557484323717914,
        626769727155832314,
        1899885130221156377,
        298197818644716470,
        245633358844767262,
        976528303012507716,
        1559467768505760317,
        44483913443667677,
        1725590647900820418,
    ],
    [
        661193250562452382,
        1043865391393629623,
        607105552343853366,
        1082113160102069187,
        1808845276619195583,
        2140932414351717404,
        852505500160497633,
        1975416625595591040,
        1230888852709529781,
        523574463061188471,
        176632315048861766,
        735704925835568324,
        821608995808733894,
        2236293070564040013,
        1772032514358201802,
        1816015347932134306,
        1317895217324711671,
    ],
    [
        620079138192442768,
        1413852673829294104,
        139241503827318837,
        890594028110926166,
        977736786367659953,
        742757423130029561,
        23963199723774009,
        951576424938956514,
        237070967500044720,
        709309519512381854,
        1217415377215609678,
        2010236133988708113,
        1922119966574907981,
        2116003072079736225,
        1058675710956122435,
        2291324972928830992,
        1412559516146851951,
    ],
    [
        1404725383601173936,
        1109185334627610252,
        2177895433232027303,
        1454696689462755000,
        559089231530930806,
        2113042216679845782,
        2212705677772153535,
        58693736675718155,
        745821794534861716,
        2071125041882463283,
        1291815773765129120,
        1446497568380515246,
        328849553205520604,
        1572459871580105595,
        770521934324544238,
        2162594601603215507,
        1790149936455713032,
    ],
    [
        843356255388265439,
        2303128525550447885,
        1403805614380726589,
        1646775209363665369,
        1656489684429545237,
        136962756483849705,
        578486145972684395,
        1186788852536058806,
        744606065061259291,
        129235498787499648,
        987271125576440333,
        2226012411473073226,
        1064839699140268466,
        1194475380559801816,
        1591779369875059359,
        1167968540715842034,
        2137709594798027207,
    ],
    [
        1972101396566503791,
        2051641146438145039,
        533173487516893382,
        545730390508709012,
        1620905962071234677,
        1870596525612362239,
        2112925984824714998,
        868498104066040073,
        1148789919911044037,
        804191351341168970,
        730557211712521605,
        1517348237025551691,
        850529375445951650,
        717602406832799554,
        2195246085431969997,
        366890617173684729,
        1334285358249915285,
    ],
    [
        327007261501768800,
        1668507619896500323,
        570663528499880382,
        849832518716040141,
        167129205132215511,
        2029511566136161224,
        1622479814610858826,
        136068578967795545,
        465570142067313353,
        443857812536418275,
        1789415231111220632,
        722738813040608858,
        1536033290496530295,
        1803331387035840455,
        1337554956082899184,
        1262248399564869336,
        1964812692109901433,
    ],
    [
        1383941110888204807,
        1660111263786065031,
        635055952618544659,
        1275158498626512324,
        897543831095138544,
        1662589174551273379,
        1474986199793476754,
        683298560121080184,
        1364215338737587469,
        2093515129575752146,
        1526800205169804711,
        2260714998735099815,
        2305114977222035959,
        2007484595445621488,
        1660849571985359945,
        413454328436564294,
        1737301931246044859,
    ],
    [
        34041530522269951,
        989904588027666109,
        1424411770440999775,
        2233548010652125537,
        478899446812575694,
        1782594547279649086,
        1366273273270145235,
        771683963354479181,
        1042197874804458977,
        47908093511613389,
        41949539543464032,
        2092607021509534694,
        1553302290867055291,
        736326122801287713,
        1087073336441938432,
        1684732520189930480,
        1523208273561944990,
    ],
    [
        66724810607261494,
        2243218526996946637,
        1995426020901038827,
        1273288642662700447,
        1060296628596275187,
        1033117315917689018,
        1301796675675435740,
        1169744488162372532,
        1921449921433723255,
        1937495610468571249,
        993465857439554119,
        16019142551863611,
        2219970175835711455,
        1006443666620132302,
        924296262882437568,
        2208413115601089136,
        264904656702381098,
    ],
    [
        86455199185446885,
        747612738103107959,
        1407924622912058943,
        198804599081851031,
        25930318115067972,
        412979531876469963,
        715204133791195588,
        1086783985939604076,
        804486103308535628,
        137063081601334719,
        822767962759173950,
        1264697699255051205,
        1772814692491721476,
        2011485527677807227,
        399961272695350873,
        373953638541612193,
        11760493978029076,
    ],
    [
        1708255563439812290,
        774500572940894871,
        99516351798717781,
        231082297833576357,
        2060166097000388081,
        1966934547960002972,
        1891076714188958609,
        37390213127216786,
        1584006554430522049,
        890578844293885408,
        1954994145879248940,
        2062392734010795911,
        1020363785670620886,
        2100645351290474524,
        1589146708828507128,
        1348093719303308944,
        180639173108775193,
    ],
    [
        2202229385283005142,
        1571862202556900423,
        368225414353564040,
        119278548399139866,
        952453769547029233,
        1720919561984421856,
        1430751071426459339,
        1307940815395225146,
        1377228865567960952,
        1990790714390702490,
        677872157032915967,
        390217894376331041,
        738924010463110026,
        224005659670894829,
        2067169282375116419,
        1185837644893947828,
        628166474873784103,
    ],
    [
        909535552350355334,
        1984196621966487204,
        1203870959107599346,
        1777238049233603922,
        717242522487790393,
        345373014150664306,
        285368686010812695,
        1609391685732340886,
        739922468389050788,
        234977750553506588,
        1990140987404995049,
        1604923589910519614,
        2278523533111581954,
        180283447256677998,
        1068630207724457873,
        304098965825306400,
        2062051512419880902,
    ],
    [
        1422118065708901004,
        836642887813982710,
        1888061742065740657,
        274692856042566127,
        1356911043266965940,
        976761670652707464,
        531998614533163061,
        1912546600455521925,
        17045467741578191,
        1654907078127935108,
        2030183713003610570,
        352607571297516120,
        355680148820201396,
        739542852109405450,
        709616794056412549,
        1336256981892236022,
        2069426500410221871,
    ],
    [
        2286984874264166784,
        2166094372062800745,
        879293621066552442,
        2278089229880928468,
        1994187841709719990,
        2093302072134305223,
        361239617152043736,
        2130415391463653737,
        1938772631940986057,
        805733948963353760,
        1424018902747763999,
        773476430100927128,
        2124442120407728160,
        121106965413699034,
        475872283540920661,
        1044131603014983440,
        1923430800318934442,
    ],
    [
        451619876923000687,
        71097883912081557,
        1628549684595849797,
        41885163308977346,
        362429429797508426,
        971288200022972701,
        1418629218562813163,
        752202962754801305,
        2086557484323717914,
        626769727155832314,
        1899885130221156377,
        298197818644716470,
        245633358844767262,
        976528303012507716,
        1559467768505760317,
        44483913443667677,
        1725590647900820418,
    ],
    [
        661193250562452382,
        1043865391393629623,
        607105552343853366,
        1082113160102069187,
        1808845276619195583,
        2140932414351717404,
        852505500160497633,
        1975416625595591040,
        1230888852709529781,
        523574463061188471,
        176632315048861766,
        735704925835568324,
        821608995808733894,
        2236293070564040013,
        1772032514358201802,
        1816015347932134306,
        1317895217324711671,
    ],
    [
        620079138192442768,
        1413852673829294104,
        139241503827318837,
        890594028110926166,
        977736786367659953,
        742757423130029561,
        23963199723774009,
        951576424938956514,
        237070967500044720,
        709309519512381854,
        1217415377215609678,
        2010236133988708113,
        1922119966574907981,
        2116003072079736225,
        1058675710956122435,
        2291324972928830992,
        1412559516146851951,
    ],
    [
        1404725383601173936,
        1109185334627610252,
        2177895433232027303,
        1454696689462755000,
        559089231530930806,
        2113042216679845782,
        2212705677772153535,
        58693736675718155,
        745821794534861716,
        2071125041882463283,
        1291815773765129120,
        1446497568380515246,
        328849553205520604,
        1572459871580105595,
        770521934324544238,
        2162594601603215507,
        1790149936455713032,
    ],
    [
        843356255388265439,
        2303128525550447885,
        1403805614380726589,
        1646775209363665369,
        1656489684429545237,
        136962756483849705,
        578486145972684395,
        1186788852536058806,
        744606065061259291,
        129235498787499648,
        987271125576440333,
        2226012411473073226,
        1064839699140268466,
        1194475380559801816,
        1591779369875059359,
        1167968540715842034,
        2137709594798027207,
    ],
    [
        1603310798512369292,
        1601100715499888665,
        1072379507175706577,
        1213079833732518096,
        2001446403012491729,
        419527915496992901,
        1536979475280290440,
        301233389316038205,
        1272639558366606573,
        1557584358952356730,
        220291452971595398,
        2054031811190278086,
        2102439523114162782,
        2043691020304126326,
        931868944669449976,
        1695482396558042263,
        569471942240552669,
    ],
    [
        694964758088788671,
        1038219072747505180,
        1918123200517156133,
        2205661282197833706,
        273525348357470113,
        2121494391145814980,
        284345776977803033,
        246418780931512795,
        2283919448816535183,
        1825147280778299237,
        1690836788966579901,
        668223754085848325,
        1157160865829805369,
        948478357813085584,
        254359412607383270,
        1135224685685690391,
        2000460937859975771,
    ],
    [
        707900843551449518,
        824721369516455951,
        1253414446894639545,
        2054381172664740386,
        2262749485189772637,
        209547095357435226,
        1486285367660805487,
        1066212727062337001,
        138467865375837985,
        502055498606674686,
        1631322852248214521,
        2099497052471449902,
        653339673532849910,
        1209350158362883423,
        801199054076992563,
        114429634116085316,
        1756481675028861979,
    ],
    [
        1233556889083976385,
        2271316709310852100,
        444255811069641330,
        1920402384108973078,
        1871076587612059977,
        848766558121280411,
        183916672241888601,
        861892142014135059,
        1662987476773654122,
        886432897949545677,
        952076682704709636,
        2191551637233270217,
        1294187932820418508,
        691043740128684730,
        1580114045565902135,
        871127380587503701,
        1982811410939956429,
    ],
    [
        1026536596545411590,
        672006064876756356,
        187865928299266937,
        1333180393912329299,
        449284825643474283,
        1880633394205436049,
        2130751968571655233,
        1666493909078411462,
        1226379145018641537,
        925073470477791740,
        1893981127322123262,
        1681869522295985310,
        563899279485833974,
        1238859114847806213,
        1910620991062119273,
        573004330562037184,
        2305020818618415056,
    ],
    [
        1649536426891372496,
        1773496319691523181,
        1901485078913394927,
        723526059786718651,
        1070993657770771508,
        561598982654578851,
        514075350227526556,
        1093001656896099183,
        1131031199773443838,
        1567007165576640890,
        1942807866844241746,
        1377184921804775466,
        525984842795669092,
        1310907835907728052,
        1744624221468767441,
        1158478818307902549,
        178013669142197363,
    ],
    [
        348651919527231585,
        1956322108668706848,
        2200609920440762066,
        1722438664965020838,
        2287184028526454944,
        1503209766683435682,
        644324372269330068,
        839244854502053928,
        822962662761844987,
        967028104753664904,
        2219260946362335472,
        643036698992345239,
        767336458352242411,
        1856163925694015939,
        310679063625611629,
        1008184523753394922,
        1557936276546470810,
    ],
    [
        1358639256207203914,
        380878389447600344,
        266118891978454911,
        1913066921087211778,
        1415556213238755288,
        1041859734009513367,
        913429674643150097,
        1281916063862733602,
        649226289994219281,
        2012458594853276202,
        806265659335377279,
        488631900297180252,
        1174216778495607541,
        66231138464254448,
        102026227548073028,
        6202600785403592,
        1876707885948960382,
    ],
    [
        302708009212573930,
        640510323803542134,
        390948779667971340,
        300092081045236432,
        1387655392749610816,
        973283902592015763,
        1559965596923378008,
        1181369898028308254,
        914040186129074102,
        900388984624711127,
        409569314904400659,
        423416875652327473,
        1646832545602879919,
        1122149242942110197,
        1318635598796713194,
        27210210178262476,
        851295614470433169,
    ],
    [
        1991387975529265218,
        1843351940963099488,
        2271846316599365716,
        2132592891396609803,
        2302343546822788973,
        2129811222706328476,
        1254066806718120786,
        2268141813723625983,
        1450651840911637570,
        726351353821067277,
        2126412864417917107,
        952525217820829278,
        1116731824858793647,
        2018141859931748691,
        1074656985216627374,
        536391227666488033,
        346611616870907958,
    ],
]


text:
@testset "tsptwstaterepresentation.jl" begin
    @testset "Constructor" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        n_city = 3
        grid_size = 10
        max_tw_gap = 3
        max_tw = 8

        generator = SeaPearl.TsptwGenerator(n_city, grid_size, max_tw_gap, max_tw)

        SeaPearl.fill_with_generator!(model, generator; seed=42)

        dist, timeWindows, pos = model.adhocInfo

        # This condition is there because of the way random are generated can change from one version to another
        if VERSION >= v"1.6.0"
            @test dist ==  [0 8 10;
                            8 0 4;
                            10 4 0]
            @test timeWindows == [0 10; 8 16; 12 15]
            @test pos == [5.331830160438614 1.7293302893695128; 4.540291355871425 9.589258763297348; 0.17686826714964354 9.735659798036858]

            sr = SeaPearl.TsptwStateRepresentation(model)

            @test sr.dist == [  0.0 0.8 1.0
                                0.8 0.0 0.4
                                1.0 0.4 0.0]
            @test trunc.(sr.timeWindows; digits=2) == [0.0  0.62
                                                        0.5 1.0
                                                        0.75 0.93]
            @test trunc.(sr.pos; digits=2) == [ 0.53 0.17
                                                0.45 0.95
                                                0.01 0.97]
            @test trunc.(sr.nodeFeatures; digits=2) == transpose(Float32[0.53 0.17 0.0 0.62 0.0 0.0; 0.45 0.95 0.5 1.0 0.0 0.0; 0.01 0.97 0.75 0.93 0.0 0.0])
        end
    end
end


text:
using ExtensibleEffects
using DataTypesBasic
using Monadic

f(a, b) = @syntax_eff begin
    a = a
    b = b
    @pure a, b
end

f([1,2,3], Identity(99))

@code_warntype f([1,2,3], Identity(99))
@code_native f([1,2,3], Identity(99))

codetest(x) = @syntax_eff begin
  v = [1,3,4]
  @pure a = x+v
  o = isodd(a) ? Option(100) : Option()
  @pure b = x + a + o
  @pure v, a, o, b
end

@code_native codetest(1)
comparef(x) = [isodd(v+x) ? (v, v+x, 100, v+x+100) : nothing for v in [1,3,4]]
@code_native comparef(1)  # at least same code

@code_warntype codetest(1)  # inference does not work yet...

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
using Documenter
using GarishPrint
using DocThemeIndigo

indigo = DocThemeIndigo.install(GarishPrint)

makedocs(;
    modules = [GarishPrint],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical="https://Roger-luo.github.io/GarishPrint.jl",
        assets=String[indigo],
    ),
    pages = [
        "Home" => "index.md",
        "References" => "ref.md",
    ],
    repo = "https://github.com/Roger-luo/GarishPrint.jl",
    sitename = "GarishPrint.jl",
)

deploydocs(; repo = "github.com/Roger-luo/GarishPrint.jl")


text:
# The default constructors for Periods work well in almost all cases
# P(x) = new((convert(Int64,x))
# The following definitions are for Period-specific safety
for period in (:Year, :Semester, :Quarter, :Month, :Week, :Day)
    period_str = string(period)
    accessor_str = lowercase(period_str)
    # Convenience method for show()
    @eval _units(x::$period) = " " * $accessor_str * (abs(value(x)) == 1 ? "" : "s")
    # periodisless
    @eval periodisless(x::$period, y::$period) = value(x) < value(y)
    # Period accessors
    typs = [:YearDate, :SemesterDate, :QuarterDate, :MonthDate, :WeekDate, :DayDate]
    reference =
        period === :Week ?
        " For details see [`$accessor_str(::Union{Date, DateTime})`](@ref)." : ""
    for typ_str in typs
        @eval begin
            @doc """
                $($period_str)(dt::$($typ_str)) -> $($period_str)

            The $($accessor_str) part of a $($typ_str) as a `$($period_str)`.$($reference)
            """ $period(dt::$(Symbol(typ_str))) = $period($(Symbol(accessor_str))(dt))
        end
    end
    @eval begin
        @doc """
            $($period_str)(v)

        Construct a `$($period_str)` object with the given `v` value. Input must be
        losslessly convertible to an [`Int64`](@ref).
        """ $period(v)
    end
end

# AbstractString parsing (mainly for IO code)
Semester(x::AbstractString) = Semester(Base.parse(Int64, x))
Undated(x::AbstractString) = Undated(Base.parse(Int64, x))
# The period type is printed when output, thus it already implies its own typeinfo
Base.typeinfo_implicit(::Type{Semester}) = true
Base.typeinfo_implicit(::Type{Undated}) = true

_units(u::Undated) = ""

#Print/show/traits
Base.print(io::IO, x::Semester) = print(io, value(x), _units(x))
Base.print(io::IO, x::Undated) = print(io, value(x), _units(x))

# See https://en.wikipedia.org/wiki/ISO_week_date
# There are 20871 weeks in 400 years
weeks(c::Day) = div(value(c), 7)
weeks(c::Week) = value(c)
weeks(c::Month) = 4.348125 * value(c)
weeks(c::Quarter) = 13.044375 * value(c)
weeks(c::Semester) = 26.08875 * value(c)
weeks(c::Year) = 52.1775 * value(c)
months(c::Day) = div(value(c), 30.436875)
months(c::Week) = div(value(c), 4.34812)
months(c::Month) = value(c)
months(c::Quarter) = 3 * value(c)
months(c::Semester) = 6 * value(c)
months(c::Year) = 12 * value(c)
quarters(c::Day) = div(value(c), 91.310625)
quarters(c::Week) = div(value(c), 13.044385)
quarters(c::Month) = div(value(c), 3)
quarters(c::Quarter) = value(c)
quarters(c::Semester) = 2 * value(c)
quarters(c::Year) = 4 * value(c)
semesters(c::Day) = div(value(c), 182.62125)
semesters(c::Week) = div(value(c), 26.08877)
semesters(c::Month) = div(value(c), 6)
semesters(c::Quarter) = div(value(c), 2)
semesters(c::Semester) = value(c)
semesters(c::Year) = 2 * value(c)
years(c::Day) = div(value(c), 365.2425)
years(c::Week) = div(value(c), 52.1775)
years(c::Month) = div(value(c), 12)
years(c::Quarter) = div(value(c), 4)
years(c::Semester) = div(value(c), 2)
years(c::Year) = value(c)

# Default values (as used by SimplePeriods)
"""
    default(p::Period) -> Period

Returns a sensible "default" value for the input Period by returning `T(1)` for Year,
Semeter, Quarter, Month, Day and Undated.
"""
function default end

default(p::Union{T,Type{T}}) where {T<:Dates.DatePeriod} = T(1)

# like div but throw an error if remainder is nonzero
function divexact(x, y)
    q, r = divrem(x, y)
    r == 0 || throw(InexactError(:divexact, Int, x / y))
    return q
end

const zero_or_fixedperiod_seed = UInt === UInt64 ? 0x5b7fc751bba97516 : 0xeae0fdcb
const nonzero_otherperiod_seed = UInt === UInt64 ? 0xe1837356ff2d2ac9 : 0x170d1b00
otherperiod_seed(x::OtherPeriod) =
    iszero(value(x)) ? zero_or_fixedperiod_seed : nonzero_otherperiod_seed

Base.hash(x::Semester, h::UInt) = hash(6 * value(x), h + otherperiod_seed(x))


text:
@testset "COOM (Co-occurence Matrix)" begin
    doc_raw = "This is a document. It has two sentences."
    doc = prepare(doc_raw, strip_punctuation|strip_whitespace|strip_case)
    sd = StringDocument{String}(doc)
    td = TokenDocument{String}(doc)
    nd = NGramDocument{String}(doc)
    crps = Corpus([sd, td])
    T = Float16
    # Results for window = 5, all terms in document used
    expected_result = [ # for window == 5
        0.0 2.0 1.0 2/3 0.5 0.4 0.0 0.0
        2.0 0.0 2.0 1.0 2/3 0.5 0.4 0.0
        1.0 2.0 0.0 2.0 1.0 2/3 0.5 0.4
        2/3 1.0 2.0 0.0 2.0 1.0 2/3 0.5
        0.5 2/3 1.0 2.0 0.0 2.0 1.0 2/3
        0.4 0.5 2/3 1.0 2.0 0.0 2.0 1.0
        0.0 0.4 0.5 2/3 1.0 2.0 0.0 2.0
        0.0 0.0 0.4 0.5 2/3 1.0 2.0 0.0]
    # Verify untyped constructor
    terms = tokens(td)
    for d in [doc, sd, td, crps]
        C = CooMatrix(d, terms)
        @test C isa CooMatrix{StringAnalysis.DEFAULT_FLOAT_TYPE}
        if !(d isa Corpus)
            @test coom(C) == expected_result
        else
            @test coom(C) == length(crps) * expected_result
        end
    end
    @test_throws ErrorException CooMatrix(nd)

    # Verify typed constructor
    terms = tokens(td)
    for d in [doc, sd, td, crps]
        C = CooMatrix{T}(d, terms)
        @test C isa CooMatrix{T}
        if !(d isa Corpus)
            @test coom(C) == T.(expected_result)
        else
            @test coom(C) == length(crps) * T.(expected_result)
        end
    end
    @test_throws ErrorException CooMatrix{T}(nd)

    # Results for window = 1, custom terms
    terms = ["this", "document", "it"]
    expected_result = [0.0 0.0 0.0; # document
                       0.0 0.0 2.0; # it
                       0.0 2.0 0.0] # this
    # Verify untyped constructor
    for d in [doc, sd, td, crps]
        C = CooMatrix{T}(d, terms, window=1)
        @test C isa CooMatrix{T}
        if !(d isa Corpus)
            @test coom(C) == T.(expected_result)
        else
            @test coom(C) == length(crps) * T.(expected_result)
        end
    end
    @test_throws ErrorException CooMatrix{T}(nd)
end


text:
using DiffEqBenchmarks
if length(ARGS) == 1
    println("Benchmarking the $(ARGS[1]) folder")
    DiffEqBenchmarks.weave_folder(ARGS[1])
elseif length(ARGS) == 2
    println("Benchmarking $(ARGS[1])/$(ARGS[2])")
    DiffEqBenchmarks.weave_file(ARGS[1],ARGS[2])
else
    error("Only 1 or 2 arguments allowed!")
end


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

mutable struct IoK8sApiCoreV1Probe <: SwaggerModel
    exec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ExecAction } # spec name: exec
    failureThreshold::Any # spec type: Union{ Nothing, Int32 } # spec name: failureThreshold
    httpGet::Any # spec type: Union{ Nothing, IoK8sApiCoreV1HTTPGetAction } # spec name: httpGet
    initialDelaySeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: initialDelaySeconds
    periodSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: periodSeconds
    successThreshold::Any # spec type: Union{ Nothing, Int32 } # spec name: successThreshold
    tcpSocket::Any # spec type: Union{ Nothing, IoK8sApiCoreV1TCPSocketAction } # spec name: tcpSocket
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiCoreV1Probe(;exec=nothing, failureThreshold=nothing, httpGet=nothing, initialDelaySeconds=nothing, periodSeconds=nothing, successThreshold=nothing, tcpSocket=nothing, timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Probe, Symbol("exec"), exec)
        setfield!(o, Symbol("exec"), exec)
        validate_property(IoK8sApiCoreV1Probe, Symbol("failureThreshold"), failureThreshold)
        setfield!(o, Symbol("failureThreshold"), failureThreshold)
        validate_property(IoK8sApiCoreV1Probe, Symbol("httpGet"), httpGet)
        setfield!(o, Symbol("httpGet"), httpGet)
        validate_property(IoK8sApiCoreV1Probe, Symbol("initialDelaySeconds"), initialDelaySeconds)
        setfield!(o, Symbol("initialDelaySeconds"), initialDelaySeconds)
        validate_property(IoK8sApiCoreV1Probe, Symbol("periodSeconds"), periodSeconds)
        setfield!(o, Symbol("periodSeconds"), periodSeconds)
        validate_property(IoK8sApiCoreV1Probe, Symbol("successThreshold"), successThreshold)
        setfield!(o, Symbol("successThreshold"), successThreshold)
        validate_property(IoK8sApiCoreV1Probe, Symbol("tcpSocket"), tcpSocket)
        setfield!(o, Symbol("tcpSocket"), tcpSocket)
        validate_property(IoK8sApiCoreV1Probe, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiCoreV1Probe

const _property_map_IoK8sApiCoreV1Probe = Dict{Symbol,Symbol}(Symbol("exec")=>Symbol("exec"), Symbol("failureThreshold")=>Symbol("failureThreshold"), Symbol("httpGet")=>Symbol("httpGet"), Symbol("initialDelaySeconds")=>Symbol("initialDelaySeconds"), Symbol("periodSeconds")=>Symbol("periodSeconds"), Symbol("successThreshold")=>Symbol("successThreshold"), Symbol("tcpSocket")=>Symbol("tcpSocket"), Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiCoreV1Probe = Dict{Symbol,String}(Symbol("exec")=>"IoK8sApiCoreV1ExecAction", Symbol("failureThreshold")=>"Int32", Symbol("httpGet")=>"IoK8sApiCoreV1HTTPGetAction", Symbol("initialDelaySeconds")=>"Int32", Symbol("periodSeconds")=>"Int32", Symbol("successThreshold")=>"Int32", Symbol("tcpSocket")=>"IoK8sApiCoreV1TCPSocketAction", Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiCoreV1Probe }) = collect(keys(_property_map_IoK8sApiCoreV1Probe))
Swagger.property_type(::Type{ IoK8sApiCoreV1Probe }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Probe[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Probe }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Probe[property_name]

function check_required(o::IoK8sApiCoreV1Probe)
    true
end

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


text:
using BinDeps

@BinDeps.setup

# define dependencies
libaswing = library_dependency("libaswing")

# define folders
srcdir = joinpath(BinDeps.depsdir(libaswing),"src")
libdir = joinpath(BinDeps.depsdir(libaswing),"usr","lib")
downloaddir = joinpath(BinDeps.depsdir(libaswing),"downloads")
aswingbin = joinpath(srcdir,"aswing","bin")
arpackdir = joinpath(srcdir,"ARPACK")
arpacklib = joinpath(arpackdir,"libarpack_aswing.a")

# get suffix
suffix = Sys.isapple() ? "dylib" : "so"

# create necessary paths
for path in [libdir,downloaddir,arpackdir]
  if !isfile(path)
    mkpath(path)
  end
end

# check if library exists and delete if necessary
if isfile(joinpath(libdir,"libaswing.$suffix"))
    rm(joinpath(libdir,"libaswing.$suffix"))
end

arpackbuildvars = [
  "home=$arpackdir",
  "PLAT=aswing",
  "FC=gfortran",
  "FFLAGS=-O3 -shared -fPIC",
  "MAKE=make",
  "PRECISIONS=complex16",
]

aswingbuildvars = [
  "INSTALL_DIR=$libdir",
  "SUFFIX=$suffix",
  "ARLIB=$arpacklib",
]

# build process
provides(SimpleBuild,
  (@build_steps begin
    ChangeDirectory(aswingbin)
    @build_steps begin
      ChangeDirectory(arpackdir)
      download_cmd("http://www.caam.rice.edu/software/ARPACK/SRC/arpack96.tar.gz",joinpath(downloaddir,"arpack96.tar.gz"))
      unpack_cmd(joinpath(downloaddir,"arpack96.tar.gz"),srcdir,".gz",".tar")
      download_cmd("http://www.caam.rice.edu/software/ARPACK/SRC/patch.tar.gz",joinpath(downloaddir,"patch.tar.gz"))
      unpack_cmd(joinpath(downloaddir,"arpack96.tar.gz"),srcdir,".gz",".tar")
      `make lib $arpackbuildvars`
    end
    `make $aswingbuildvars`
    `make install $aswingbuildvars`
  end),libaswing, os = :Unix)

@BinDeps.install Dict(:libaswing => :libaswing)


text:
using Colors
using FreeTypeAbstraction
using LaTeXStrings
using Luxor
using MathTeXEngine

"""
    latex_text_size(lstr::LaTeXString)
Returns the width and height of a latex string.
"""
function latex_text_size(lstr::LaTeXString)
    sentence = generate_tex_elements(lstr)
    els = filter(x -> x[1] isa TeXChar, sentence)
    chars = [x[1].char for x in els]
    fonts = [x[1].font for x in els]
    pos_x = [x[2][1] for x in els]
    pos_y = [x[2][2] for x in els]
    scales = [x[3] for x in els]
    extents = [FreeTypeAbstraction.get_extent(f, c) for (f, c) in zip(fonts, chars)]

    textw = []
    texth = []
    for i = 1:length(extents)
        textw = push!(
            textw,
            FreeTypeAbstraction.height_insensitive_boundingbox(
                extents[i],
                fonts[i],
            ).widths[1] * scales[i] + pos_x[i],
        )
        texth = push!(
            texth,
            FreeTypeAbstraction.height_insensitive_boundingbox(
                extents[i],
                fonts[i],
            ).widths[2] * scales[i] - pos_y[i],
        )
    end
    return (maximum(textw), maximum(texth))
end

"""
    Luxor.text(lstr::LaTeXString, valign=:baseline, halign=:left; kwargs...)
Draws LaTeX string using `MathTexEngine.jl`. Hence, uses ModernCMU as font family.
Note that `valign` is not perfect.
This function assumes that the axis are in the standard Luxor directions,
i.e. ↓→.
"""
function Luxor.text(
    lstr::LaTeXString,
    pt::Point;
    valign = :baseline,
    halign = :left,
    angle = 0,
)

    # Function from MathTexEngine
    sentence = generate_tex_elements(lstr)

    # Get current font size.
    textsize = get_fontsize()

    textw, texth = latex_text_size(lstr)

    if halign === :left
        translate_x = 0
    elseif halign === :right
        translate_x = -textw * textsize
    elseif halign === :center
        translate_x = -textw / 2 * textsize
    end

    if valign === :baseline
        translate_y = 0
    elseif valign === :top
        translate_y = -textsize
    elseif valign === :bottom
        translate_y = texth * textsize
    elseif valign === :middle
        translate_y = textsize / 2
    end

    # Writes text using ModernCMU font.
    for text in sentence
        if text[1] isa TeXChar
            @layer begin
                translate(translate_x, translate_y)
                fontface(text[1].font.family_name)
                fontsize(textsize * text[3])
                Luxor.text(
                    string(text[1].char),
                    pt + Point(text[2]...) * textsize * (1, -1),
                    angle = angle,
                )
            end
        elseif text[1] isa HLine
            @layer begin
                translate(translate_x, translate_y)
                pointstart = pt + Point(text[2]...) * textsize * (1, -1)
                pointend = pointstart + Point(text[1].width, 0) * textsize
                setline(1 * text[1].thickness * textsize)
                line(pointstart, pointend, :stroke)
            end
        end
    end
end

# Auxiliary drawing functions
function make_drawing(width, height, img_path, bkg_color, origin_p)
    d = Drawing(width, height, img_path)
    background(bkg_color)
    origin(origin_p)
    return d
end

cart2luxor(p) = Point(p.x, -p.y)
cart2luxor(p::Array) = Point(p[1], -p[2])
cart2luxor(x, y) = Point(x, -y)

luxor2cart(p::Point) = [p.x, -p.y]

parametric_point(t, x, y) = x + y * t

# Constants
## Canvas constants
width = 500
height = 500
path = "simple_diagram.png"
color = RGBA(0.0, 0.0, 0.0)
op = Point(width / 2, height / 2)

# Setting up drawing
my_draw = make_drawing(width, height, path, color, op)

sethue(RGBA(1.0, 1.0, 1.0))

p1 = [-125, 125] # Top left quadrant
p2 = [125, 125] # Top right quadrant
p3 = [-125, -125] # Bottom left quadrant
p4 = [125, -125] # Bottom right quadrant

function morphtext(
    label1::Any,
    label2::Any;
    offset1::Real = 20,
    offset2::Real = 20,
    angle = 0,
    size = nothing,
)
    fontsize(size)
    text(label1, cart2luxor(0, offset1), valign = :middle, halign = :center, angle = angle)
    text(label2, cart2luxor(0, -offset2), valign = :middle, halign = :center, angle = angle)
end

function morphtext(label::Any; offset::Real = 20, angle = 0, size = nothing)
    fontsize(size)
    text(label, cart2luxor(0, offset), valign = :middle, halign = :center, angle = angle)
end

function morphism(
    initial_point::Array,
    terminal_point::Array;
    label1::Any = "",
    label2::Any = "",
    offset1::Real = 10,
    offset2::Real = 10,
    size::Real = 16,
    padding::Real = 0.15,
)

     v⃗ = terminal_point .- initial_point
       θ = atan(v⃗[2] / v⃗[1]) + (v⃗[1] < 0 ? π : 0)

    padded_initial_point =
        between(cart2luxor(initial_point), cart2luxor(terminal_point), padding)
    padded_terminal_point =
        between(cart2luxor(initial_point), cart2luxor(terminal_point), 1 - padding)

    arrow(
        padded_initial_point,
        padded_terminal_point,
        arrowheadlength = 15,
        decoration = 0.5,
        decorate = () -> morphtext(
            label1,
            label2,
            offset1 = offset1,
            offset2 = offset2,
            angle = θ,
            size = size,
        ),
        linewidth = 3,
    )

end

function morphism(
    initial_point_1::Array,
    terminal_point_1::Array,
    initial_point_2::Array,
    terminal_point_2::Array;
    label1::Any = "",
    label2::Any = "",
    style1::String = "solid",
    style2::String = "solid",
    offset1::Real = 10,
    offset2::Real = 10,
    size::Real = 16,
    padding::Real = 0.15,
)

    Set([initial_point_1, terminal_point_1]) != Set([initial_point_2, terminal_point_2]) &&
        @error "Only two different points allowed"

    padded_is = []
    padded_ts = []
    angles = []
    for (i, t) in [[initial_point_1, terminal_point_1], [initial_point_2, terminal_point_2]]
         v⃗ = t .- i
           θ = atan(v⃗[2] / v⃗[1]) + (v⃗[1] < 0 ? π : 0)

        padded_i = between(cart2luxor(i), cart2luxor(t), padding)
        padded_t = between(cart2luxor(i), cart2luxor(t), 1 - padding)

         v⃗ₚ = luxor2cart(padded_i) .- luxor2cart(padded_t)
        x = 1
        y = 1
          if (-1 * v⃗ₚ[1] * x) / v⃗ₚ[2] |> isinf
              x = (v⃗ₚ[2] * y) / (-1 * v⃗ₚ[1])
        else
              y = (-1 * v⃗ₚ[1] * x) / v⃗ₚ[2]
        end
        i_new = parametric_point(
            length(padded_is) == 0 ? 10 : -10,
            luxor2cart(padded_i),
            [x; y],
        )
        t_new = parametric_point(
            length(padded_ts) == 0 ? 10 : -10,
            luxor2cart(padded_t),
            [x; y],
        )
        push!(padded_is, i_new)
        push!(padded_ts, t_new)
        push!(angles, θ)
    end

    setdash(style1)

    arrow(
        cart2luxor(padded_is[1]),
        cart2luxor(padded_ts[1]),
        arrowheadlength = 15,
        decoration = 0.5,
        decorate = () ->
            morphtext(label1, offset = offset1 + 20, angle = angles[1], size = size),
        linewidth = 3,
    )

    setdash(style2)

    arrow(
        cart2luxor(padded_is[2]),
        cart2luxor(padded_ts[2]),
        arrowheadlength = 15,
        decoration = 0.5,
        decorate = () ->
            morphtext(label2, offset = -offset2 - 10, angle = angles[2], size = size),
        linewidth = 3,
    )

    setdash("solid")

end

morphism( p1, p2, p2, p1, label1 = L"id_{A}", offset1 = -0, offset2 = 0, style1 = "longdashed", style2 = "dotted", padding = 0.15)
morphism(p4, p2, label1 = L"\pi", label2 = L"f", offset2 = 0, offset1 = 30, padding = 0.15)
morphism(p1, p3, label2 = L"f", offset2 = 20, padding = 0.15)
morphism(p3, p4, p3, p4, label1 = L"\phi_{12}", label2 = L"id_{B}", offset2 = 20, padding = 0.15)
morphism(p2, p3, label2 = L"f", offset2 = 20, padding = 0.15)

sethue("white")
fontsize(24)
text(L"A", cart2luxor(p1), valign = :middle, halign = :center)
text(L"A", cart2luxor(p2), valign = :middle, halign = :center)
text(L"B", cart2luxor(p3), valign = :middle, halign = :center)
text(L"b \in B", cart2luxor(p4), valign = :middle, halign = :center)

finish()


text:
function sc(x)
  if(abs(x)>1.0)
	a=trunc(log10(abs(x)))
	b=x/(10.0^a)
  else
	a=trunc(log10(abs(x)))-1
	b=x/(10.0^a)
  end
  return a,b
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Persistence
using Test

include("lib.jl")
include("transactions.jl")
include("persist.jl")