text:
############################################################################################
"""
Representation of a trajectory, ie a Hamiltonian with a discrete integrator that
also checks for divergence.
"""
struct TrajectoryNUTS{S<:Hamiltonian,F,C}
    "Hamiltonian"
    H::S
    "Log density of negative log energy at initial point."
    ℓH₀::F
    "Stepsize for leapfrog."
    ϵ::F
    "Maximum tree depth."
    max_depth::Int64
    "Smallest decrease allowed in the log density."
    min_Δ::Float64
    "Turn statistic configuration."
    turn_statistic_configuration::C
    function TrajectoryNUTS(
        H::S,
        ℓH₀::F,
        ϵ::F,
        max_depth=DEFAULT_MAX_TREE_DEPTH,
        min_Δ=-1000.0,
        turn_statistic_configuration=Val{:generalized}(),
    ) where {S<:Hamiltonian,F<:AbstractFloat}
        @argcheck 0 < max_depth ≤ MAX_DIRECTIONS_DEPTH
        @argcheck min_Δ < 0
        C = typeof(turn_statistic_configuration)
        return new{S,F,C}(H, ℓH₀, ϵ, max_depth, min_Δ, turn_statistic_configuration)
    end
end

"""
$(SIGNATURES)
Move along the trajectory in the specified direction. Return the new position.

# Examples
```julia
```

"""
function move(trajectory::T, phasepoint::PhasePoint, fwd::Bool) where {T<:TrajectoryNUTS}
    @unpack H, ϵ = trajectory
    return leapfrog(H, phasepoint, fwd ? ϵ : -ϵ)
end

############################################################################################
# Proposal
"""
$(SIGNATURES)
Given (relative) log probabilities `ω₁` and `ω₂`, return the log probabiliy of drawing a sample from the second (`logprob2`).
When `bias`, biases towards the second argument, introducing anti-correlations.

# Examples
```julia
```

"""
function biased_progressive_logprob2(
    bias::Bool, ω₁::T, ω₂::T, ω=BaytesCore.logaddexp(ω₁, ω₂)
) where {T<:Real}
    return ω₂ - (bias ? ω₁ : ω)
end

"""
$(SIGNATURES)
Random boolean which is `true` with the given probability `exp(logprob)`, which can be `≥ 1`
in which case no random value is drawn.

# Examples
```julia
```

"""
function rand_bool_logprob(_rng::Random.AbstractRNG, logprob::T) where {T<:Real}
    return logprob ≥ 0 || (randexp(_rng, T) > -logprob)
end

"""
$(SIGNATURES)
Calculate the log probability if selecting the subtree corresponding to `ω₂`. Being the log
of a probability, it is always `≤ 0`, but implementations are allowed to return and accept
values `> 0` and treat them as `0`. When `is_doubling`, the tree corresponding to `ω₂` was obtained from a doubling step (this
can be relevant eg for biased progressive sampling). The value `ω = logaddexp(ω₁, ω₂)` is provided for avoiding redundant calculations.
See [`biased_progressive_logprob2`](@ref) for an implementation.

# Examples
```julia
```

"""
function calculate_logprob2(::TrajectoryNUTS, is_doubling, ω₁, ω₂, ω)
    return biased_progressive_logprob2(is_doubling, ω₁, ω₂, ω)
end

"""
$(SIGNATURES)
Combine two proposals `ζ₁, ζ₂` on `trajectory`, with log probability `logprob2` for
selecting `ζ₂`.
 `ζ₁` is before `ζ₂` iff `is_forward`.

# Examples
```julia
```

"""
function combine_proposals(
    _rng::Random.AbstractRNG,
    trajectory::TrajectoryNUTS,
    z₁::PhasePoint,
    z₂::PhasePoint,
    logprob2::T,
    is_forward,
) where {T<:Real}
    return rand_bool_logprob(_rng, logprob2) ? z₂ : z₁
end

############################################################################################
# statistics for visited nodes
struct AcceptanceStatistic{F<:AbstractFloat}
    "
    Logarithm of the sum of metropolis acceptances probabilities over the whole trajectory
    (including invalid parts).
    "
    log_sum_α::F
    "Total number of leapfrog steps."
    steps::Int64
end

function combine_acceptance_statistics(A::AcceptanceStatistic, B::AcceptanceStatistic)
    return AcceptanceStatistic(
        BaytesCore.logaddexp(A.log_sum_α, B.log_sum_α), A.steps + B.steps
    )
end

"""
$(SIGNATURES)
Acceptance statistic for a leaf. The initial leaf is considered not to be visited.

# Examples
```julia
```

"""
function leaf_acceptance_statistic(Δ, is_initial)
    return if is_initial
        AcceptanceStatistic(oftype(Δ, -Inf), 0)
    else
        AcceptanceStatistic(min(Δ, 0), 1)
    end
end

"""
$(SIGNATURES)
Return the acceptance rate (a `Real` betwen `0` and `1`).

# Examples
```julia
```

"""
acceptance_rate(A::AcceptanceStatistic) = min(exp(A.log_sum_α) / A.steps, 1)

"""
$(SIGNATURES)
Combine visited node statistics for adjacent trees trajectory. Implementation should be
invariant to the ordering of `v₁` and `v₂` (ie the operation is commutative).

# Examples
```julia
```

"""
combine_visited_statistics(::TrajectoryNUTS, v, w) = combine_acceptance_statistics(v, w)

############################################################################################
# abstract trajectory interface
"""
$(TYPEDEF)
Information about an invalid (sub)tree, using positions relative to the starting node.
1. When `left < right`, this tree was *turning*.
2. When `left == right`, this is a *divergent* node.
3. `left == 1 && right == 0` is used as a sentinel value for reaching maximum depth without
encountering any invalid trees (see [`REACHED_MAX_DEPTH`](@ref). All other `left > right`
values are disallowed.

# Fields
$(TYPEDFIELDS)
"""
struct InvalidTree
    left::Int64
    right::Int64
end
InvalidTree(i::Integer) = InvalidTree(i, i)
is_divergent(invalid_tree::InvalidTree) = invalid_tree.left == invalid_tree.right

############################################################################################
# Directions

"""
$(TYPEDEF)
Internal type implementing random directions.

# Fields
$(TYPEDFIELDS)
"""
struct Directions
    flags::UInt32
end

Base.rand(_rng::Random.AbstractRNG, ::Type{Directions}) = Directions(rand(_rng, UInt32))

"""
$(SIGNATURES)
Return the next direction flag and the new state of directions. Results are undefined for
more than [`MAX_DIRECTIONS_DEPTH`](@ref) updates.

# Examples
```julia
```

"""
function next_direction(directions::Directions)
    @unpack flags = directions
    return Bool(flags & 0x01), Directions(flags >>> 1)
end

"Sentinel value for reaching maximum depth."
const REACHED_MAX_DEPTH = InvalidTree(1, 0)

################################################################################
# turn analysis

"""
$(TYPEDEF)
Statistics for the identification of turning points. See Betancourt (2017, appendix), and
subsequent discussion of improvements at
<https://discourse.mc-stan.org/t/nuts-misses-u-turns-runs-in-circles-until-max-treedepth/9727/>.
Momenta `p₋` and `p₊` are kept so that they can be added to `ρ` when combining turn
statistics.
Turn detection is always done by [`combine_turn_statistics`](@ref), which returns `nothing`
in case of turning. A `GeneralizedTurnStatistic` should always correspond to a trajectory
that is *not* turning (or a leaf node, where the concept does not apply).

# Fields
$(TYPEDFIELDS)
"""
struct GeneralizedTurnStatistic{T}
    "momentum at the left edge of the trajectory"
    ρ₋::T
    "p♯ at the left edge of the trajectory"
    ρ♯₋::T
    "momentum at the right edge of the trajectory"
    ρ₊::T
    "p♯ at the right edge of the trajectory"
    ρ♯₊::T
    "sum of momenta along trajectory"
    ρₛ::T
end

function leaf_turn_statistic(::Val{:generalized}, H::Hamiltonian, phasepoint::PhasePoint)
    @unpack ρ = phasepoint
    ρ♯ = calculate_ρ♯(H, phasepoint)
    return GeneralizedTurnStatistic(ρ, ρ♯, ρ, ρ♯, ρ)
end

"""
Internal test for turning. See Betancourt (2017, appendix).
"""
_is_turning(ρ♯₋, ρ♯₊, ρ) = dot(ρ♯₋, ρ) < 0 || dot(ρ♯₊, ρ) < 0

"""
$(SIGNATURES)
Combine turn statistics on trajectory. Implementation can assume that the trees that
correspond to the turn statistics have the same ordering.
When
```julia
τ = combine_turn_statistics(trajectory, τ₁, τ₂)
is_turning(trajectory, τ)
```
the combined turn statistic `τ` is guaranteed not to escape the caller, so it can eg change
type.

# Examples
```julia
```

"""
function combine_turn_statistics(
    trajectory::TrajectoryNUTS, x::GeneralizedTurnStatistic, y::GeneralizedTurnStatistic
)
    _is_turning(x.ρ♯₋, y.ρ♯₋, x.ρₛ + y.ρ₋) && return nothing
    _is_turning(x.ρ♯₊, y.ρ♯₊, x.ρ₊ + y.ρₛ) && return nothing
    ρₛᵤₘ = x.ρₛ + y.ρₛ
    #!NOTE: Most left and Right Momenta
    _is_turning(x.ρ♯₋, y.ρ♯₊, ρₛᵤₘ) && return nothing
    return GeneralizedTurnStatistic(x.ρ₋, x.ρ♯₋, y.ρ₊, y.ρ♯₊, ρₛᵤₘ)
end

"""
Internal test for turning. See Betancourt (2017, appendix).
"""
is_turning(::TrajectoryNUTS, ::GeneralizedTurnStatistic) = false
is_turning(::TrajectoryNUTS, ::Nothing) = true

############################################################################################
# utilities

"""
$(SIGNATURES)
Combine turn statistics with the given direction. When `is_forward`, `τ₁` is before `τ₂`,
otherwise after. Internal helper function.

# Examples
```julia
```

"""
@inline function combine_turn_statistics_in_direction(
    trajectory::TrajectoryNUTS,
    τ₁::GeneralizedTurnStatistic,
    τ₂::GeneralizedTurnStatistic,
    is_forward::Bool,
)
    if is_forward
        combine_turn_statistics(trajectory, τ₁, τ₂)
    else
        combine_turn_statistics(trajectory, τ₂, τ₁)
    end
end

function combine_proposals_and_logweights(
    _rng::Random.AbstractRNG,
    trajectory::TrajectoryNUTS,
    ζ₁::PhasePoint,
    ζ₂::PhasePoint,
    ω₁::R,
    ω₂::R,
    is_forward::Bool,
    is_doubling::Bool,
) where {R<:Real}
    ω = BaytesCore.logaddexp(ω₁, ω₂)
    logprob2 = calculate_logprob2(trajectory, is_doubling, ω₁, ω₂, ω)
    ζ = combine_proposals(_rng, trajectory, ζ₁, ζ₂, logprob2, is_forward)
    return ζ, ω
end

############################################################################################
# leafs

"""
$(SIGNATURES)
Information for a tree made of a single node. When `is_initial == true`, this is the first
node.
The first value is either
1. `nothing` for a divergent node,
2. a tuple containing the proposal `ζ`, the log weight (probability) of the node `ω`, the
turn statistics `τ` (never tested with `is_turning` for leaf nodes).
The second value is the visited node information.

# Examples
```julia
```

"""
function leaf(trajectory::TrajectoryNUTS, phasepoint::PhasePoint, is_initial)
    @unpack H, ℓH₀, min_Δ, turn_statistic_configuration = trajectory
    Δ = is_initial ? zero(ℓH₀) : ℓdensity(trajectory.H, phasepoint) - ℓH₀
    isdiv = Δ < min_Δ
    v = leaf_acceptance_statistic(Δ, is_initial)
    if isdiv
        nothing, v
    else
        τ = leaf_turn_statistic(turn_statistic_configuration, H, phasepoint)
        (phasepoint, Δ, τ), v
    end
end

############################################################################################
"""
$(SIGNATURES)
Traverse the tree of given `depth` adjacent to point `z` in `trajectory`.
`is_forward` specifies the direction, `rng` is used for random numbers in
[`combine_proposals`](@ref). `i` is an integer position relative to the initial node (`0`).
The *first value* is either
1. an `InvalidTree`, indicating the first divergent node or turning subtree that was
encounteted and invalidated this tree.
2. a tuple of `(ζ, ω, τ, z′, i′), with
    - `ζ`: the proposal from the tree.
    - `ω`: the log weight of the subtree that corresponds to the proposal
    - `τ`: turn statistics
    - `z′`: the last node of the tree
    - `i′`: the position of the last node relative to the initial node.
The *second value* is always the visited node statistic.

# Examples
```julia
```

"""
function adjacent_tree(
    _rng::Random.AbstractRNG,
    trajectory::TrajectoryNUTS,
    z::PhasePoint,
    i::Int64,
    depth::Int64,
    is_forward::Bool,
)
    i′ = i + (is_forward ? 1 : -1)
    if depth == 0
        z′ = move(trajectory, z, is_forward)
        ζωτ, v = leaf(trajectory, z′, false)
        if ζωτ ≡ nothing
            InvalidTree(i′), v
        else
            (ζωτ..., z′, i′), v
        end
    else
        ## “left” tree
        t₋, v₋ = adjacent_tree(_rng, trajectory, z, i, depth - 1, is_forward)
        t₋ isa InvalidTree && return t₋, v₋
        ζ₋, ω₋, τ₋, z₋, i₋ = t₋

        ## “right” tree — visited information from left is kept even if invalid
        t₊, v₊ = adjacent_tree(_rng, trajectory, z₋, i₋, depth - 1, is_forward)
        v = combine_visited_statistics(trajectory, v₋, v₊)
        t₊ isa InvalidTree && return t₊, v
        ζ₊, ω₊, τ₊, z₊, i₊ = t₊

        ## turning invalidates
        τ = combine_turn_statistics_in_direction(trajectory, τ₋, τ₊, is_forward)
        is_turning(trajectory, τ) && return InvalidTree(i′, i₊), v

        ## valid subtree, combine proposals
        ζ, ω = combine_proposals_and_logweights(
            _rng, trajectory, ζ₋, ζ₊, ω₋, ω₊, is_forward, false
        )
        (ζ, ω, τ, z₊, i₊), v
    end
end

"""
$(SIGNATURES)
Sample a `trajectory` starting at `z`, up to `max_depth`. `directions` determines the tree
expansion directions.
Return the following values
- `ζ`: proposal from the tree
- `v`: visited node statistics
- `termination`: an `InvalidTree` (this includes the last doubling step turning, which is
  technically a valid tree) or `REACHED_MAX_DEPTH` when all subtrees were valid and no
  turning happens.
- `depth`: the depth of the tree that was sampled from. Doubling steps that lead to an
  invalid adjacent tree do not contribute to `depth`.

# Examples
```julia
```

"""
function sample_trajectory(
    _rng::Random.AbstractRNG,
    trajectory::TrajectoryNUTS,
    z::PhasePoint,
    directions::Directions,
)
    @unpack max_depth = trajectory
    @argcheck max_depth ≤ MAX_DIRECTIONS_DEPTH
    (ζ, ω, τ), v = leaf(trajectory, z, true)
    z₋ = z₊ = z
    depth = 0
    termination = REACHED_MAX_DEPTH
    i₋ = i₊ = 0
    while depth < max_depth
        is_forward, directions = next_direction(directions)
        t′, v′ = adjacent_tree(
            _rng, trajectory, is_forward ? z₊ : z₋, is_forward ? i₊ : i₋, depth, is_forward
        )
        v = combine_visited_statistics(trajectory, v, v′)

        ## invalid adjacent tree: stop
        t′ isa InvalidTree && (termination = t′; break)

        ## extract information from adjacent tree
        ζ′, ω′, τ′, z′, i′ = t′

        ## update edges and combine proposals
        if is_forward
            z₊, i₊ = z′, i′
        else
            z₋, i₋ = z′, i′
        end

        ## tree has doubled successfully
        ζ, ω = combine_proposals_and_logweights(
            _rng, trajectory, ζ, ζ′, ω, ω′, is_forward, true
        )
        depth += 1

        ## when the combined tree is turning, stop
        τ = combine_turn_statistics_in_direction(trajectory, τ, τ′, is_forward)
        is_turning(trajectory, τ) && (termination = InvalidTree(i₋, i₊);
        break)
    end
    return ζ, v, termination, depth
end

############################################################################################
#export
export TrajectoryNUTS,
    biased_progressive_logprob2,
    rand_bool_logprob,
    calculate_logprob2,
    combine_proposals,
    AcceptanceStatistic,
    combine_acceptance_statistics,
    leaf_acceptance_statistic,
    acceptance_rate,
    combine_visited_statistics,
    InvalidTree,
    is_divergent,
    Directions,
    next_direction,
    REACHED_MAX_DEPTH,
    GeneralizedTurnStatistic,
    leaf_turn_statistic,
    _is_turning,
    combine_turn_statistics,
    combine_turn_statistics_in_direction,
    combine_proposals_and_logweights,
    leaf,
    adjacent_tree,
    sample_trajectory

ext: jl


text:
"""
    InverseGaussian(μ,λ)

The *inverse Gaussian distribution* with mean `μ` and shape `λ` has probability density function

```math
f(x; \\mu, \\lambda) = \\sqrt{\\frac{\\lambda}{2\\pi x^3}}
\\exp\\!\\left(\\frac{-\\lambda(x-\\mu)^2}{2\\mu^2x}\\right), \\quad x > 0
```

```julia
InverseGaussian()              # Inverse Gaussian distribution with unit mean and unit shape, i.e. InverseGaussian(1, 1)
InverseGaussian(μ),            # Inverse Gaussian distribution with mean μ and unit shape, i.e. InverseGaussian(μ, 1)
InverseGaussian(μ, λ)          # Inverse Gaussian distribution with mean μ and shape λ

params(d)           # Get the parameters, i.e. (μ, λ)
mean(d)             # Get the mean parameter, i.e. μ
shape(d)            # Get the shape parameter, i.e. λ
```

External links

* [Inverse Gaussian distribution on Wikipedia](http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution)

"""
struct InverseGaussian{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    λ::T
    InverseGaussian{T}(μ::T, λ::T) where {T<:Real} = new{T}(μ, λ)
end

function InverseGaussian(μ::T, λ::T; check_args::Bool=true) where {T<:Real}
    @check_args InverseGaussian (μ, μ > zero(μ)) (λ, λ > zero(λ))
    return InverseGaussian{T}(μ, λ)
end

InverseGaussian(μ::Real, λ::Real; check_args::Bool=true) = InverseGaussian(promote(μ, λ)...; check_args=check_args)
InverseGaussian(μ::Integer, λ::Integer; check_args::Bool=true) = InverseGaussian(float(μ), float(λ); check_args=check_args)
InverseGaussian(μ::Real; check_args::Bool=true) = InverseGaussian(μ, one(μ); check_args=check_args)
InverseGaussian() = InverseGaussian{Float64}(1.0, 1.0)

@distr_support InverseGaussian 0.0 Inf

#### Conversions

function convert(::Type{InverseGaussian{T}}, μ::S, λ::S) where {T <: Real, S <: Real}
    InverseGaussian(T(μ), T(λ))
end
function convert(::Type{InverseGaussian{T}}, d::InverseGaussian{S}) where {T <: Real, S <: Real}
    InverseGaussian(T(d.μ), T(d.λ), check_args=false)
end

#### Parameters

shape(d::InverseGaussian) = d.λ
params(d::InverseGaussian) = (d.μ, d.λ)
partype(::InverseGaussian{T}) where {T} = T

#### Statistics

mean(d::InverseGaussian) = d.μ

var(d::InverseGaussian) = d.μ^3 / d.λ

skewness(d::InverseGaussian) = 3sqrt(d.μ / d.λ)

kurtosis(d::InverseGaussian) = 15d.μ / d.λ

function mode(d::InverseGaussian)
    μ, λ = params(d)
    r = μ / λ
    μ * (sqrt(1 + (3r/2)^2) - (3r/2))
end

#### Evaluation

function pdf(d::InverseGaussian{T}, x::Real) where T<:Real
    if x > 0
        μ, λ = params(d)
        return sqrt(λ / (twoπ * x^3)) * exp(-λ * (x - μ)^2 / (2μ^2 * x))
    else
        return zero(T)
    end
end

function logpdf(d::InverseGaussian{T}, x::Real) where T<:Real
    if x > 0
        μ, λ = params(d)
        return (log(λ) - (log2π + 3log(x)) - λ * (x - μ)^2 / (μ^2 * x))/2
    else
        return -T(Inf)
    end
end

function cdf(d::InverseGaussian, x::Real)
    μ, λ = params(d)
    y = max(x, 0)
    u = sqrt(λ / y)
    v = y / μ
    z = normcdf(u * (v - 1)) + exp(2λ / μ) * normcdf(-u * (v + 1))

    # otherwise `NaN` is returned for `+Inf`
    return isinf(x) && x > 0 ? one(z) : z
end

function ccdf(d::InverseGaussian, x::Real)
    μ, λ = params(d)
    y = max(x, 0)
    u = sqrt(λ / y)
    v = y / μ
    z = normccdf(u * (v - 1)) - exp(2λ / μ) * normcdf(-u * (v + 1))

    # otherwise `NaN` is returned for `+Inf`
    return isinf(x) && x > 0 ? zero(z) : z
end

function logcdf(d::InverseGaussian, x::Real)
    μ, λ = params(d)
    y = max(x, 0)
    u = sqrt(λ / y)
    v = y / μ

    a = normlogcdf(u * (v - 1))
    b = 2λ / μ + normlogcdf(-u * (v + 1))
    z = logaddexp(a, b)

    # otherwise `NaN` is returned for `+Inf`
    return isinf(x) && x > 0 ? zero(z) : z
end

function logccdf(d::InverseGaussian, x::Real)
    μ, λ = params(d)
    y = max(x, 0)
    u = sqrt(λ / y)
    v = y / μ

    a = normlogccdf(u * (v - 1))
    b = 2λ / μ + normlogcdf(-u * (v + 1))
    z = logsubexp(a, b)

    # otherwise `NaN` is returned for `+Inf`
    return isinf(x) && x > 0 ? oftype(z, -Inf) : z
end

@quantile_newton InverseGaussian

#### Sampling

# rand method from:
#   John R. Michael, William R. Schucany and Roy W. Haas (1976)
#   Generating Random Variates Using Transformations with Multiple Roots
#   The American Statistician , Vol. 30, No. 2, pp. 88-90
function rand(rng::AbstractRNG, d::InverseGaussian)
    μ, λ = params(d)
    z = randn(rng)
    v = z * z
    w = μ * v
    x1 = μ + μ / (2λ) * (w - sqrt(w * (4λ + w)))
    p1 = μ / (μ + x1)
    u = rand(rng)
    u >= p1 ? μ^2 / x1 : x1
end

#### Fit model

"""
Sufficient statistics for `InverseGaussian`, containing the weighted
sum of observations, the weighted sum of inverse points and sum of weights.
"""
struct InverseGaussianStats <: SufficientStats
    sx::Float64      # (weighted) sum of x
    sinvx::Float64   # (weighted) sum of 1/x
    sw::Float64      # sum of sample weight
end

function suffstats(::Type{<:InverseGaussian}, x::AbstractVector{<:Real})
    sx = sum(x)
    sinvx = sum(inv, x)
    InverseGaussianStats(sx, sinvx, length(x))
end

function suffstats(::Type{<:InverseGaussian}, x::AbstractVector{<:Real}, w::AbstractVector{<:Real})
    n = length(x)
    if length(w) != n
        throw(DimensionMismatch("Inconsistent argument dimensions."))
    end
    T = promote_type(eltype(x), eltype(w))
    sx = zero(T)
    sinvx = zero(T)
    sw = zero(T)
    @inbounds @simd for i in eachindex(x)
        sx += w[i]*x[i]
        sinvx += w[i]/x[i]
        sw += w[i]
    end
    InverseGaussianStats(sx, sinvx, sw)
end

function fit_mle(::Type{<:InverseGaussian}, ss::InverseGaussianStats)
    mu = ss.sx / ss.sw
    invlambda = ss.sinvx / ss.sw  -  inv(mu)
    InverseGaussian(mu, inv(invlambda))
end


text:
"""
BlobLoG stores information about the location of peaks as discovered by `blob_LoG`.
It has fields:

- location: the location of a peak in the filtered image (a CartesianIndex)
- σ: the value of σ which lead to the largest `-LoG`-filtered amplitude at this location
- amplitude: the value of the `-LoG(σ)`-filtered image at the peak

Note that the radius is equal to σ√2.

See also: [`blob_LoG`](@ref).
"""
struct BlobLoG{T,S,N}
    location::CartesianIndex{N}
    σ::S
    amplitude::T
end
BlobLoG(; location, σ, amplitude) = BlobLoG(location, σ, amplitude)

function Base.show(io::IO, bl::BlobLoG)
    print(io, "BlobLoG(location=", bl.location, ", σ=", bl.σ, ", amplitude=", bl.amplitude, ")")
end

"""
    blob_LoG(img, σscales; edges::Tuple=(true, false, ...), σshape::Tuple=(1, ...), rthresh=0.001) -> Vector{BlobLoG}

Find "blobs" in an N-D image using the negative Lapacian of Gaussians
with the specifed vector or tuple of σ values. The algorithm searches for places
where the filtered image (for a particular σ) is at a peak compared to all
spatially- and σ-adjacent voxels, where σ is `σscales[i] * σshape` for some i.
By default, `σshape` is an ntuple of 1s.

The optional `edges` argument controls whether peaks on the edges are
included. `edges` can be `true` or `false`, or a N+1-tuple in which
the first entry controls whether edge-σ values are eligible to serve
as peaks, and the remaining N entries control each of the N dimensions
of `img`.

`rthresh` controls the minimum amplitude of peaks in the -LoG-filtered image,
as a fraction of `maximum(abs, img)` and the volume of the Gaussian.

# Examples

While most images are 2- or 3-dimensional, it will be easier to illustrate this with
a one-dimensional "image" containing two Gaussian blobs of different sizes:

```jldoctest; setup=:(using ImageFiltering), filter=r"amplitude=.*"]
julia> σs = 2.0.^(1:6);

julia> img = zeros(100); img[20:30] = [exp(-x^2/(2*4^2)) for x=-5:5]; img[50:80] = [exp(-x^2/(2*8^2)) for x=-15:15];

julia> blob_LoG(img, σs; edges=false)
2-element Vector{BlobLoG{Float64, Tuple{Float64}, 1}}:
 location=CartesianIndex(25,), σ=(4.0,), amplitude=0.10453155018303673
 location=CartesianIndex(65,), σ=(8.0,), amplitude=0.046175719034527364
```

The other two are centered in their corresponding "features," and the width `σ`
reflects the width of the feature itself.

`blob_LoG` tends to work best for shapes that are "Gaussian-like" but does
generalize somewhat.

# Citation:

Lindeberg T (1998), "Feature Detection with Automatic Scale Selection",
International Journal of Computer Vision, 30(2), 79–116.

See also: [`BlobLoG`](@ref).
"""
function blob_LoG(img::AbstractArray{T,N}, σscales;
                  edges::Union{Bool,Tuple{Bool,Vararg{Bool,N}}}=(true, ntuple(d->false, Val(N))...),
                  σshape::NTuple{N,Real}=ntuple(d->1, Val(N)),
                  rthresh::Real=1//1000) where {T<:Union{AbstractGray,Real},N}
    if edges isa Bool
        edges = (edges, ntuple(d->edges,Val(N))...)
    end
    sigmas = sort(σscales)
    img_LoG = multiLoG(img, sigmas, σshape)
    maxima = findlocalmaxima(img_LoG; edges=edges)
    # The "density" should not be much smaller than 1/volume of the Gaussian
    if !iszero(rthresh)
        athresh = rthresh./(sigmas.^N .* prod(σshape))
        imgmax = maximum(abs, img)
        [BlobLoG(CartesianIndex(tail(x.I)), sigmas[x[1]].*σshape, img_LoG[x]) for x in maxima if img_LoG[x] > athresh[x[1]]*imgmax]
    else
        [BlobLoG(CartesianIndex(tail(x.I)), sigmas[x[1]].*σshape, img_LoG[x]) for x in maxima]
    end
end

function multiLoG(img::AbstractArray{T,N}, sigmas, σshape) where {T,N}
    issorted(sigmas) || error("sigmas must be sorted")
    img_LoG = similar(img, float(eltype(T)), (Base.OneTo(length(sigmas)), axes(img)...))
    colons = ntuple(d->Colon(), Val(N))
    @inbounds for (isigma, σ) in enumerate(sigmas)
        LoG_slice = @view img_LoG[isigma, colons...]
        imfilter!(LoG_slice, img, Kernel.LoG(ntuple(i->σ*σshape[i], Val(N))), "reflect")
        LoG_slice .*= -σ
    end
    return img_LoG
end

default_window(img) = (cs = coords_spatial(img); ntuple(d -> d ∈ cs ? 3 : 1, ndims(img)))

"""
    findlocalmaxima(img; window=default_window(img), edges=true) -> Vector{CartesianIndex}

Returns the coordinates of elements whose value is larger than all of
their immediate neighbors.  `edges` is a boolean specifying whether to include the
first and last elements of each dimension, or a tuple-of-Bool
specifying edge behavior for each dimension separately.

The `default_window` is 3 for each spatial dimension of `img`, and 1 otherwise, implying
that maxima are detected over nearest-neighbors in each spatial "slice" by default.
"""
findlocalmaxima(img::AbstractArray; window=default_window(img), edges=true) =
    findlocalextrema(>, img, window, edges)

"""
    findlocalminima(img; window=default_window(img), edges=true) -> Vector{CartesianIndex}

Like [`findlocalmaxima`](@ref), but returns the coordinates of the smallest elements.
"""
findlocalminima(img::AbstractArray; window=default_window(img), edges=true) =
    findlocalextrema(<, img, window, edges)

findlocalextrema(f, img::AbstractArray{T,N}, window, edges::Bool) where {T,N} = findlocalextrema(f, img, window, ntuple(d->edges,Val(N)))

function findlocalextrema(f::F, img::AbstractArray{T,N}, window::Dims{N}, edges::NTuple{N,Bool}) where {F,T<:Union{Gray,Number},N}
    extrema = Vector{CartesianIndex{N}}(undef, 0)
    Iedge = CartesianIndex(map(!, edges))
    R0 = CartesianIndices(img)
    R = clippedinds(R0, Iedge)
    halfwindow = CartesianIndex(map(x -> x >> 1, window))
    Rinterior = clippedinds(R0, halfwindow)
    Rwindow = _colon(-halfwindow, halfwindow)
    z = zero(halfwindow)
    for i in R
        isextrema = true
        img_i = img[i]
        if i ∈ Rinterior
            # If i is in the interior, we don't have to worry about i+j being out-of-bounds
            for j in Rwindow
                j == z && continue
                if !f(img_i, img[i+j])
                    isextrema = false
                    break
                end
            end
        else
            for j in Rwindow
                (j == z || i+j ∉ R0) && continue
                if !f(img_i, img[i+j])
                    isextrema = false
                    break
                end
            end
        end
        isextrema && push!(extrema, i)
    end
    extrema
end

clippedinds(Router, Iclip) = _colon(first(Router)+Iclip, last(Router)-Iclip)


text:
"""
A struct that represents an agent storing it's accuracy and balance.
The accurancy is set up once during the agent creation and remain the same during the simulation.
The balance is updated accoring to agent's perfromance in the curation process.
"""
mutable struct Agent{TI<:Integer, TF<:Float64}
    accuracy::TI
    balance::TF
end

# constructor
Agent(a::Integer) = Agent(min(100,a), 50.0)

"""
Creates a group of X (count) agents all of them having the same accuracy (acc)
"""
function setupAgentsWithFixedAccuracy(count, acc)
    Agent[Agent(acc) for i in 1:count]
end

"""
Creates a group of X (count) agents all with accuracy distributed according to the
normal distribution with passed parameters (mean and std)
"""
function setupRandomAgents(count, mean, std)
    agents = Agent[Agent(trunc(Int, rand(TruncatedNormal(mean, std, 0, 100)))) for i in 1:count]
end

function setupMixedAgents(lowCount, lowAcc, highCount, highAcc)
    agents = vcat(setupAgentsWithFixedAccuracy(lowCount, lowAcc), setupAgentsWithFixedAccuracy(highCount, highAcc))
end

"""
This function represents the perception of an item by the agent.
Agents with lower accuracy tend to overrate an item.
This variation of the function rates the item with a fixed positive bias.
"""
function evaluateCandidateByAgentFixed(candidate, agent)
    evaluation = candidate + (100-agent.accuracy)/3
end

"""
This function represents the perception of an item by the agent.
Agents with lower accuracy tend to overrate an item.
This variation of the function rates the item with a randomized bias based on the normal distribution.
"""
function evaluateCandidateByAgent(candidate, agent)
    if (agent.accuracy == 100)
        candidate
    else
        evaluation = rand(Normal(candidate, (100-agent.accuracy)/3))
        candidate + abs(evaluation - candidate)
    end
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
const type_to_label = Dict{String,String}([
    "article" => "journal",
    "book" => "book",
    "booklet" => "booklet",
    "eprint" => "eprint",
    "inbook" => "book chapter",
    "incollection" => "book section",
    "inproceedings" => "conference",
    "manual" => "manual",
    "mastersthesis" => "master's thesis",
    "misc" => "other",
    "phdthesis" => "doctoral thesis",
    "proceedings" => "proceedings",
    "techreport" => "report",
    "unpublished" => "other",
])

"""
    xtitle(entry

Format the title of an `Entry` for web export.
"""
function xtitle(entry)
    return isdefined(entry, :title) ? entry.title : get(entry.fields, "title", "")
end

"""
    xnames(entry, editors = false; names = :full)

Format the name of an `Entry` for web export.

# Arguments:
- `entry`: an entry
- `editors`: `true` if the name describes editors
- `names`: :full (full names) or :last (last names + first name abbreviation)
"""
function xnames(
    entry,
    editors=false;
    names=:full, # Current options: :last, :full
)
    # forces the names to be editors' name if the entry are Proceedings
    if !editors && entry.type ∈ ["proceedings"]
        return xnames(entry, true)
    end
    entry_names = editors ? entry.editors : entry.authors

    if names == :last
        parts = map(s -> [s.particle, s.last, s.junior], entry_names)
    else
        parts = map(s -> [s.first, s.middle, s.particle, s.last, s.junior], entry_names)
    end

    entry_names = map(parts) do s
        return join(filter(!isempty, s), " ")
    end
    str = join(entry_names, ", ")
    return replace(str, r"[\n\r ]+" => " ") # TODO: make it cleaner (is replace still necessary)
end

"""
    xin(entry)

Format the appears-`in` field of an `Entry` for web export.
"""
function xin(entry)
    str = ""
    temp_last = false
    temp = []
    if entry.type == "article"
        temp = [
            entry.in.journal,
            entry.in.volume * (entry.in.number != "" ? "($(entry.in.number))" : ""),
            entry.in.pages,
            entry.date.year,
        ]
    elseif entry.type == "book"
        temp = [entry.in.publisher, entry.in.address, entry.date.year]
    elseif entry.type == "booklet"
        temp = [entry.access.howpublished, entry.date.year]
    elseif entry.type == "eprint"
        if isempty(entry.eprint.archive_prefix)
            str *= entry.eprint.eprint
        else
            str *= "$(entry.eprint.archive_prefix):$(entry.eprint.eprint) [$(entry.eprint.primary_class)]"
        end
    elseif entry.type == "inbook"
        temp = [
            entry.booktitle,
            isempty(entry.in.chapter) ? entry.in.pages : entry.in.chapter,
            entry.in.publisher,
            entry.in.address,
        ]
    elseif entry.type == "incollection"
        # TODO: check if this new or the old format is/was correct, that "editors" seems out of place (and the title was switched with the names)?
        temp = [
            "In $(entry.booktitle)",
            "editors",
            xnames(entry, true),
            entry.in.pages * ". " * entry.in.publisher, # TODO: conditional ". " if one of the strings is empty?
            entry.in.address,
            entry.date.year,
        ]
    elseif entry.type == "inproceedings"
        temp_last = entry.in.publisher != ""
        temp = [
            " In " * entry.booktitle,
            entry.in.series,
            entry.in.pages,
            entry.in.address,
            entry.date.year,
            entry.in.publisher,
        ]
    elseif entry.type == "manual"
        temp = [entry.in.organization, entry.in.address, entry.date.year]
    elseif entry.type ∈ ["mastersthesis", "phdthesis"]
        temp = [
            (entry.type == "mastersthesis" ? "Master's" : "PhD") * " thesis",
            entry.in.school,
            entry.in.address,
            entry.date.year,
        ]
    elseif entry.type == "misc"
        temp_last =
            get(entry.fields, "note", "") != "" &&
            entry.access.howpublished != "" &&
            entry.date.year != ""
        temp = [
            entry.access.howpublished
            entry.date.year
            get(entry.fields, "note", "")
        ]
    elseif entry.type == "proceedings"
        temp_last = entry.in.publisher != ""
        temp = [
            (entry.in.volume != "" ? "Volume $(entry.in.volume) of " : "") *
            entry.in.series,
            entry.in.address,
            entry.date.year,
            entry.in.publisher,
        ]
        # TODO: check if this old line here was a hidden bug
        # str *= entry.in.publisher != "" ? ". $(entry.in.address)" : ""
    elseif entry.type == "techreport"
        temp = [
            entry.in.number != "" ? "Technical Report $(entry.in.number)" : "",
            entry.in.institution,
            entry.in.address,
            entry.date.year,
        ]
    elseif entry.type == "unpublished"
        temp = [get(entry.fields, "note", ""), entry.date.year]
    end
    if temp_last
        str *= join(filter!(!isempty, temp), ", ", ". ")
    else
        str *= join(filter!(!isempty, temp), ", ")
    end
    if !isempty(str)
        str *= "."
    end
    return str
end

"""
    xyear(entry)

Format the year of an `Entry` for web export.
"""
xyear(entry) = entry.date.year

"""
    xlink(entry)

Format the download link of an `Entry` for web export.
"""
function xlink(entry)
    if entry.access.doi != ""
        return "https://doi.org/" * entry.access.doi
    elseif entry.access.url != ""
        return entry.access.url
    end
    return ""
end

"""
    xfile(entry)

Format the downloadable path of an `Entry` file for web export.
"""
xfile(entry) = "files/$(entry.id).pdf"

"""
    xcite(entry)

Format the BibTeX cite output of an `Entry` for web export.
"""
xcite(entry) = string(entry)

"""
    xlabels(entry)

Format the labels of an `Entry` for web export.
"""
function xlabels(entry)
    str = get(entry.fields, "swp-labels", "")
    if isempty(str)
        str = get(entry.fields, "labels", "")
    end
    if isempty(str)
        return [entry.type]
    end
    return split(str, r"[\n\r ]*,[\n\r ]*")
end

"""
    Publication

A structure to store all the information necessary to web export.
"""
struct Publication
    id::String
    type::String
    title::String
    names::String
    in::String
    year::String
    link::String
    file::String
    cite::String
    labels::Vector{String}
end

"""
    Publication(entry)

Construct a `Publication` (compatible with web export) from an `Entry`.
"""
function Publication(entry)
    id = entry.id
    type = entry.type
    title = entry.title
    names = xnames(entry)
    in_ = xin(entry)
    year = xyear(entry)
    link = xlink(entry)
    file = xfile(entry)
    cite = export_bibtex(entry)
    labels = xlabels(entry)
    return Publication(id, type, title, names, in_, year, link, file, cite, labels)
end

"""
    export_web(bibliography::DataStructures.OrderedDict{String,BibInternal.Entry})
Export a bibliography in internal format to the web format of the [StaticWebPages.jl](https://github.com/Humans-of-Julia/StaticWebPages.jl) package. Also used by [DocumenterCitations.jl](https://github.com/ali-ramadhan/DocumenterCitations.jl).
"""
function export_web(bibliography)
    entries = Vector{Publication}()
    for entry in values(bibliography)
        p = Publication(entry)
        push!(entries, p)
    end
    return entries
end


text:
# You barely reach the safety of the cave when the whale smashes into the cave
# mouth, collapsing it. Sensors indicate another exit to this cave at a much
# greater depth, so you have no choice but to press on.
#
# As your submarine slowly makes its way through the cave system, you notice
# that the four-digit seven-segment displays in your submarine are
# malfunctioning; they must have been damaged during the escape. You'll be in a
# lot of trouble without them, so you'd better figure out what's wrong.
#
# Each digit of a seven-segment display is rendered by turning on or off any of
# seven segments named a through g:
#
#   0:      1:      2:      3:      4:
#  aaaa    ....    aaaa    aaaa    ....
# b    c  .    c  .    c  .    c  b    c
# b    c  .    c  .    c  .    c  b    c
#  ....    ....    dddd    dddd    dddd
# e    f  .    f  e    .  .    f  .    f
# e    f  .    f  e    .  .    f  .    f
#  gggg    ....    gggg    gggg    ....
#
#   5:      6:      7:      8:      9:
#  aaaa    aaaa    aaaa    aaaa    aaaa
# b    .  b    .  .    c  b    c  b    c
# b    .  b    .  .    c  b    c  b    c
#  dddd    dddd    ....    dddd    dddd
# .    f  e    f  .    f  e    f  .    f
# .    f  e    f  .    f  e    f  .    f
#  gggg    gggg    ....    gggg    gggg
#
# So, to render a 1, only segments c and f would be turned on; the rest would
# be off. To render a 7, only segments a, c, and f would be turned on.
#
# The problem is that the signals which control the segments have been mixed up
# on each display. The submarine is still trying to display numbers by
# producing output on signal wires a through g, but those wires are connected
# to segments randomly. Worse, the wire/segment connections are mixed up
# separately for each four-digit display! (All of the digits within a display
# use the same connections, though.)
#
# So, you might know that only signal wires b and g are turned on, but that
# doesn't mean segments b and g are turned on: the only digit that uses two
# segments is 1, so it must mean segments c and f are meant to be on. With just
# that information, you still can't tell which wire (b/g) goes to which segment
# (c/f). For that, you'll need to collect more information.
#
# For each display, you watch the changing signals for a while, make a note of
# all ten unique signal patterns you see, and then write down a single four
# digit output value (your puzzle input). Using the signal patterns, you should
# be able to work out which pattern corresponds to which digit.
#
# For example, here is what you might see in a single entry in your notes:
#
# acedgfb cdfbe gcdfa fbcad dab cefabd cdfgeb eafb cagedb ab |
# cdfeb fcadb cdfeb cdbaf
#
# (The entry is wrapped here to two lines so it fits; in your notes, it will
# all be on a single line.)
#
# Each entry consists of ten unique signal patterns, a | delimiter, and finally
# the four digit output value. Within an entry, the same wire/segment
# connections are used (but you don't know what the connections actually are).
# The unique signal patterns correspond to the ten different ways the submarine
# tries to render a digit using the current wire/segment connections. Because 7
# is the only digit that uses three segments, dab in the above example means
# that to render a 7, signal lines d, a, and b are on. Because 4 is the only
# digit that uses four segments, eafb means that to render a 4, signal lines e,
# a, f, and b are on.
#
# Using this information, you should be able to work out which combination of
# signal wires corresponds to each of the ten digits. Then, you can decode the
# four digit output value. Unfortunately, in the above example, all of the
# digits in the output value (cdfeb fcadb cdfeb cdbaf) use five segments and
# are more difficult to deduce.
#
# For now, focus on the easy digits. Consider this larger example:
#
# be cfbegad cbdgef fgaecd cgeb fdcge agebfd fecdb fabcd edb |
# fdgacbe cefdb cefbgd gcbe
# edbfga begcd cbg gc gcadebf fbgde acbgfd abcde gfcbed gfec |
# fcgedb cgb dgebacf gc
# fgaebd cg bdaec gdafb agbcfd gdcbef bgcad gfac gcb cdgabef |
# cg cg fdcagb cbg
# fbegcd cbd adcefb dageb afcb bc aefdc ecdab fgdeca fcdbega |
# efabcd cedba gadfec cb
# aecbfdg fbg gf bafeg dbefa fcge gcbea fcaegb dgceab fcbdga |
# gecf egdcabf bgf bfgea
# fgeab ca afcebg bdacfeg cfaedg gcfdb baec bfadeg bafgc acf |
# gebdcfa ecba ca fadegcb
# dbcfg fgd bdegcaf fgec aegbdf ecdfab fbedc dacgb gdcebf gf |
# cefg dcbef fcge gbcadfe
# bdfegc cbegaf gecbf dfcage bdacg ed bedf ced adcbefg gebcd |
# ed bcgafe cdgba cbgef
# egadfb cdbfeg cegd fecab cgb gbdefca cg fgcdab egfdb bfceg |
# gbdfcae bgc cg cgb
# gcafb gcf dcaebfg ecagb gf abcdeg gaef cafbge fdbac fegbdc |
# fgae cfgab fg bagce
#
# Because the digits 1, 4, 7, and 8 each use a unique number of segments, you
# should be able to tell which combinations of signals correspond to those
# digits. Counting only digits in the output values (the part after | on each
# line), in the above example, there are 26 instances of digits that use a
# unique number of segments (highlighted above).
#
# In the output values, how many times do digits 1, 4, 7, or 8 appear?

cnt = sum(eachline("input.txt")) do l
  mapreduce(+, split(split(l, " | ")[2])) do w
    len = length(w)
    len == 2 || len == 4 || len == 3 || len == 7
  end
end

println("Result: ", cnt)


text:
 # at the moment these are just Float64, so most of the precision below is lost
 # would be better to make these "Irrational" but that can wait a little as there
 # are worse pieces of the code
 const stieltjes_n = [γ,
                      -0.0728158454836767248605863758749013191377363383, 	# A082633
                      -0.0096903631928723184845303860352125293590658061, 	# A086279
                      +0.0020538344203033458661600465427533842857158044, 	# A086280
                      +0.0023253700654673000574681701775260680009044694, 	# A086281
                      +0.0007933238173010627017533348774444448307315394, 	# A086282
                      -0.0002387693454301996098724218419080042777837151, 	# A183141
                      -0.0005272895670577510460740975054788582819962534, 	# A183167
                      -0.0003521233538030395096020521650012087417291805, 	# A183206
                      -0.0000343947744180880481779146237982273906207895, 	# A184853
                      +0.0002053328149090647946837222892370653029598537, 	# A184854
                      ]

# we should probably implement some code to actually calculate these at least for checking,
# but crude sum is very very slow, and a lookup table is all we really need
# here are some more values if needed http://fredrikj.net/math/stieltjes100k20d.txt

# calculation
# https://hal.archives-ouvertes.fr/hal-01026562/document
# https://mathoverflow.net/questions/105535/how-to-check-numerical-precision-of-my-computation-of-stieltjes-constants
# https://arxiv.org/pdf/1309.2877.pdf
#
#
# Lambert-W for calculations https://github.com/jlapeyre/LambertW.jl

# old values from Ainsworth and Howell
# gamma = [0.57721566490153
#          -0.07281584548368,
#          -0.00969036319287,
#          0.00205383442030,
#          0.00232537006546,
#          0.00079332381728,
#          -0.0002387693455,
#          -0.0005272895671,
#          -0.0003521233539,
#          -0.0000343947747,
#          ]

"""
    stieltjes(n)

 Provides the first 10 Stieltjes (generalized Euler-Mascheroni) constants (see
 Abramowitz and Stegunm, 23.2.5) or [https://en.wikipedia.org/wiki/Stieltjes_constants](https://en.wikipedia.org/wiki/Stieltjes_constants).

 There is a table at "The Generalized Euler-Mascheroni Constants", O.R. Ainsworth and L.W.Howell
 NASA Technical Paper 2264, Jan 1984
[https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19840007812.pdf](https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19840007812.pdf)
 but the OEIS has more accurate values, which will be useful when I get around to
 bit float versions of the code.

 Note that stieltjes(0) = γ, the Euler–Mascheroni constant, also called just Euler's constant.
 [https://en.wikipedia.org/wiki/Euler-Mascheroni_constant](https://en.wikipedia.org/wiki/Euler%E2%80%93Mascheroni_constant)

## Arguments
* `n::Integer`: the number of elements to compute.

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> stieltjes(0)
0.5772156649015329
```
"""
function stieltjes(n::Integer)
    if n<0
        throw(DomainError(n))
    elseif n>10
        throw(DomainError(n, "Only the first 11 Stieltjes numbers are defined so far."))
    end
    return stieltjes_n[n+1]
end


text:
module Ex5_13
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.13 - Solution of an EMMP using Procedure EMMS
println("Example 5.13")

# define s as an improper transfer function
s = rtf('s');
# enter Gu(s), Gf(s) and Mr(s)
Gu = [s/(s^2+3*s+2) 1/(s+2);
     s/(s+1) 0;
      0 1/(s+2)];
Gf = [s/(s^2+3*s+2) 1/(s+2);
     s/(s+1) 0;
      0 1/(s+2)];
Mr = dss(eye(2));                  # enter Mr(s)
p, mf = size(Gf); mu = size(Gu,2);

# compute left nullspace basis as Q1(s) = [ I -Gu(s) ]
Q1 = dss([I -Gu]); Rf = dss(Gf)

# check solvability condition
rank(evalfr(Rf,rand())) == mf || error("No solution exist")

# check for unstable or infinite zeros
gzero(Rf)   # zeros at infinity and in the origin exist

atol = 1.e-7                # tolerance for rank tests
sdeg = -1                   # set stability degree

# using minimal dynamic covers, compute Q2 such that Q2*Q1 is a lower order
# left annihilator and Q2*Rf full row rank;
# select basis vectors [2 3] of Q1 to combine them with basis vector 1
# to obtain Rf_Q = Q2*[Rf Q1], with Q2*Rf in the leading mf columns and
# Q2*Q1 in the trailing p+mu columns
cinp = [ 2, 3, 1]
Rf_Q = glmcover1([Rf[cinp,:] Q1[cinp,:]], mf; atol)[1]

# compute the irreducible realization of Qtilde = Mr*(inv(Rf)*Q) by
# first solving the linear rational matrix equation Rf*X = Q
X = grsol(Rf_Q, p+mu; atol)[1]
Qtilde = gir(Mr*X; atol)

# compute stable and proper Q = Q4*Qtilde with suitable diagonal Q4 = M
Q = dss(zeros(0,p+mu)); M = dss(zeros(0,0));
for i = 1:mf
    Qi, Mi = glcf(Qtilde[i,:]; atol, sdeg);
    # scale with gain to fit example
    sc = zpk(dss2rm(Mi; atol1 = atol, atol2 = atol)[1,1])[3]
    global Q = [Q; Qi/sc];
    global M = append(M, Mi/sc);
end

# convert to standard state space representation
Q = gss2ss(Q)[1]; M = gss2ss(M)[1]

# display results
println("Q = $(dss2rm(Q, atol = 1.e-7))")
println("M = $(dss2rm(M, atol = 1.e-7))")

# check solution
G = dss([Gu Gf;eye(mu,mu+mf)]); F = [zeros(mf,mu) M*Mr];
@test iszero(Q*G-F; atol)

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
# Autogenerated wrapper script for LibYAML_jll for x86_64-apple-darwin
export libyaml

JLLWrappers.@generate_wrapper_header("LibYAML")
JLLWrappers.@declare_library_product(libyaml, "@rpath/libyaml-0.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libyaml,
        "lib/libyaml-0.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module A1R1S2
export GameEnv, GameSpec
include("game.jl")
end


text:
export AffineChartHomotopy, on_affine_chart

"""
    AffineChartHomotopy(H::AbstractHomotopy, v::PVector{T,N})

Given a homotopy ``H(x,t): (ℙ^{m_1} × ⋯ × ℙ^{m_N}) × ℂ → ℂⁿ`` this creates a new affine
homotopy ``H̄`` which operates on the affine chart defined by the vector
``v ∈ ℙ^{m_1} × ⋯ × ℙ^{m_N}`` and the augmented conditions ``vᵀx = 1``.
"""
struct AffineChartHomotopy{H<:AbstractHomotopy,N} <: AbstractHomotopy
    homotopy::H
    chart::PVector{ComplexF64,N}
end

"""
    on_affine_chart(H::Union{Homotopy,AbstractHomotopy}, proj_dims)

Construct an `AffineChartHomotopy` on a randomly generated chart `v`. Each entry is drawn
idepdently from a univariate normal distribution.
"""
function on_affine_chart(H::Homotopy, proj_dims = nothing)
    on_affine_chart(ModelKitHomotopy(H), proj_dims)
end
function on_affine_chart(H::AbstractHomotopy, proj_dims = nothing)
    if proj_dims === nothing
        proj_dims = [size(H, 2) - 1]
    end
    N = length(proj_dims)
    chart = PVector(randn(ComplexF64, sum(proj_dims) + N), tuple(proj_dims...))
    AffineChartHomotopy(H, chart)
end

function Base.size(H::AffineChartHomotopy{<:AbstractHomotopy,N}) where {N}
    m, n = size(H.homotopy)
    (m + N, n)
end

function set_solution!(x::AbstractVector, H::AffineChartHomotopy, y::AbstractVector, t)
    x .= y
    on_chart!(x, H.chart)
    x
end

start_parameters!(H::AffineChartHomotopy, p) = start_parameters!(H.homotopy, p)
target_parameters!(H::AffineChartHomotopy, p) = target_parameters!(H.homotopy, p)

function on_chart!(x::AbstractVector, v::PVector)
    ranges = dimension_indices(v)
    for range in ranges
        λ = zero(eltype(x))
        @inbounds for i in range
            λ += v[i] * x[i]
        end
        λ⁻¹ = @fastmath inv(λ)
        for i in range
            x[i] *= λ⁻¹
        end
    end
    x
end

function evaluate!(u, H::AffineChartHomotopy{<:Any,N}, x::AbstractVector, t) where {N}
    evaluate!(u, H.homotopy, x, t)
    evaluate_chart!(u, H.chart, x)
    u
end

function evaluate_and_jacobian!(
    u,
    U,
    H::AffineChartHomotopy{<:Any,N},
    x::AbstractVector,
    t,
) where {N}
    evaluate_and_jacobian!(u, U, H.homotopy, x, t)
    evaluate_chart!(u, H.chart, x)
    jacobian_chart!(U, H.chart, x)
    nothing
end

function taylor!(u, v::Val, H::AffineChartHomotopy, tx, t)
    u .= zero(eltype(u))
    taylor!(u, v, H.homotopy, tx, t)
    # affine chart part is always zero since it is an affine linear form
    u
end


text:
module Stack

using GenServers

export start, push, pop, info

# Client (interface)

start(default) = genserver(@__MODULE__, default)
start(name, default) = genserver(@__MODULE__, default, name=name)

push(srv, element) = cast(srv, Val(:push), element)

pop(srv) = call(srv, Val(:pop))

info(srv) = first(query(srv))

reset(srv) = cast(srv, Val(:empty))

stop(srv) = exit!(srv)

# Server (callbacks)

init(stack) = stack

oncast(stack, ::Val{:push}, element) = push!(stack, element)
oncast(stack, ::Val{:empty}) = empty!(stack)

oncall(stack, ::Val{:pop}) = pop!(stack)

end


text:
function sum_border(A::Matrix{Float64})
    sum = 0.
    w, h = size(A)
    for i = 1:w, j = 1:h
        if i == 1 || j == 1 || i == w || j == h
            sum += abs(A[i, j])
        end
    end
    sum
end

function test_equidistance(df::DisplacementField)
    @testset "equidistance" begin
        w, h = size(df.delta_X)
        for i = 1:w, j = 1:h
            @test_approx_eq df.x[j] (j-1)/(h-1)
            @test_approx_eq df.y[i] (i-1)/(w-1)
        end
    end
end

@testset "uniform_displacement" begin
    df = uniform_displacement(5, 8)
    show(df); println()
    showcompact(df); println()
    @test typeof(df) <: DisplacementField
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)

    srand(123)
    df = uniform_displacement(5, 8, .2)
    @test typeof(df) <: DisplacementField
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)
    #@plottest "uniform_displacement" plot(df)

    srand(123)
    df = uniform_displacement(5, 8, scale = .2)
    @test typeof(df) <: DisplacementField
    @test typeof(df) <: ImageTransformation
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)
    #@plottest "uniform_displacement" plot(df)

    srand(123)
    df = uniform_displacement(5, 8, .2, true, true)
    @test typeof(df) <: DisplacementField
    @test length(df.x) == 5
    @test length(df.y) == 8
    @test size(df.delta_X) == (8, 5)
    @test size(df.delta_Y) == (8, 5)
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)
    #@plottest "uniform_displacement" plot(df)

    df = uniform_displacement(10, 9, static_border = false)
    @test typeof(df) <: DisplacementField
    @test sum_border(df.delta_X) > 0.
    @test sum_border(df.delta_Y) > 0.
    test_equidistance(df)

    df = uniform_displacement(10, 9, .1, false, true)
    @test typeof(df) <: DisplacementField
    @test length(df.x) == 10
    @test length(df.y) == 9
    @test size(df.delta_X) == (9, 10)
    @test size(df.delta_Y) == (9, 10)
    @test sum_border(df.delta_X) > 0.
    @test sum_border(df.delta_Y) > 0.
    test_equidistance(df)
end

@testset "gaussian_displacement" begin
    df = gaussian_displacement(10, 12)
    show(df); println()
    showcompact(df); println()
    @test typeof(df) <: DisplacementField
    @test typeof(df) <: ImageTransformation
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)

    srand(123)
    df = gaussian_displacement(10, 12, .3, 2)
    @test typeof(df) <: DisplacementField
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)
    #@plottest "gaussian_displacement" plot(df)

    srand(123)
    df = gaussian_displacement(10, 12, scale = .3, sigma = 2)
    @test typeof(df) <: DisplacementField
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)
    #@plottest "gaussian_displacement" plot(df)

    srand(123)
    df = gaussian_displacement(10, 12, scale = .3, sigma = 1)
    @testset "iterations" begin
        srand(123)
        df2 = gaussian_displacement(10, 12, scale = .3, sigma = 1, iterations = 2)
        for i = 2:11, j = 2:9
            if i == 2 || j == 2 || i == 11 || j == 9
                @test abs(df2.delta_X[i,j]) != abs(df.delta_X[i,j])
                @test abs(df2.delta_Y[i,j]) != abs(df.delta_Y[i,j])
            end
        end
    end

    srand(123)
    df = gaussian_displacement(10, 12, .3, 2, 1, true, true)
    @test typeof(df) <: DisplacementField
    @test length(df.x) == 10
    @test length(df.y) == 12
    @test size(df.delta_X) == (12, 10)
    @test size(df.delta_Y) == (12, 10)
    @test sum_border(df.delta_X) == 0.
    @test sum_border(df.delta_Y) == 0.
    test_equidistance(df)
    #@plottest "gaussian_displacement" plot(df)

    df = gaussian_displacement(10, 12, static_border = false)
    @test typeof(df) <: DisplacementField
    @test sum_border(df.delta_X) > 0.
    @test sum_border(df.delta_Y) > 0.
    test_equidistance(df)

    df = gaussian_displacement(10, 12, .3, 2, 1, false, true)
    @test typeof(df) <: DisplacementField
    @test sum_border(df.delta_X) > 0.
    @test sum_border(df.delta_Y) > 0.
    test_equidistance(df)
end


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ 3d1e8af8-16ef-11ec-18d8-e7fe83a40de4
begin
    using Pkg
    Pkg.activate("..")
    using Javis, Colors, Plots
end

# ╔═╡ 985f6ae4-c559-4756-b9d7-98696c7cbf7f
begin
    function ground(args...)
        background("white")
        sethue("black")
    end
end

# ╔═╡ c25e5278-b7fa-4ec4-abcc-08617fa08872
begin
    mutable struct HarmonicSum
        α::Float64
        values::Vector{Point}
    end

    function HarmonicSum(α::Float64, start_pos::Point; scaling = 20, n_points = 100)
        sequence = map(k -> 1 / k^α, 1:n_points) .* scaling
        return HarmonicSum(α, map(1:n_points) do k
            start_pos + (0, -reduce(+, sequence[1:k]))
        end)
    end
end

# ╔═╡ ef565322-cc4d-4f9d-8558-c722d490c611
function show_sum(S::HarmonicSum, firstframe, lastframe; color = colorant"red")
    points = S.values
    n_frames = lastframe - firstframe + 1
    pointiter = range(1, length(points) - 1, step = (length(points) - 1) ÷ n_frames)
    for idx in 1:n_frames
        Object(
            (firstframe + idx):lastframe,
            @JShape begin
                setline(3)
                sethue(color)
                line(points[pointiter[idx]], points[pointiter[idx + 1]], :stroke)
            end
        )
    end
end

# ╔═╡ 80489189-55c9-46c3-a665-cb58b8f0152d
function dash_sum(S::HarmonicSum, firstframe, lastframe; color = colorant"red")
    points = S.values
    n_frames = lastframe - firstframe + 1
    pointiter = range(1, length(points) - 1, step = (length(points) - 1) ÷ n_frames)
    for idx in 1:n_frames
        Object(
            (firstframe + idx):lastframe,
            @JShape begin
                sethue(color)
                line(
                    points[pointiter[idx + 1]] - (2, 0),
                    points[pointiter[idx + 1]] + (2, 0),
                    :stroke,
                )
            end
        )
    end
end

# ╔═╡ f9c60c0d-9150-42e7-98ed-83d971366292
begin
    function ζ(;
        n_sums = 100,
        height = 500,
        width = 500,
        firstframe = 1,
        lastframe = 100,
        margin = 10,
        llim = -200,
        rlim = 200,
    )

        n_frames = lastframe - firstframe
        colorpalette = range(colorant"red", colorant"blue", length = n_sums)
        Object(
            @JShape begin
                tickline(
                    Point(-200, 200 - margin),
                    Point(200, 200 - margin),
                    startnumber = 0.0,
                    finishnumber = 2.0,
                )
            end
        )
        Object(JLine(Point(0, -height ÷ 2), Point(0, 200 - 20), color = "black"))

        Xpoints = range(llim, rlim, length = n_sums)
        sums = map(enumerate(Xpoints)) do (idx, x)
            s = HarmonicSum(1.0 + x / rlim, Point(x, 200), n_points = 10000)
            dash_sum(s, firstframe, lastframe, color = colorpalette[idx])
            s
        end
    end
    series_video = Video(500, 500)
    Background(1:250, (args...) -> begin
        background("white")
        sethue("black")
    end)

    mypoint = Point(range(-200, 200, length = 100)[50], 200)
    l1 = @JLayer 1:250 begin
        s = HarmonicSum(1.0 + mypoint.x / 200, mypoint, n_points = 10000)
        dash_sum(s, 1, 250)
    end
    act!(l1, Action(50:100, anim_translate(mypoint, Point(200, 0))))
    act!(l1, Action(150:200, anim_translate(Point(200, 0), mypoint)))

    l2 = @JLayer 101:250 begin
        ζ(n_sums = 100, firstframe = 1, lastframe = 150)
    end

    act!(l2, Action(1:50, appear(:fade)))

    render(series_video, pathname = "../output/series_sum.gif")
end

# ╔═╡ Cell order:
# ╟─3d1e8af8-16ef-11ec-18d8-e7fe83a40de4
# ╟─985f6ae4-c559-4756-b9d7-98696c7cbf7f
# ╟─c25e5278-b7fa-4ec4-abcc-08617fa08872
# ╟─ef565322-cc4d-4f9d-8558-c722d490c611
# ╟─80489189-55c9-46c3-a665-cb58b8f0152d
# ╟─f9c60c0d-9150-42e7-98ed-83d971366292


text:
using Convex, SCS, LightGraphs, LinearAlgebra

"""
    find_weights(A)

Takes in Adjacency matrix and gives weights of edges, P
"""
function find_weights(A)
    # the maximal value in sources and destinations is the number of nodes
    n = maximum([size(A,1); size(A,2)])
    idx = findall(x-> (x==0),A)
    # P = Semidefinite(n)
    P = Variable(n,n)

    constraints = [P*ones(n,1) == ones(n,1),P>=0,P==P']

    idx = findall(x -> x == 0, A)
    for ii = 1:length(idx)
        xpos = idx[ii][1]
        ypos = idx[ii][2]
        constraint3 = (P[xpos,ypos]==0)
        push!(constraints,constraint3)
    end
    problem = minimize(opnorm(P - (1/n)*ones(n,n)),constraints)
    solve!(problem, SCSSolver(verbose=false,max_iters=10^4,eps=1e-9))

    Pval = copy(P.value)
    Pval[(Pval) .< 1e-4] .= 0
    Pval = Symmetric(Pval)
    return Pval
end


text:
@deprecate gaussweights quadweights

quadweights(dmeasure::DiscreteWeight, measure::Weight) =
    quadweights(points(dmeasure), weights(dmeasure), measure)

quadweights(grid::AbstractGrid, weights, measure::Weight) =
    quadweights(grid, measure) ./ weights

quadweights(grid::AbstractGrid, ::Ones, measure::Weight) =
    quadweights(grid, measure)

function quadweights(grid::MappedGrid, measure::MappedWeight)
    @assert iscompatible(forward_map(grid), forward_map(measure))
    sgrid = supergrid(grid)
    quadweights(sgrid, supermeasure(measure))
end

quadweights(grid::ProductGrid, measure::ProductWeight) =
    OuterProductArray(map(quadweights, components(grid), components(measure))...)

quadweights(grid::ChebyshevTNodes{T}, ::ChebyshevTWeight{T}) where {T} =
    GridArrays.ChebyshevTWeights{T}(length(grid))

quadweights(grid::ChebyshevUNodes{T}, ::ChebyshevUWeight{T}) where {T} =
    GridArrays.ChebyshevUWeights{T}(length(grid))

quadweights(grid::LegendreNodes{T}, ::LegendreWeight{T}) where {T} =
    weights(GaussLegendre{T}(length(grid)))

function quadweights(grid::LaguerreNodes{T}, measure::LaguerreWeight{T}) where {T}
	@assert grid.α ≈ measure.α
    weights(GaussLaguerre(length(grid),grid.α))
end

quadweights(grid::HermiteNodes{T}, ::HermiteWeight{T}) where {T} =
    weights(GaussHermite{T}(length(grid)))

function quadweights(grid::JacobiNodes{T}, measure::JacobiWeight{T}) where {T}
    @assert (grid.α ≈ measure.α) & (grid.β ≈ measure.β)
	weights(GaussJacobi(length(grid), grid.α, grid.β))
end

quadweights(grid::FourierGrid{T}, measure::FourierWeight{T}) where {T} =
    NormalizedArray{T}(length(grid))

function quadweights(grid::PeriodicEquispacedGrid{T}, measure::FourierWeight{T}) where {T}
    @assert coverdomain(grid)≈support(measure)
	NormalizedArray{T}(length(grid))
end

function quadweights(grid::MidpointEquispacedGrid{T}, measure::FourierWeight{T}) where {T}
    @assert coverdomain(grid)≈support(measure)
	NormalizedArray{T}(length(grid))
end

function quadweights(grid::PeriodicEquispacedGrid{T}, measure::LebesgueMeasure{T}) where {T}
    @assert coverdomain(grid)≈support(measure)
	DomainSets.width(support(measure))*NormalizedArray{T}(length(grid))
end

function quadweights(grid::MidpointEquispacedGrid{T}, measure::LebesgueMeasure{T}) where {T}
    @assert coverdomain(grid)≈support(measure)
	DomainSets.width(support(measure))*NormalizedArray{T}(length(grid))
end

function quadweights(grid::ChebyshevNodes{T}, ::LegendreWeight{T}) where {T}
	x, w = fejer_first_rule(length(grid), T)
	w
end

function quadweights(grid::ChebyshevExtremae{T}, ::LegendreWeight{T}) where {T}
	x, w = clenshaw_curtis(length(grid)-1, T)
	w
end

function quadweights(grid::ChebyshevNodes{T}, measure::LebesgueMeasure{T}) where {T}
    @assert ≈(map(support, (grid,measure))...)
	quadweights(grid, LegendreWeight{T}())
end

function quadweights(grid::ChebyshevExtremae{T}, measure::LebesgueMeasure{T}) where {T}
    @assert ≈(map(support, (grid,measure))...)
	quadweights(grid, LegendreWeight{T}())
end

quadweights(grid::AbstractGrid, measure::Weight) =
	default_quadratureweights(grid, measure)

function default_quadratureweights(grid, measure)
	@debug "No known quadrature normalization available for grid $(typeof(grid)) with measure $(typeof(measure)) \n choosing Riemann sum normalization"
	riemann_normalization(grid, measure)
end

# Fall back in 1d. We don't know any special properties of the grid: we use a
# Riemann sum and multiply by the weight function of the measure on the grid.
function riemann_normalization(grid::AbstractGrid1d, measure)
	a = infimum(support(measure))
	b = supremum(support(measure))
	riemannsum_weights(grid, a, b) .*  measuresampling(grid, measure)
end

function riemannsum_weights(grid::AbstractGrid, a = leftendpoint(coverdomain(grid)), b = rightendpoint(coverdomain(grid)))
	M = length(grid)
	L = b-a
	t = [grid[end] - L; collect(grid); grid[1]+L]
	T = prectype(grid)
	weights = zeros(T, M)
	for m in 1:M
		weights[m] = (t[m+2]-t[m]) / 2
	end
	weights
end

# For any Lebesgue measure this will be the identity operator
measuresampling(grid::AbstractGrid, measure::LebesgueMeasure) =
	Ones{prectype(grid)}(size(grid))

"""
	measuresampling(grid::AbstractGrid, ::Array, measure::Weight)

Return an array with a weight function of a measure evaluated on a grid.
"""
measuresampling(grid::AbstractGrid, measure::Weight) =
    [unsafe_weightfun(measure, xi) for xi in grid]


text:
Base.@kwdef mutable struct LabelIntervalDisplayProperties{T₁ <: NTuple} <: AbstractDisplayProperties
    id::String
    caption::String
    col::ImVec4 = ImVec4(0.9, 0, 0, 0.3)
    textcol::ImVec4 = ImVec4(0.4, 0.4, 0.4, 1)
    createwindow::Bool = true
    layout::RectangularLayout
    padding::T₁ = (0, 0, 0 ,0)
end

function get_id(p::LabelIntervalDisplayProperties)
    p.id
end

function set_id!(p::LabelIntervalDisplayProperties, id::String)
    p.id = id
end

function get_caption(p::LabelIntervalDisplayProperties)
    p.caption
end

function set_caption!(p::LabelIntervalDisplayProperties, caption::String)
    p.caption = caption
end

function get_col(p::LabelIntervalDisplayProperties)
    p.col
end

function set_col!(p::LabelIntervalDisplayProperties, col::ImVec4)
    p.col = col
end

function get_textcol(p::LabelIntervalDisplayProperties)
    p.textcol
end

function set_textcol!(p::LabelIntervalDisplayProperties, textcol::ImVec4)
    p.textcol = textcol
end

function is_new_window(p::LabelIntervalDisplayProperties)
    p.createwindow
end

function get_layout(p::LabelIntervalDisplayProperties)
    p.layout
end

function set_layout!(p::LabelIntervalDisplayProperties, layout::RectangularLayout)
    p.layout = layout
end

function get_padding(p::LabelIntervalDisplayProperties)
    p.padding
end

function set_padding!(p::LabelIntervalDisplayProperties, padding::NTuple)
    p.padding = padding
end


text:
module hinton

export hintondiag

using Tk
using Base.Graphics
using Cairo

function hintondiag(matrix::Matrix{Float64}; default_size = 400)
    # Makes hinton diagram of the given matrix
    #
    # Parameters
    # ----------
    # matrix : the weight matrix to be visualized
    # default_size : default size along max dimension in pixels

    height, width = size(matrix);

    if height < width
        w_frame = default_size;
        h_frame = w_frame * height / width;
    else
        h_frame = default_size;
        w_frame = h_frame * width / height;
    end

    win = Toplevel("weight matrix", int(w_frame), int(h_frame));
    can = Canvas(win);
    pack(can, expand = true, fill = "both");
    draw_diag(can, matrix);

    # For resizing
    can.resize = can->draw_diag(can, matrix);

end

function draw_diag(can::Canvas, matrix::Matrix{Float64})
    # Actual drawing function, separated to support resizing
    #
    # Parameters
    # ----------
    # can : Canvas
    # matrix : the weight matrix to be visualized

    context = getgc(can);

    w_frame = Tk.width(can);
    h_frame = Tk.height(can);

    set_coords(context, 0, 0, w_frame, h_frame, 0, w_frame, 0, h_frame);
    set_source_rgb(context, 0.5, 0.5, 0.5); # Set background to gray
    paint(context);

    height, width = size(matrix);

    w_cell = w_frame / width;
    h_cell = h_frame / height;

    cell_area = w_cell * h_cell * 0.7; # 0.7 to stop merging of cells
    max_weight = 2 ^ (ceil(log(maximum(abs(matrix))) / log(2)));
    area_multiplier = cell_area / max_weight;

    for x in 1:width
        for y in 1:height
            weight = matrix[y, x]
            if weight > 0
                draw_block(context, w_cell * (x - 0.5), h_cell * (y - 0.5), min(1, weight) * area_multiplier, 1);
            elseif weight < 0
                draw_block(context, w_cell * (x - 0.5), h_cell * (y - 0.5), min(1, -weight) * area_multiplier, 0);
            end
        end
    end

    reveal(can);
    Tk.update();
end

function draw_block(ctx::CairoContext, x, y, area, color)
    # Draws block of given size on the context
    #
    # Parameters
    # ----------
    # ctx : the context of canvas
    # x, y : the center of the block
    # area : area of the square
    # color : 1 = white, 0 = black

    side = sqrt(area);
    rectangle(ctx, x - side / 2, y - side / 2, side, side);

    if color == 1
        set_source_rgb(ctx, 1, 1, 1);
    else
        set_source_rgb(ctx, 0, 0, 0);
    end
    fill(ctx);
end

end # module


text:
using ProjectMoana
using ColorTypes

@testset "Colors: Basic operations" begin
    @test RGB(.1,.3,.5)+RGB(.4,.6,.2) ≈ RGB(.5,.9,.7)
    @test RGB(.8,.7,.5)-RGB(.4,.6,.2) ≈ RGB(.4,.1,.3)
    @test RGB(1,.4,.4)*2 ≈ RGB(2,.8,.8)
    @test 3*RGB(.1,.2,.3) ≈ RGB(.3,.6,.9)
end

col1 = RGB(1.0, 2.0, 3.0)
col2 = RGB(9.0, 5.0, 7.0)

@testset "Colors: Luminosity function" begin
    @test isapprox(luminosity(col1), 2.0)
    @test isapprox(luminosity(col2), 7.0)
end


text:
using IntervalUnionArithmetic
using Test

@testset "constructor" begin

    x = 1..2
    a = intervalUnion(x)
    @test length(a.v) == 1
    @test a.v[1] == x

    a2 = intervalUnion(1, 2)
    @test length(a.v) == 1
    @test a.v[1] == x

    x2 = 3 .. 4
    b = intervalUnion([x, x2])
    @test length(b.v) == 2
    @test b.v == [x, x2]

    b2 = x ∪ x2
    @test b2.v == b.v

    x3 = interval(5, 6)
    b3 = b2 ∪ x3
    @test b3.v == [x, x2, x3]

    # Remove empties
    c = x ∪ ∅
    @test length(c.v) == 1
    @test c.v[1] == x

    # Remove empties
    c2 = x ∪ ∅ ∪ ∅ ∪ ∅ ∪ ∅ ∪ ∅
    @test length(c2.v) == 1
    @test c2.v[1] == x

    # Don't remove if IntervalUnion is called
    c3 = IntervalUnion([x, ∅])
    @test length(c3.v) == 2
    @test c3.v == [x, ∅]

    # Condense
    d = interval(1, 2) ∪ interval(1.5, 3)
    @test length(d.v) == 1
    @test d.v[1] == interval(1, 3)

    d2 = interval(1, 2) ∪ interval(2, 3)
    @test length(d2.v) == 1
    @test d2.v[1] == interval(1, 3)

end

@testset "close gaps" begin

end


text:
# sqrt example

using Distributions
using KernelDensityEstimate, KernelDensityEstimatePlotting
using IncrementalInference
using Gadfly, DataFrames

# requried for overloading with new factors
import IncrementalInference: getSample

include(joinpath(@__DIR__, "SquareRootTypes.jl"))

## Direct example

fg = initfg()

addVariable!(fg, :x, ContinuousScalar)
# addVariable!(fg, :y, x0,N=N)

# TODO perhaps make Mixture instead for multiple computation
pts = rand(Distributions.Normal(4.0,0.05),100)   #;rand(Distributions.Normal(144.0,0.05),N)]
md = kde!(pts)

npx = NumbersPrior(md)
addVariable!(fg, :xy, ContinuousScalar)

addFactor!(fg, [:xy;], npx)

xty = Square(Distributions.Normal(0.0,0.01))
addFactor!(fg, [:x,:xy], xty)

# drawGraph(fg)

# initialize from prior
doautoinit!(fg, :xy)
# initialize any random numbers for the square root initial value
manualinit!(fg, :x, randn(1,100))

# find solution
tree, smt, hist = solveTree!(fg)

## plot the result
plotKDE(map(x->getKDE(fg,x), [:x; :xy]))


text:
using Test
using GeophysicalModelGenerator

function main_gravity()
############# Input #############
# survey
x        = [-20.,   100.]
y        = [100., 200.]
z        = [-100., 10.]
nx       = 111
ny       = 101
nz       = 91

# spheres
centerX  = [  50.,   10.,   80.]
centerY  = [ 150.,  150.,  180.]
centerZ  = [ -30.,  -10.,  -50.]
radius   = [  10.,    3.,   10.]
rho      = [2650., 2600., 2750.]

# background density model
background = [0. 2700.]
#################################

######### Set things up #########
# constants
G        = 6.67408e-11
mGal     = 1e5

# coordinates
x_vec    = LinRange(x[1],x[2],nx)
y_vec    = LinRange(y[1],y[2],ny)
z_vec    = LinRange(z[1],z[2],nz)
Y,X,Z    = meshgrid(y_vec,x_vec,z_vec)

# reference model
RefMod   = zeros(nz)
RefMod[findall(x->x>=0, z_vec)] .= background[1]
RefMod[findall(x->x<0,  z_vec)] .= background[2]

# densities
RHO      = zeros(nx,ny,nz)
RHO[findall(x->x>=0, Z)] .= background[1]
RHO[findall(x->x<0,  Z)] .= background[2]

# put in spheres
if !(length(centerX) == length(centerY) && length(centerX) == length(centerZ) && length(centerX) == length(radius))
    error("Sphere definition is wrong")
end

numSpheres = length(rho)

for i = 1 : numSpheres
    d = (X .- centerX[i]).^2 .+ (Y .- centerY[i]).^2 .+ (Z .- centerZ[i]).^2
    RHO[findall(x->x<radius[i]^2, d)] .= rho[i]
end
#################################

############ Compute ############
# Analytical Solution
ana = zeros(nx,ny)
for iS = 1 : numSpheres
    for iX = 1 : nx
        for iY = 1 : ny
            d          = ((x_vec[iX] - centerX[iS])^2 + (y_vec[iY] - centerY[iS])^2 + (centerZ[iS])^2)^0.5
            depth      = -centerZ[iS]
            ana[iX,iY] = ana[iX,iY] + (4 * π * G * (radius[iS]^3) * (rho[iS] - background[2]) * depth) / (3 * (d^3)) * mGal
        end
    end
end

# Voxel Code
dg1, gradX1, gradY1 = voxGrav(X, Y, Z, RHO, refMod = RefMod, outName = "Benchmark1", printing = false)

# Other options (same result)
dg2, gradX2, gradY2 = voxGrav(X, Y, Z, RHO, rhoTol = 1, refMod = "SW", outName = "Benchmark2", printing = false)

# Other options (different result)
dg3, gradX3, gradY3 = voxGrav(X, Y, Z, RHO, rhoTol = 70, refMod = "NW", outName = "Benchmark3", printing = false)
#################################

############# Test ##############
maxVal = maximum(broadcast(abs, ana))
maxErr = maxVal / 20

# check differences
diff1 = broadcast(abs, ana-dg1)
diff2 = broadcast(abs, ana-dg2)
diff3 = broadcast(abs, ana-dg3)

# check
@test maximum(diff1) < maxErr
@test maximum(diff2) < maxErr
@test maximum(diff3) > maxErr
#################################
end

main_gravity()


text:
using Pkg

dependencies = [
    "Statistics",
    "FFTW",
    "LinearAlgebra"
    "Statistics",
    "SparseArrays",
    "BSON",
    "Plots",
    "ArgParse",
    "DelimitedFiles",
    "Interpolations"
]

Pkg.add(dependencies)

max_stars_repo_name: jquetzalcoatl/DynamicStructureFactor

max_issues_repo_name: jquetzalcoatl/DynamicStructureFactor

max_forks_repo_name: jquetzalcoatl/DynamicStructureFactor


text:
module SXTMerger
using DataFrames, FITSIO, FITSIO.Libcfitsio

export sxt_l2evtlist_merge

include("sxt_l2evtlist_merge.jl")
include("sxt_df2vec.jl")
include("sxt_gti_merge.jl")
include("sxt_badpix_merge.jl")
include("sxt_evt2df.jl")

end

max_stars_repo_name: gulabd/SXTMerger

max_issues_repo_name: gulabd/SXTMerger


text:
function correct_grid_size(FT)
    grid = RegularCartesianGrid(FT; size=(4, 6, 8), length=(2π, 4π, 9π))

    # Checking ≈ as the grid could be storing Float32 values.
    return (grid.Nx ≈ 4  && grid.Ny ≈ 6  && grid.Nz ≈ 8 &&
            grid.Lx ≈ 2π && grid.Ly ≈ 4π && grid.Lz ≈ 9π)
end

function correct_halo_size(FT)
    grid = RegularCartesianGrid(FT; size=(4, 6, 8), length=(2π, 4π, 9π), halo=(1, 2, 3))
    return (grid.Hx == 1  && grid.Hy == 2  && grid.Hz == 3)
end

function faces_start_at_zero(FT)
    grid = RegularCartesianGrid(FT; size=(10, 10, 10), length=(2π, 2π, 2π))
    return grid.xF[1] == 0 && grid.yF[1] == 0 && grid.zF[end] == 0
end

function end_faces_match_grid_length(FT)
    grid = RegularCartesianGrid(FT; size=(12, 13, 14), length=(π, π^2, π^3))
    return (grid.xF[end] - grid.xF[1] ≈ π   &&
            grid.yF[end] - grid.yF[1] ≈ π^2 &&
            grid.zF[end] - grid.zF[1] ≈ π^3)
end

function ranges_have_correct_length(FT)
    Nx, Ny, Nz = 8, 9, 10
    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1, 1, 1))
    return (length(grid.xC) == Nx && length(grid.xF) == Nx+1 &&
            length(grid.yC) == Ny && length(grid.yF) == Ny+1 &&
            length(grid.zC) == Nz && length(grid.zF) == Nz+1)
end

# See: https://github.com/climate-machine/Oceananigans.jl/issues/480
function no_roundoff_error_in_ranges(FT)
    Nx, Ny, Nz = 1, 1, 64
    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1, 1, π/2))
    return length(grid.zC) == Nz
end

@testset "Grids" begin
    println("Testing grids...")

    @testset "Grid initialization" begin
        println("  Testing grid initialization...")
        for FT in float_types
            @test correct_grid_size(FT)
            @test correct_halo_size(FT)
            @test faces_start_at_zero(FT)
            @test end_faces_match_grid_length(FT)
            @test ranges_have_correct_length(FT)
            @test no_roundoff_error_in_ranges(FT)
        end
    end

    @testset "Grid dimensions" begin
        println("  Testing grid constructor errors...")
        L = (100, 100, 100)
        for FT in float_types
            @test isbitstype(typeof(RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 1, 1))))

            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32,), length=L)
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 64), length=L)
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32, 16), length=L)

            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32.0), length=(1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(20.1, 32, 32), length=(1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, nothing, 32), length=(1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, "32", 32), length=(1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, nothing, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, "1", 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, 1, 1), halo=(1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, 1, 1), halo=(1.0, 1, 1))

            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=2)
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), y=[1, 2])
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), z=(-π, π))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=1, y=2, z=3)
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(0, 2), z=4)
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=("0", "1"))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=(1, 2, 3))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(1, 0), y=(1//7, 5//7), z=(1, 2))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 2, 3), x=(0, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 2, 3), x=(0, 1), y=(1, 5), z=(-π, π))
        end
    end
end


text:
using PageDown

include("page.jl")
include("macro.jl")


text:
using ITensors,
      ITensorLocalMPO

import ITensors: op

function op(::OpName"Givens", ::SiteType"Electron", s1::Index, s2::Index; θ, θup = θ, θdn = θ)
  ampo = AutoMPO()
  ampo += (θup, "Cdagup", 1, "Cup", 2)
  ampo -= (θup, "Cdagup", 2, "Cup", 1)
  ampo += (θdn, "Cdagdn", 1, "Cdn", 2)
  ampo -= (θdn, "Cdagdn", 2, "Cdn", 1)
  return exp(prod(MPO(ampo, [s1, s2])))
end

function main(N; θ, pattern)
  s = siteinds("Electron", N; conserve_qns = true)

  U = 1.0
  a = AutoMPO()
  for n in 1:N-1
    a .+= -1.0, "Cdagup", n, "Cup", n+1
    a .+= -1.0, "Cdagup", n+1, "Cup", n
    a .+= -1.0, "Cdagdn", n, "Cdn", n+1
    a .+= -1.0, "Cdagdn", n+1, "Cdn", n
  end
  for n in 1:N
    a .+= U, "Nupdn", n
  end
  H = MPO(a, s)

  if pattern == "staircase"
    gates = [("Givens", (n, n+1), (θ = θ,)) for n in 1:N-1]
  elseif pattern == "brick"
    gates = [("Givens", (n, n+1), (θ = θ,)) for n in 1:2:N-1]
    append!(gates, [("Givens", (n, n+1), (θ = θ,)) for n in 2:2:N-1])
  else
    error("Gate pattern $pattern not supported")
  end
  U = ops(s, gates)

  # Without basis change
  UH1 = H
  mpo_evolution_succeeded = false
  UH1 = try
    UH1 = apply(U, UH1; apply_dag = true, cutoff = 1e-8)
    UH1 = apply(U, UH1; apply_dag = true, cutoff = 1e-8)
    mpo_evolution_succeeded = true
    UH1
  catch
    println("apply(U, H::MPO, apply_dag = true, [...]) failed, likely because of overflow")
  end

  # With basis change
  basis_transformation = scale_bases(H; normalize_op = "Id")

  UX, HX = combine_and_transform(U, H, basis_transformation...)
  UHX = HX
  UHX = apply(UX, UHX; cutoff = 1e-8)
  UHX = apply(UX, UHX; cutoff = 1e-8)

  # Fix normalization
  normsUHX = norm.(UHX)
  UHX ./= normsUHX
  lognormsUHX = log.(normsUHX)
  lognormUHX = sum(lognormsUHX)
  UHX .*= exp(lognormUHX / N)

  UH2 = inv_transform_and_uncombine(UHX, basis_transformation...)

  @show maxlinkdim(H), maxlinkdim(HX), maxlinkdim(UHX), maxlinkdim(UH2)
  mpo_evolution_succeeded && @show maxlinkdim(UH1)
  @show norm(H), norm(HX), norm(UHX), norm(UH2)
  mpo_evolution_succeeded && @show norm(UH1)
  @show maximum(norm, H), maximum(norm, HX), maximum(norm, UHX), maximum(norm, UH2)
  mpo_evolution_succeeded && @show maximum(norm, UH1)
  mpo_evolution_succeeded && @show inner(UH1, UH2) / (norm(UH1) * norm(UH2))

  return nothing
end

max_stars_repo_name: mtfishman/ITensorLocalMPO

max_issues_repo_name: mtfishman/ITensorLocalMPO

max_forks_repo_name: mtfishman/ITensorLocalMPO


text:
export RDFT

"""
`RDFT([domainType=Float64::Type,] dim_in::Tuple [,dims=1])`

`RDFT(dim_in...)`

`RDFT(x::AbstractArray [,dims=1])`

Creates a `LinearOperator` which, when multiplied with a real array `x`, returns the DFT over the dimension `dims`, exploiting Hermitian symmetry.

```julia
julia> RDFT(Float64,(10,10))
ℱ  ℝ^(10, 10) -> ℂ^(6, 10)

julia> RDFT((10,10,10),2)
ℱ  ℝ^(10, 10, 10) -> ℂ^(10, 6, 10)

```

"""
struct RDFT{T <:Number,
	       N,
	       T1<:AbstractFFTs.Plan,
	       T2<:AbstractFFTs.Plan,
	       T3<:AbstractArray{Complex{T},N}
	      } <: LinearOperator
	dim_in::NTuple{N,Int}
	dim_out::NTuple{N,Int}
	A::T1
	At::T2
	b2::T3
	y2::T3
	Zp::ZeroPad{Complex{T},N}
end

# Constructors
#standard constructor

function RDFT(x::AbstractArray{T,N}, dims::Int=1) where {T<:Real,N}
	A = plan_rfft(x,dims)
	b2 = zeros(complex(T),size(x))
	y2 = zeros(complex(T),size(x))
	At = plan_bfft(y2,dims)
	dim_in = size(x)
	dim_out = ()
	for i = 1:N
		dim_out = i == dims ? (dim_out...,div(dim_in[i],2)+1) : (dim_out...,dim_in[i])
	end
	Z = ZeroPad(Complex{T},dim_out,size(b2))
	RDFT{T,N,typeof(A),typeof(At),typeof(b2)}(dim_in,dim_out,A,At,b2,y2,Z)
end

RDFT(T::Type,dim_in::NTuple{N,Int},dims::Int=1) where {N} = RDFT(zeros(T,dim_in),dims)
RDFT(dim_in::NTuple{N,Int},dims::Int=1) where {N} = RDFT(zeros(dim_in),dims)
RDFT(dim_in::Vararg{Int}) = RDFT(dim_in)
RDFT(T::Type,dim_in::Vararg{Int}) = RDFT(T,dim_in)

# Mappings

function mul!(y::T3,
              L::RDFT{T,N,T1,T2,T3},
              b::T4) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}}
	mul!(y,L.A,b)
end

function mul!(y::T4,
              L::AdjointOperator{RDFT{T,N,T1,T2,T3}},
              b::T3) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}}
    A = L.A
	mul!(A.b2,A.Zp,b)
	mul!(A.y2,A.At,A.b2)
	y .= real.(A.y2)
	return y

end

# Properties

size(L::RDFT) = (L.dim_out,L.dim_in)

fun_name(A::RDFT) = "ℱ"

domainType(L::RDFT{T,N}) where {T,N} = T
codomainType(L::RDFT{T,N}) where {T,N} = Complex{T}

is_AAc_diagonal(    L::RDFT) = false #TODO but might be true?
is_invertible(      L::RDFT) = true
is_full_row_rank(   L::RDFT) = true


text:
"""
    StateVecs

Provides a state vector given a `NamedTuple` containing
tuples of the variable name and the number of its components.
"""
module StateVecs

using ..Grids
export StateVec, over_sub_domains, Cut

struct FieldsPerElement{T}
  vals::Vector{T}
end

"""
    StateVec{T, I, NT}

A state vector containing the number of subdomains,
`n_subdomains`, a `NamedTuple` variable mapper, a tuple
of the variable names, and a vector of vectors, containing
the values for all of the variables.
"""
struct StateVec{T, I, NT}
  """
  Number of subdomains.
  """
  n_subdomains::I

  """
  A `NamedTuple` that maps the variable name to its unique index.
  """
  var_mapper::NT

  """
  A tuple of all variable names.
  """
  var_names::Tuple{Vararg{Symbol}}

  """
  A vector of vectors containing all values for all variables.
  """
  fields::Vector{FieldsPerElement{T}}
end

"""
    get_var_mapper(vars)

Get a `NamedTuple` that maps a variable name
and subdomain ID to a unique index. Example:
```
vars = ((:ρ_0, 1), (:w, 3), (:a, 3), (:α_0, 1))
var_mapper = get_var_mapper(vars)
((:ρ_0, :w, :a, :α_0), (ρ_0 = [1], w = [2, 3, 4], a = [5, 6, 7], α_0 = [8]))
```
"""
function get_var_mapper(vars)
  var_names = tuple([v for (v, nsd) in vars]...)
  end_index = cumsum([nsd for (v, nsd) in vars])
  start_index = [1,[1+x for x in end_index][1:end-1]...]
  vals = [collect(a:b) for (a,b) in zip(start_index, end_index)]
  var_mapper = NamedTuple{var_names}(vals)
  return var_names, var_mapper
end

"""
    StateVec(vars::Tuple{Vararg{Tuple{Symbol, I}}}, grid::Grid{T}) where {I, T}

Return a state vector, given a tuple of tuples of variable
names and the number of their subdomains.
"""
function StateVec(vars::Tuple{Vararg{Tuple{Symbol, I}}}, grid::Grid{T}) where {I, T}
  n_subdomains = max([nsd for (v, nsd) in vars]...)
  n_vars = sum([nsd for (v, nsd) in vars])
  var_names, var_mapper = get_var_mapper(vars)
  all_vars = [eltype(grid.dz)(0) for v in 1:n_vars]
  fields = [FieldsPerElement(deepcopy(all_vars)) for k in over_elems(grid)]
  return StateVec(n_subdomains, var_mapper, var_names, fields)
end

function Base.show(io::IO, sv::StateVec)
  println(io, "----------------------- State Vector")
  println(io, "n_subdomains    = ", sv.n_subdomains)
  println(io, "var_mapper      = ", sv.var_mapper)
  headers = [ length(over_sub_domains(sv, name))==1 ? string(name) : string(name)*"_"*string(i_sd)
             for name in sv.var_names for i_sd in over_sub_domains(sv, name)]
  n_vars = length(headers)
  n_elem = length(sv.fields)
  data = reshape([sv[name, k, i_sd] for name in sv.var_names for
    i_sd in over_sub_domains(sv, name) for k in 1:n_elem], n_elem, n_vars)
  Base.print_matrix(io, [reshape(Text.(headers), 1, n_vars); data])
  println(io, "-----------------------")
end

"""
    over_sub_domains(state_vec::StateVec)

Get list of indexes from 1 to the maximum subdomain size.
"""
over_sub_domains(state_vec::StateVec) = 1:state_vec.n_subdomains

"""
    over_sub_domains(state_vec::StateVec, j::Int)

Get list of indexes from 1 to the maximum subdomain size, except the given index.
"""
over_sub_domains(state_vec::StateVec, j::Int) = [i for i in 1:state_vec.n_subdomains if !(i==j)]

"""
    over_sub_domains(state_vec::StateVec, name::Symbol)

Get list of indexes over all subdomains for variable `name`.
"""
over_sub_domains(state_vec::StateVec, name::Symbol) = 1:length(state_vec.var_mapper[name])

Base.getindex(sv::StateVec, name::Symbol, k, i_sd=1) = sv.fields[k].vals[sv.var_mapper[name][i_sd]]
function Base.setindex!(sv::StateVec, val, name::Symbol, k, i_sd = 1)
  sv.fields[k].vals[sv.var_mapper[name][i_sd]] = val
end

abstract type AbstractCut{I} end

"""
    Cut{I} <: AbstractCut{I}

A Cut struct used to slice the state
vector along the grid-element dimension.
This is used to as an API to pass Cuts
into local derivative/interpolation routines.
"""
struct Cut{I} <: AbstractCut{I}
  e::I
end
Base.getindex(sv::StateVec, name::Symbol, Cut::Cut, i_sd=1) = [sv[name, k, i_sd] for k in Cut.e-1:Cut.e+1]

Base.isnan(sv::StateVec) = any([any(isnan.(fpe.vals)) for fpe in sv.fields])
Base.isinf(sv::StateVec) = any([any(isinf.(fpe.vals)) for fpe in sv.fields])

end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
module XLA

@nospecialize

using IRTools, IRTools.All, PyCall, Mjolnir, MacroTools
using IRTools: block
using IRTools.Inner: entry
using MacroTools: @capture
import Mjolnir: AType, Partial, Shape, Multi, Basic, Const, KwFunc, abstract, instead,
  ptuple, widen, @abstract

export @code_xla, @code_hlo, xla, isxla

function __init__()
  global xlaclient = pyimport("jaxlib.xla_client")
  PyCall.npyinitialize()
  requires()
end

include("ir/reloop.jl")
include("ir/builder.jl")
include("ir/ops.jl")

include("compile/passes.jl")
include("compile/convert.jl")
include("compile/rt.jl")

include("lib.jl")

macro code_typed(ex)
  @capture(ex, f_(args__)) || error("@code_typed f(args...)")
  quote
    trace(Const($(esc(f))), xtypeof.(($(esc.(args)...),))...) |> renumber
  end
end

macro code_xla(ex)
  @capture(ex, f_(args__)) || error("@code_xla f(args...)")
  quote
    tr = trace(Const($(esc(f))), xtypeof.(($(esc.(args)...),))...)
    deletearg!(tr, 1)
    convert_xla!(tr, xtypeof(($(esc.(args)...),))) |> renumber
  end
end

macro code_hlo(ex)
  @capture(ex, f_(args__)) || error("@code_hlo f(args...)")
  quote
    tr = trace(Const($(esc(f))), xtypeof.(($(esc.(args)...),))...)
    deletearg!(tr, 1)
    ir = convert_xla!(tr, xtypeof(($(esc.(args)...),)))
    build(controlflow(ir)).as_hlo_text() |> print
  end
end

end # module


text:
g = FactorGraph(3, 4)
@test nv(g) == 7
@test ne(g) == 0
@test_throws ArgumentError add_edge!(g, 1, 2)
@test_throws ArgumentError add_edge!(g, 6, 7)
add_edge!(g, 1, 4)
@test ne(g) == 1

g = FactorGraph(3, 4, 8)
@test g.nvars == 3
@test g.nfacts == 4
@test nv(g) == 7
@test ne(g) == 8

g = BeliefPropagation.random_bipartite_regular_graph(5, 6, 2)
degs = degree(g)
@test all(degs[1:5] .== 2)

g = BeliefPropagation.random_bipartite_regular_graph(5, 6, 2, first_regular=false)
degs = degree(g)
@test all(degs[6:11] .== 2)