text:
export detf, quad, laguerre, estimate_mult, convergence, starting_points!, bounds
export iteration_one, merge12, merge2, zrst, correct

"""
    detf(t, s, x)

Determinant quotient and derivatives.
"""
function detf(t::AbstractMatrix, s::AbstractMatrix, x::T, y::T=T(0)) where T<:Real
    n = size(t, 2)
    S = real(promote_type(eltype(t), eltype(s), T))
    ϵ = eps(S)^2
    ξ = oneunit(S)
    η = zero(S)
    ζ = zero(S)
    κ = 0
    if n == 0
        return ξ, η, ζ, κ
    end
    ξ_1, η_1, ζ_1 = ξ, η, ζ
    tii, sii = real(t[1,1]), real(s[1,1])
    ξ = (tii - sii * x ) - sii * y
    if n == 2000
        println("1 $ξ")
    end
    if iszero(ξ)
        ξ = max(abs(tii), abs(sii * (x + y)), ξ_1) * ϵ
    end
    η = sii / ξ
    ζ = zero(S)
    κ += ξ < 0
    for i = 2:n
        η_0, ζ_0 = η_1, ζ_1
        ξ_1, η_1, ζ_1 = ξ, η, ζ
        tii, tim, sii, sim = ((real(t[1,i]), t[2,i-1], real(s[1,i]), s[2,i-1]))
        f, fp, fpp = quad(tim, sim, x, y)
        tsi = (tii - sii * x) - sii * y
        ξ = tsi - f / ξ_1
        if n == 2000
            println("$i $ξ")
        end
        if iszero(ξ)
            ξ = (abs(tim) + abs(sim) * abs(x + y))^2 * ϵ / ξ_1
        end
        η = (tsi * η_1 + sii - (f * η_0 - fp ) / ξ_1) / ξ
        ζ = (tsi * ζ_1 + 2 * sii * η_1 - (ζ_0 * f - 2 * η_0 * fp + fpp) / ξ_1) / ξ
        κ += (ξ < 0)
    end
    return ξ, η, ζ, κ
end

# real(a * b')
realprod(a::T, b::T) where T<:Real = a * b
realprod(a::T, b::T) where T<:Complex = real(a) * real(b) + imag(a) * imag(b)

"""
    quad(a, b, x::Real)

`f(x) = abs2(a - b * x)` and `f'(x), f''(x)`
"""
function quad(a::T, b::T, x::R, y::R=0) where {T<:Union{Real,Complex},R<:Real}
    f = abs2((widen(a) - b * x) - b * y)
    bb = abs2(b) * 2
    fp = bb * x - realprod(a, b) * 2 + bb * y
    oftype(a, f), fp, bb
end

"""
    laguerre(t, s, x, r)

Calculate one iteration step of Laguerre's iteration method for polynomial roots.
The polynomial is given by `x -> det(t - s * x)`.
"""
function laguerre(t::AbstractMatrix, s::AbstractMatrix, x::Real, r::Integer)
    n = size(t, 2)
    _, η, ζ, κ = detf(t, s, x)
    laguerre(η, ζ, n, r), κ
end
# in case of negative discriminant fall back to Newton step (adapted for multiplicity r)
function laguerre(η::T, ζ::T, n::Integer, r::Integer) where T<:Real
    disc = (η^2 * (n-1) - ζ * n) * (n - r) / r
    disc >= 0 ? n / ( η + copysign(sqrt(disc), η) ) : r / η
end

"""
    estmlt

Estimate multiplicity of eigenvalue number j close to x
"""
function estimate_mult(x, sign::Integer, λ::Vector, j::Integer)
    n = size(λ, 1)
    1 <= j <= n || throw(ArgumentError("j = $j not in range 1:$n"))
    m = j
    lj = λ[j]
    tol = abs(x - lj) * 0.01
    mlt = 1
    while true
        m += sign
        1 <= m <= n || break
        abs(λ[m] - lj) >= tol && break
        mlt += 1
    end
    mlt
end

function convergence(x2, x1, x0)
    ϵ = eps(x2)
    dx2 = abs(x2 - x1)
    #dx2 > ϵ * 1000 && return false
    dx2 <= ϵ && return true
    dx1 = abs(x1 - x0)
    #dx2 >= dx1 && return true
    dx2 < dx1 && dx2^2 / (dx1 - dx2) <= ϵ && return true
    return false
end

κ(x, t, s) = detf(t, s, x, zero(x))[4]

"""
    starting_points!(estev, t, s, a, b)

Select starting points for the eigenvalues of pencil(t, s) in interval `[a,b)`.
`estev` are the eigenvalues of the split system in the same interval, which have been
calculated accurately.
"""
function starting_points!(λ0::Vector{<:Real}, t, s, a::Real, b::Real)
    @assert issorted(λ0)
    m = length(λ0)  # number of ev in input data
    @assert  m == 0 || prevfloat(a) <= λ0[1] <= λ0[end] < nextfloat(b)

    ka = κ(a, t, s)
    kb = κ(b, t, s)
    offseta = false
    offsetb = false
    if m == 0
        if ka < kb
            push!(λ0, a)
            push!(λ0, b)
            kb = ka + 1
            m += 2
            offseta = true
        end
    else
        if m >= 1
            c = λ0[1]
            kc = κ(c, t, s)
            if ka < kc
                pushfirst!(λ0, a)
                m += 1
            end
            offseta = kc <= ka + 1
        end
        if m >= 2
            d = λ0[end]
            kd = κ(d, t, s)
            if kb > kd
                push!(λ0, b)
                m += 1
            end
            offsetb = m - offseta - kb + ka
        end
    end
    #@assert length(λ0) - offseta - offsetb == kb - ka "$m $(length(λ0)) $offseta $offsetb $kb $ka \n $(λ0)"
    λ0, ka, offseta, offsetb
end

"""
    interval(λ, j, i)

The interval for eigenvalue j in interval. Starting-point is the midpoint for inner
indices and one of the boundaries for first and last index.
"""
function interval(λ::Vector, j::Integer)
    n = size(λ, 1)
    λ[max(j-1, 1)], λ[min(j+1, n)]
end

"""
    iteration_one(a, x, b, j, n, r, t, s)

Iterate in interval `[a,b]` for the eigenvalue number `j` of `n`.
The eigenvalue number `j` with multiplicity `r` must be contained in [a, b].
Algorithm: First bisection loops until indication for Laguerre loops.
Test for convergence during Laguerre loop.
"""
function iteration_one(a::T, x::T, b::T, j::Integer, n::Integer, r::Integer, t, s) where T
    @assert a <= x <= b
    @assert 1 <= j <= n
    @assert 1 <= r
    mult = r
    a0, x0, b0 = a, x, b

    xx = ( a0 + b0 ) / 2
    bb = ( b0 - xx ) / 2
    aa = ( a0 - xx ) / 2
    x =  x - xx
    bcount, lcount = 0, 0
    _, η, ζ, κ = detf(t, s, xx, x)
    sig = κ < j
    x1 = NaN
    while a < b
        κ1 = κ
        x0 = x1
        x1 = x
        dx = laguerre(η, ζ, n, mult)

        if sig
            bx = ( bb - x ) / 2
            if 0 < dx <= bx
                x += dx
                lcount += 1
            else
                x += bx
                bcount += 1
            end
            a = x1
        else
            ax = ( aa - x ) / 2
            if ax <= dx < 0
                x += dx
                lcount += 1
            else
                x += ax
                bcount += 1
            end
            b = x1
        end

        convergence(xx+x, xx+x1, xx+x0) && break

        if (lcount == 10 || bcount == 10 ) && n == 1000
            println("bcount=$bcount lcount=$lcount ($a0, $x0, $b0, $j, $n, $r)")
        end
        if lcount > 50
            println("no convergence for $j/$n: $a <= $x <= $b ")
            break
        end

        while x != x1
            _, η, ζ, κ = detf(t, s, xx, x)
            dκ = abs(κ - κ1)
            ( mult <= 1 || dκ <= 1 ) && break
            mult = dκ
            x += (x1 - x) / 2
        end
        sig = κ < j
    end
    xx + x, bcount, lcount
end

function pushres!(res, ev, a, b)
    if a <= ev < b
        push!(res, ev)
    end
    res
end

function merge12(t::AbstractMatrix{T}, s::AbstractMatrix{T}, a, b) where T
    n = size(t, 2)
    res = real(T)[]
    if n == 1
        ev = real(t[1,1]) / real(s[1,1])
        return pushres!(res, ev, a, b), [], []
    elseif n == 2
        s11, s12, s22 = s[1,1], s[2,1], s[1,2]
        ss = s11 * s22 - abs2(s12)
        t11, t12, t22 = t[1,1], t[2,1], t[1,2]
        tt = t11 * t22 - abs2(t12)
        st = t11 * s22 + s11 * t22 - realprod(t12, s12) * 2
        iszero(ss) && return pushres!(res, tt / st, a, b)
        if iszero(tt)
            pushres!(res, T(0), a, b)
            return pushres!(res, st / ss, a, b)
        end
        u = st / ss / 2
        v = tt / ss
        d = u^2 - v
        d = sqrt(d)
        s = u + copysign(d, u)
        if u >= 0
            pushres!(res, v / s, a, b)
            pushres!(res, s, a, b)
        else
            pushres!(res, s, a, b)
            pushres!(res, v / s, a, b)
        end
    end
    res, [], []
end

function merge2(t, s, a, b)
    n = size(t, 2)
    n <= 2 && return merge12(t, s, a, b)
    #println("merge($n)")
    T = real(eltype(t))
    n2 = (n + 1) ÷ 2
    res1, = merge2(view(t, :, 1:n2), s, a, b)
    res2, = merge2(view(t, :, n2+1:n), view(s, :, n2+1:n), a, b)
    μ = sort([res1; res2])

    μ, ka, offseta, offsetb = starting_points!(μ, t, s, a, b)
    p = length(μ) - offseta - offsetb
    λ = Vector{T}(undef, p)
    bcount = zeros(Int, p)
    lcount = zeros(Int, p)

    for i = 1 : p
        j = i + ka
        jj = i + offseta
        x = μ[jj]
        sig = κ(x, t, s) < j ? 1 : -1
        a, b = interval(μ, jj)
        mult = estimate_mult(x, sig, μ, jj)
        λ[i], bcount[i], lcount[i] = iteration_one(a, x, b, j, n, mult, t, s)
        #println("merge($n): $(λ[i]) <- iterate ($a <= $x < $b) j = $j isig = $sig, mult = $mult" )
    end
    λ, bcount, lcount
end

function bisect(t, s, a, b, x, dx, j)
    while b - a > dx
        k = κ(x, t, s)
        if k >= j
            b = x
        else
            a = x
        end
        if isfinite(a) && isfinite(b)
            x = a + (b - a) / 2
        else
            x = 2 * x
        end
    end
    a, b
end

function bounds(t, s, scale::T, dx, j) where T
    x = zero(T)
    if κ(x, t, s) >= j
        a = T(-Inf)
        b = x
        x = -scale
        lb, ub = bisect(t, s, a, b, x, dx, j)
    else
        a = x
        b = T(Inf)
        x = scale
        lb, ub = bisect(t, s, a, b, x, dx, j)
    end
    lb, ub
end

function bounds(t, s)
    T = real(eltype(t))
    n = size(t, 2)
    scale = norm(t) / norm(s)
    if !isfinite(scale)
        scale = oneunit(T)
    end
    dx = scale / 4
    lb, _ = bounds(t, s, scale, dx, 1)
    _, ub = bounds(t, s, scale, dx, n)
    lb, ub
end

function zrst(t, s, a=-Inf, b=Inf)
    lb, ub = bounds(t, s)
    aa, bb = max(a, lb), min(b, ub)
    merge2(t, s, aa, bb)
end

function correct(t, s, x, dx)
    n = size(t, 2)
    r = 1
    a = x - dx
    b = x + dx
    j1 = κ(a, t, s)
    j2 = κ(b, t, s)
    j = j2
    iteration_one(a, x, b, j, n, r, t, s)
end


text:
module PSDMatrices

using LinearAlgebra
import LinearAlgebra: cholesky, diag, det, logdet

import Base: +, *, /, \, ==, inv, Matrix, kron, copy

# Write your package code here.
include("psdmatrix.jl")

export PSDMatrix

end

ext: jl


text:
# Autogenerated wrapper script for YASM_jll for armv6l-linux-gnueabihf
export vsyasm, yasm, ytasm

using NASM_jll
JLLWrappers.@generate_wrapper_header("YASM")
JLLWrappers.@declare_executable_product(vsyasm)
JLLWrappers.@declare_executable_product(yasm)
JLLWrappers.@declare_executable_product(ytasm)
function __init__()
    JLLWrappers.@generate_init_header(NASM_jll)
    JLLWrappers.@init_executable_product(
        vsyasm,
        "bin/vsyasm",
    )

    JLLWrappers.@init_executable_product(
        yasm,
        "bin/yasm",
    )

    JLLWrappers.@init_executable_product(
        ytasm,
        "bin/ytasm",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#######################################################################
#
# An example of creating Excel Stock charts with XlsxWriter.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_stock.xlsx")
	ws = add_worksheet!(wb)
	bold = add_format!(wb, Dict("bold"=> 1))

	date_format = add_format!(wb, Dict("num_format"=> "dd/mm/yyyy"))

	chart = add_chart!(wb, Dict("type"=> "stock"))

	# Add the worksheet data that the charts will refer to.
	headings = ["Date", "High", "Low", "Close"]
	data = [
		[Date(2007, 1, 1) Date(2007, 1, 2) Date(2007, 1, 3) Date(2007, 1, 4) Date(2007, 1, 5)];
		[27.2 25.03 19.05 20.34 18.5];
		[23.49 19.55 15.12 17.84 16.34];
		[25.45 23.05 17.32 20.45 17.34];
	]

	write_row!(ws, "A1", headings, fmt=bold)

	for row in 1:5
		write!(ws, row, 0, data[1, row], fmt=date_format)
		write!(ws, row, 1, data[2, row])
		write!(ws, row, 2, data[3, row])
		write!(ws, row, 3, data[4, row])
	end

	set_column!(ws, "A:D", 11)

	# Add a series for each of the High-Low-Close columns.
	add_series!(chart, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$6",
		"values"=> "=Sheet1!\$B\$2:\$B\$6",
	))

	add_series!(chart, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$6",
		"values"=>     "=Sheet1!\$C\$2:\$C\$6",
	))

	add_series!(chart, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$6",
		"values"=> "=Sheet1!\$D\$2:\$D\$6",
	))

	# Add a chart title and some axis labels.
	set_title!(chart, Dict("name"=> "High-Low-Close"))
	set_x_axis!(chart, Dict("name"=> "Date"))
	set_y_axis!(chart, Dict("name"=> "Share price"))

	insert_chart!(ws, "E9", chart)

	close(wb)
	isfile("chart_stock.xlsx")
end

test()


text:
##############################################################################
##
## Reshaping
##
## Also, see issue # ??
##
##############################################################################

##############################################################################
##
## stack()
## melt()
##
##############################################################################

"""
Stacks a DataFrame; convert from a wide to long format

```julia
stack(df::AbstractDataFrame, measure_vars, id_vars)
stack(df::AbstractDataFrame, measure_vars)
stack(df::AbstractDataFrame)
melt(df::AbstractDataFrame, id_vars, measure_vars)
melt(df::AbstractDataFrame, id_vars)
```

### Arguments

* `df` : the AbstractDataFrame to be stacked

* `measure_vars` : the columns to be stacked (the measurement
  variables), a normal column indexing type, like a Symbol,
  Vector{Symbol}, Int, etc.; for `melt`, defaults to all
  variables that are not `id_vars`

* `id_vars` : the identifier columns that are repeated during
  stacking, a normal column indexing type; for `stack` defaults to all
  variables that are not `measure_vars`

If neither `measure_vars` or `id_vars` are given, `measure_vars`
defaults to all floating point columns.

### Result

* `::DataFrame` : the long-format dataframe with column `:value`
  holding the values of the stacked columns (`measure_vars`), with
  column `:variable` a Vector of Symbols with the `measure_vars` name,
  and with columns for each of the `id_vars`.

See also `stackdf` and `meltdf` for stacking methods that return a
view into the original DataFrame. See `unstack` for converting from
long to wide format.

### Examples

```julia
d1 = DataFrame(a = repeat([1:3;], inner = [4]),
               b = repeat([1:4;], inner = [3]),
               c = randn(12),
               d = randn(12),
               e = map(string, 'a':'l'))

d1s = stack(d1, [:c, :d])
d1s2 = stack(d1, [:c, :d], [:a])
d1m = melt(d1, [:a, :b, :e])
```

"""
function stack(df::AbstractDataFrame, measure_vars::Vector{Int}, id_vars::Vector{Int})
    N = length(measure_vars)
    cnames = names(df)[id_vars]
    insert!(cnames, 1, :value)
    insert!(cnames, 1, :variable)
    DataFrame(Any[Compat.repeat(_names(df)[measure_vars], inner=nrow(df)),   # variable
                  vcat([df[c] for c in measure_vars]...),                    # value
                  [Compat.repeat(df[c], outer=N) for c in id_vars]...],      # id_var columns
              cnames)
end
function stack(df::AbstractDataFrame, measure_vars::Int, id_vars::Int)
    stack(df, [measure_vars], [id_vars])
end
function stack(df::AbstractDataFrame, measure_vars::Vector{Int}, id_vars::Int)
    stack(df, measure_vars, [id_vars])
end
function stack(df::AbstractDataFrame, measure_vars::Int, id_vars::Vector{Int})
    stackdf(df, [measure_vars], id_vars)
end
stack(df::AbstractDataFrame, measure_vars, id_vars) =
    stack(df, index(df)[measure_vars], index(df)[id_vars])
function stack(df::AbstractDataFrame, measure_vars)
    mv_inds = index(df)[measure_vars]
    stack(df, mv_inds, _setdiff(1:ncol(df), mv_inds))
end
function stack(df::AbstractDataFrame)
    idx = [1:length(df);][[t <: AbstractFloat for t in eltypes(df)]]
    stack(df, idx)
end

"""
Stacks a DataFrame; convert from a wide to long format; see
`stack`.
"""
melt(df::AbstractDataFrame, id_vars::@compat(Union{Int,Symbol})) = melt(df, [id_vars])
function melt(df::AbstractDataFrame, id_vars)
    id_inds = index(df)[id_vars]
    stack(df, _setdiff(1:ncol(df), id_inds), id_inds)
end
melt(df::AbstractDataFrame, id_vars, measure_vars) = stack(df, measure_vars, id_vars)
melt(df::AbstractDataFrame) = stack(df)

##############################################################################
##
## unstack()
##
##############################################################################

"""
Unstacks a DataFrame; convert from a long to wide format

```julia
unstack(df::AbstractDataFrame, rowkey, colkey, value)
unstack(df::AbstractDataFrame, colkey, value)
unstack(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame to be unstacked

* `rowkey` : the column with a unique key for each row, if not given,
  find a key by grouping on anything not a `colkey` or `value`

* `colkey` : the column holding the column names in wide format,
  defaults to `:variable`

* `value` : the value column, defaults to `:value`

### Result

* `::DataFrame` : the wide-format dataframe

### Examples

```julia
wide = DataFrame(id = 1:12,
                 a  = repeat([1:3;], inner = [4]),
                 b  = repeat([1:4;], inner = [3]),
                 c  = randn(12),
                 d  = randn(12))

long = stack(wide)
wide0 = unstack(long)
wide1 = unstack(long, :variable, :value)
wide2 = unstack(long, :id, :variable, :value)
```
Note that there are some differences between the widened results above.

"""
function unstack(df::AbstractDataFrame, rowkey::Int, colkey::Int, value::Int)
    # `rowkey` integer indicating which column to place along rows
    # `colkey` integer indicating which column to place along column headers
    # `value` integer indicating which column has values
    refkeycol = PooledDataArray(df[rowkey])
    valuecol = df[value]
    # TODO make a version with a default refkeycol
    keycol = PooledDataArray(df[colkey])
    Nrow = length(refkeycol.pool)
    Ncol = length(keycol.pool)
    # TODO make fillNA(type, length)
    payload = DataFrame(Any[DataArray(eltype(valuecol), Nrow) for i in 1:Ncol], map(Symbol, keycol.pool))
    nowarning = true
    for k in 1:nrow(df)
        j = @compat Int(keycol.refs[k])
        i = @compat Int(refkeycol.refs[k])
        if i > 0 && j > 0
            if nowarning && !isna(payload[j][i])
                warn("Duplicate entries in unstack.")
                nowarning = false
            end
            payload[j][i]  = valuecol[k]
        end
    end
    insert!(payload, 1, refkeycol.pool, _names(df)[rowkey])
end
unstack(df::AbstractDataFrame, rowkey, colkey, value) =
    unstack(df, index(df)[rowkey], index(df)[colkey], index(df)[value])

# Version of unstack with just the colkey and value columns provided
unstack(df::AbstractDataFrame, colkey, value) =
    unstack(df, index(df)[colkey], index(df)[value])

function unstack(df::AbstractDataFrame, colkey::Int, value::Int)
    # group on anything not a key or value:
    g = groupby(df, setdiff(_names(df), _names(df)[[colkey, value]]))
    groupidxs = [g.idx[g.starts[i]:g.ends[i]] for i in 1:length(g.starts)]
    rowkey = PooledDataArray(zeros(Int, size(df, 1)), [1:length(groupidxs);])
    for i in 1:length(groupidxs)
        rowkey[groupidxs[i]] = i
    end
    keycol = PooledDataArray(df[colkey])
    valuecol = df[value]
    df1 = df[g.idx[g.starts], g.cols]
    keys = unique(keycol)
    Nrow = length(g)
    Ncol = length(keycol.pool)
    df2 = DataFrame(Any[DataArray(fill(valuecol[1], Nrow), fill(true, Nrow)) for i in 1:Ncol], map(@compat(Symbol), keycol.pool))
    nowarning = true
    for k in 1:nrow(df)
        j = @compat Int(keycol.refs[k])
        i = rowkey[k]
        if i > 0 && j > 0
            if nowarning && !isna(df2[j][i])
                warn("Duplicate entries in unstack.")
                nowarning = false
            end
            df2[j][i]  = valuecol[k]
        end
    end
    hcat(df1, df2)
end

unstack(df::AbstractDataFrame) = unstack(df, :id, :variable, :value)

##############################################################################
##
## Reshaping using referencing (issue #145)
## New AbstractVector types (all read only):
##     StackedVector
##     RepeatedVector
##
##############################################################################

"""
An AbstractVector{Any} that is a linear, concatenated view into
another set of AbstractVectors

NOTE: Not exported.

### Constructor

```julia
RepeatedVector(d::AbstractVector...)
```

### Arguments

* `d...` : one or more AbstractVectors

### Examples

```julia
StackedVector(Any[[1,2], [9,10], [11,12]])  # [1,2,9,10,11,12]
```

"""
type StackedVector <: AbstractVector{Any}
    components::Vector{Any}
end

function Base.getindex(v::StackedVector,i::Real)
    lengths = [length(x)::Int for x in v.components]
    cumlengths = [0; cumsum(lengths)]
    j = searchsortedlast(cumlengths .+ 1, i)
    if j > length(cumlengths)
        error("indexing bounds error")
    end
    k = i - cumlengths[j]
    if k < 1 || k > length(v.components[j])
        error("indexing bounds error")
    end
    v.components[j][k]
end

function Base.getindex{I<:Real}(v::StackedVector,i::AbstractVector{I})
    result = similar(v.components[1], length(i))
    for idx in 1:length(i)
        result[idx] = v[i[idx]]
    end
    result
end

Base.size(v::StackedVector) = (length(v),)
Base.length(v::StackedVector) = sum(map(length, v.components))
Base.ndims(v::StackedVector) = 1
Base.eltype(v::StackedVector) = promote_type(map(eltype, v.components)...)
Base.similar(v::StackedVector, T, dims::Dims) = similar(v.components[1], T, dims)

DataArrays.PooledDataArray(v::StackedVector) = PooledDataArray(v[:]) # could be more efficient

"""
An AbstractVector that is a view into another AbstractVector with
repeated elements

NOTE: Not exported.

### Constructor

```julia
RepeatedVector(parent::AbstractVector, inner::Int, outer::Int)
```

### Arguments

* `parent` : the AbstractVector that's repeated
* `inner` : the numer of times each element is repeated
* `outer` : the numer of times the whole vector is repeated after
  expanded by `inner`

`inner` and `outer` have the same meaning as similarly named arguments
to `repeat`.

### Examples

```julia
RepeatedVector([1,2], 3, 1)   # [1,1,1,2,2,2]
RepeatedVector([1,2], 1, 3)   # [1,2,1,2,1,2]
RepeatedVector([1,2], 2, 2)   # [1,2,1,2,1,2,1,2]
```

"""
type RepeatedVector{T} <: AbstractVector{T}
    parent::AbstractVector{T}
    inner::Int
    outer::Int
end

function Base.getindex{T,I<:Real}(v::RepeatedVector{T},i::AbstractVector{I})
    N = length(v.parent)
    idx = Int[Base.fld1(mod1(j,v.inner*N),v.inner) for j in i]
    v.parent[idx]
end
function Base.getindex{T}(v::RepeatedVector{T},i::Real)
    N = length(v.parent)
    idx = Base.fld1(mod1(i,v.inner*N),v.inner)
    v.parent[idx]
end
Base.getindex(v::RepeatedVector,i::Range) = getindex(v, [i;])

Base.size(v::RepeatedVector) = (length(v),)
Base.length(v::RepeatedVector) = v.inner * v.outer * length(v.parent)
Base.ndims(v::RepeatedVector) = 1
Base.eltype{T}(v::RepeatedVector{T}) = T
Base.reverse(v::RepeatedVector) = RepeatedVector(reverse(v.parent), v.inner, v.outer)
Base.similar(v::RepeatedVector, T, dims::Dims) = similar(v.parent, T, dims)
Base.unique(v::RepeatedVector) = unique(v.parent)

function DataArrays.PooledDataArray(v::RepeatedVector)
    res = DataArrays.PooledDataArray(v.parent)
    res.refs = repeat(res.refs, inner = [v.inner], outer = [v.outer])
    res
end

##############################################################################
##
## stackdf()
## meltdf()
## Reshaping using referencing (issue #145), using the above vector types
##
##############################################################################

"""
A stacked view of a DataFrame (long format)

Like `stack` and `melt`, but a view is returned rather than data
copies.

```julia
stackdf(df::AbstractDataFrame, measure_vars, id_vars)
stackdf(df::AbstractDataFrame, measure_vars)
meltdf(df::AbstractDataFrame, id_vars, measure_vars)
meltdf(df::AbstractDataFrame, id_vars)
```

### Arguments

* `df` : the wide AbstractDataFrame

* `measure_vars` : the columns to be stacked (the measurement
  variables), a normal column indexing type, like a Symbol,
  Vector{Symbol}, Int, etc.; for `melt`, defaults to all
  variables that are not `id_vars`

* `id_vars` : the identifier columns that are repeated during
  stacking, a normal column indexing type; for `stack` defaults to all
  variables that are not `measure_vars`

### Result

* `::DataFrame` : the long-format dataframe with column `:value`
  holding the values of the stacked columns (`measure_vars`), with
  column `:variable` a Vector of Symbols with the `measure_vars` name,
  and with columns for each of the `id_vars`.

The result is a view because the columns are special AbstractVectors
that return indexed views into the original DataFrame.

### Examples

```julia
d1 = DataFrame(a = repeat([1:3;], inner = [4]),
               b = repeat([1:4;], inner = [3]),
               c = randn(12),
               d = randn(12),
               e = map(string, 'a':'l'))

d1s = stackdf(d1, [:c, :d])
d1s2 = stackdf(d1, [:c, :d], [:a])
d1m = meltdf(d1, [:a, :b, :e])
```

"""
function stackdf(df::AbstractDataFrame, measure_vars::Vector{Int}, id_vars::Vector{Int})
    N = length(measure_vars)
    cnames = names(df)[id_vars]
    insert!(cnames, 1, :value)
    insert!(cnames, 1, :variable)
    DataFrame(Any[RepeatedVector(_names(df)[measure_vars], nrow(df), 1),   # variable
                  StackedVector(Any[df[:,c] for c in measure_vars]),     # value
                  [RepeatedVector(df[:,c], 1, N) for c in id_vars]...],     # id_var columns
              cnames)
end
function stackdf(df::AbstractDataFrame, measure_vars::Int, id_vars::Int)
    stackdf(df, [measure_vars], [id_vars])
end
function stackdf(df::AbstractDataFrame, measure_vars, id_vars::Int)
    stackdf(df, measure_vars, [id_vars])
end
function stackdf(df::AbstractDataFrame, measure_vars::Int, id_vars)
    stackdf(df, [measure_vars], id_vars)
end
function stackdf(df::AbstractDataFrame, measure_vars, id_vars)
    stackdf(df, index(df)[measure_vars], index(df)[id_vars])
end
function stackdf(df::AbstractDataFrame, measure_vars)
    m_inds = index(df)[measure_vars]
    stackdf(df, m_inds, _setdiff(1:ncol(df), m_inds))
end
function stackdf(df::AbstractDataFrame)
    idx = [1:length(df);][[t <: AbstractFloat for t in eltypes(df)]]
    stackdf(df, idx)
end

"""
A stacked view of a DataFrame (long format); see `stackdf`
"""
function meltdf(df::AbstractDataFrame, id_vars)
    id_inds = index(df)[id_vars]
    stackdf(df, _setdiff(1:ncol(df), id_inds), id_inds)
end
meltdf(df::AbstractDataFrame, id_vars, measure_vars) =
    stackdf(df, measure_vars, id_vars)
meltdf(df::AbstractDataFrame) = stackdf(df)


text:
# Visualize results of Bayesian experiments (average over many different problems) that are saved in the results folder (names should match)

using Statistics
using PyPlot
using HDF5
using DataFrames
using CSV

# OPTIONAL: write results in a csv file
updateCSV = "on"

# NAME OF THE INSTANCE (should match the stored data)
PB = "Bayesian"
tweak = 0.5
T = 20000
N = 100
K = 5
nbBreaks = 6

fname="results/$(PB)_tweak_$(tweak)_"

# POLICIES TO INCLUDE
policies=["klUCB","EXP3SOpt","SWKLUCB","DKLUCB","MKLUCB2","CUSUMKLUCB2","GLRKLUCB","GLRKLUCBLoc"]
names =["klUCB","EXP3.S","SW-klUCB","D-klUCB","M-klUCB","CUSUM-klUCB","GLR-klUCB","GLR-klUCB-Local"]

lP = length(policies)

df=DataFrame(Algo= String[],Regret=Float32[],Stdev=Float32[],NbRestarts=Float32[])
clf()

# display results and plot the Bayesian regret
for imeth in 1:lP
    policy = policies[imeth]
    policyname = names[imeth]
    name = "$(fname)_$(policy)_T_$(T)_N_$(N)_K_$(K)_U_$(nbBreaks)"
    Regret=h5read(name,"Regret")
    Restarts=h5read(name,"Restarts")
    tsave=floor.(Int,collect(range(1,T,length=500)))
    # print
    regfinal = mean(Regret[:,end])
    stdevfinal = std(Regret[:,end])
    restarts = mean(Restarts)
    push!(df,("$(policy)",regfinal,stdevfinal,restarts))
    print("Mean final regret for $(policyname) is $(regfinal)\n")
    print("Standard deviation of the final regret for $(policyname) is $(stdevfinal))\n")
	print("Mean number of restarts for $(policyname) is $(restarts)\n\n")
	# plot
	plot(vec(tsave),vec(mean(Regret,dims=1)),label="$(policyname)")
end

legend()

if (updateCSV == "on")
    CSV.write("$(fname)T_$(T)_N_$(N)_K_$(K)_U_$(nbBreaks).csv",df)
end


text:
medium=Medium(:acou_homo2D,5)
update!(medium, [:vp,:rho], randn_perc=0.1)
ageom=AGeom(medium.mgrid, SSrcs(4), Srcs(1), Recs(100))
update!(ageom, SSrcs(),[0,0],990.0,[0, 2π])
update!(ageom, Recs(),[0,0],990.0,[0, 2π])

wav, tgrid=ricker(medium, 3, 0.8)

for field in [:p, :vx,:vz]
	println("############ Testing Backprop for source type ", field)
	srcwav = SrcWav(tgrid, ageom, [field])
	update!(srcwav, [field], wav)

	for sflags in [[1,-1],[2,-2]]
		pa=SeisForwExpt(FdtdAcou(),npw=1, tgrid=tgrid,
		#	pml_edges=[:null],
			gmodel_flag=false,
			sflags=[sflags[1]],
			# snaps_field=true,
			verbose=true,
			backprop_flag=1,
			illum_flag=true,ageom=[ageom], srcwav=[srcwav],
			medium=medium);

		update!(pa);
		rec1=deepcopy(pa.c.data[1])

		# change source flag and update wavelets in pa
		pa.c.sflags=[sflags[2]];
		GeoPhyInv.update!(pa,[srcwav])
		pa.c.backprop_flag=-1 # do backpropagation

		update!(pa)
		rec2=deepcopy(pa.c.data[1])

		# time reverse
		reverse!(rec2);

		# compare results
		# least-squares misfit
		paerr=GeoPhyInv.VNamedD_misfit(rec1, rec2)
		err = GeoPhyInv.func_grad!(paerr)

		# normalized error
		error = err[1]./paerr.ynorm

		# desired accuracy?
		@test error<1e-20
	end
end


text:
immutable FIRFilter
	filtobjPtr::Ptr{Void}
end

function Base.display(filt::FIRFilter)
	ccall((:firfilt_crcf_print, "libliquid"), Void, (Ptr{Void}, ), filt.filtobjPtr)
end

function createfirfilter(h::Array{Complex64})
	ms = ccall((:firfilt_crcf_create, "libliquid"), Ptr{Void}, (Ptr{Complex64},Int), h,length(h));
	return FIRFilter(ms);

end

function recreatefirfilter()
end


text:
export BrickletPiezoSpeakerV2Beep
struct BrickletPiezoSpeakerV2Beep
    frequency::Integer
    volume::Integer
    duration::Integer
    duration_remaining::Integer
end

export BrickletPiezoSpeakerV2Alarm
struct BrickletPiezoSpeakerV2Alarm
    start_frequency::Integer
    end_frequency::Integer
    step_size::Integer
    step_delay::Integer
    volume::Integer
    duration::Integer
    duration_remaining::Integer
    current_frequency::Integer
end

export BrickletPiezoSpeakerV2SPITFPErrorCount
struct BrickletPiezoSpeakerV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletPiezoSpeakerV2Identity
struct BrickletPiezoSpeakerV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletPiezoSpeakerV2
"""
Creates beep and alarm with configurable volume and frequency
"""
mutable struct BrickletPiezoSpeakerV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletPiezoSpeakerV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_piezo_speaker_v2")
        deviceInternal = package.BrickletPiezoSpeakerV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_beep
"""
    $(SIGNATURES)

Beeps with the given frequency and volume for the duration.

A duration of 0 stops the current beep if any is ongoing.
A duration of 4294967295 results in an infinite beep.
"""
function set_beep(device::BrickletPiezoSpeakerV2, frequency, volume, duration)
    device.deviceInternal.set_beep(frequency, volume, duration)
end

export get_beep
"""
    $(SIGNATURES)

Returns the last beep settings as set by :func:`Set Beep`. If a beep is currently
running it also returns the remaining duration of the beep.

If the frequency or volume is updated during a beep (with :func:`Update Frequency`
or :func:`Update Volume`) this function returns the updated value.
"""
function get_beep(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.get_beep()
end

export set_alarm
"""
    $(SIGNATURES)

Creates an alarm (a tone that goes back and force between two specified frequencies).

The following parameters can be set:

* Start Frequency: Start frequency of the alarm.
* End Frequency: End frequency of the alarm.
* Step Size: Size of one step of the sweep between the start/end frequencies.
* Step Delay: Delay between two steps (duration of time that one tone is used in a sweep).
* Duration: Duration of the alarm.

A duration of 0 stops the current alarm if any is ongoing.
A duration of 4294967295 results in an infinite alarm.

Below you can find two sets of example settings that you can try out. You can use
these as a starting point to find an alarm signal that suits your application.

Example 1: 10 seconds of loud annoying fast alarm

* Start Frequency = 800
* End Frequency = 2000
* Step Size = 10
* Step Delay = 1
* Volume = 10
* Duration = 10000

Example 2: 10 seconds of soft siren sound with slow build-up

* Start Frequency = 250
* End Frequency = 750
* Step Size = 1
* Step Delay = 5
* Volume = 0
* Duration = 10000

The following conditions must be met:

* Start Frequency: has to be smaller than end frequency
* End Frequency: has to be bigger than start frequency
* Step Size: has to be small enough to fit into the frequency range
* Step Delay: has to be small enough to fit into the duration
"""
function set_alarm(device::BrickletPiezoSpeakerV2, start_frequency, end_frequency, step_size, step_delay, volume, duration)
    device.deviceInternal.set_alarm(start_frequency, end_frequency, step_size, step_delay, volume, duration)
end

export get_alarm
"""
    $(SIGNATURES)

Returns the last alarm settings as set by :func:`Set Alarm`. If an alarm is currently
running it also returns the remaining duration of the alarm as well as the
current frequency of the alarm.

If the volume is updated during an alarm (with :func:`Update Volume`)
this function returns the updated value.
"""
function get_alarm(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.get_alarm()
end

export update_volume
"""
    $(SIGNATURES)

Updates the volume of an ongoing beep or alarm.
"""
function update_volume(device::BrickletPiezoSpeakerV2, volume)
    device.deviceInternal.update_volume(volume)
end

export update_frequency
"""
    $(SIGNATURES)

Updates the frequency of an ongoing beep.
"""
function update_frequency(device::BrickletPiezoSpeakerV2, frequency)
    device.deviceInternal.update_frequency(frequency)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletPiezoSpeakerV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletPiezoSpeakerV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletPiezoSpeakerV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletPiezoSpeakerV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletPiezoSpeakerV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletPiezoSpeakerV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletPiezoSpeakerV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletPiezoSpeakerV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using LinearAlgebra
using CVChannel
using Test
using DelimitedFiles

"""
In this script we look at the PPT relaxation of the communication value of the
Werner-Holevo channel for multiple copies and its relation to the parameters
λ and the number of copies.
"""

println("Here we look at how the PPT relaxation of cv scales in number of copies.")
println("We do this for n ≤ 8 and λ in [0,0.1,...,1].") #It takes a lot more time if n = 9 & λ != 0
println("This also allows us to generate our non-multiplicativity results.")
@testset "Generate cv as function of n and lambda" begin
    d = 3
    n_list = [1:8;]
    λ_list = [0:0.01:1;]
    data_table = zeros(length(n_list),length(λ_list))
    for n_val in n_list
        println("Now obtaining values for n = ", n_val)
        λ_ctr = 1
        for λ_val in λ_list
            if λ_ctr == Int(ceil(length(λ_list)/2))
                println("Now halfway through this value of n.")
            end
            A,B,g,a = generalWHLPConstraints(n_val,d,λ_val*ones(n_val))
            cv_ppt, opt_var = wernerHolevoCVPPT(n_val,d,A,B,g,a)
            data_table[n_val,λ_ctr] = cv_ppt
            λ_ctr += 1
        end
    end
    #This generates the table about multiplicativity
    mult_table =[
        (data_table[2,:] - data_table[1,:].^2)';
        (data_table[4,:] - data_table[2,:].^2)';
        (data_table[8,:] - data_table[4,:].^2)';
    ]
    λ_header = zeros(1,length(λ_list))
    spacing = ""
    for i in 1:length(λ_list)
        λ_header[i] = λ_list[i]
        spacing = [spacing ""]
    end
    header_p_1 = hcat("n\\lambda:", λ_header)
    header_p_2 = hcat("lambda:", λ_header)
    col_inf_p_2 = ["cv(N^2) - cv(N)^2";"cv(N^4) - cv(N^2)^2";"cv(N^8) - cv(N^4)^2"]
    data_to_save_p_1 = vcat(hcat(n_list, data_table),spacing)
    data_to_save_p_2 = vcat(header_p_2,hcat(col_inf_p_2,mult_table))
    data_to_save = vcat(header_p_1,vcat(data_to_save_p_1,data_to_save_p_2))
    println("Here are the results:")
    show(stdout, "text/plain", data_to_save)

    info_vec = Vector{Union{Nothing,String}}(nothing, size(data_to_save)[1])
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by multi-copy-werner-holevo-cv.jl"
    info_vec[3] = "d = " * string(d)
    info_vec[4] = "n_list = " * string(n_list)
    info_vec[5] = "lambda_list = " * string(λ_list)
    data_to_save = hcat(data_to_save,info_vec)
    print("\nPlease name the file you'd like to write the results to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')

    @test isapprox(data_table[2,1], 3, atol = 1e-6)
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
using StringRepFilter
using Test

@testset "StringRepFilter.jl" begin

    ## --------------------------------------------------------------
    @info("Testing has_match")

    struct Obj; A end

    for d in [Dict(:A => "AAB"), (;A = "AAB"), Obj("AAB")]

        @test has_match(d, "B")
        @test has_match(d, Dict, :A) #  dict or has property A
        @test has_match(d, :A => String)
        @test has_match(d, :A => [(x) -> length(x) == 3])

        @test has_match(d, :A => "B")
        @test has_match(d, :A => r"^A")
        @test has_match(d, :A => r"B$")
        @test !has_match(d, :A => r"C")
        @test !has_match(d, :A => r"^B")
        @test !has_match(d, :B => r"^A")

        @test has_match(d, [:A => r"B$"])
        @test !has_match(d, [:A => r"^B"])
        @test !has_match(d, [:B => r"^A"])

        if (d isa Dict)
            # Check wrong key type
            @test !has_match(d, "A" => r"^A")
            @test !has_match(d, "A" => r"^B")
        end

    end

    ## --------------------------------------------------------------
    @info("Testing finders and filters")

    col = ["AAA", "BBB", 123]
    @test all(filter_match(col, r"^[A|1]") .== filter_match(col, "A", "1"))
    @test findfirst_match(col, r"^[A|1]") == (1, first(col))
    @test findall_match(col, r"^[A|1]") == [1, 3]

    ## --------------------------------------------------------------
end


text:
module MPIFiles

using UUIDs
using Graphics: @mustimplement

using Reexport
@reexport using SparsityOperators
@reexport using AxisArrays
const axes = Base.axes
@reexport using Interpolations
@reexport using HDF5
@reexport using Dates
@reexport using DelimitedFiles
@reexport using Images
@reexport using LinearAlgebra
@reexport using Random
@reexport using Mmap
@reexport using Statistics
@reexport using Unitful
using CodecZlib
using Tar
using Pkg.PlatformEngines
using Pkg.GitTools
using Pkg.Artifacts
using Unitful
using InteractiveUtils
using UnitfulAngles
using Inflate, SHA
using StableRNGs
using REPL: fielddoc
using DocStringExtensions

if VERSION < v"1.1"
  isnothing(x) = x == nothing
end

### global import list ###

import Base: convert, get, getindex, haskey, iterate, length, ndims, range, read, show, time, write, close
import FileIO: save
import HDF5: h5read
import Interpolations: interpolate

### export list ###

export MPIFile

# general parameters
export version, uuid

# study parameters
export studyName, studyNumber, studyUuid, studyDescription, studyTime

# experiment parameters
export experimentName, experimentNumber, experimentUuid, experimentDescription, experimentSubject,
      experimentIsSimulation, experimentIsCalibration,
      experimentHasMeasurement, experimentHasReconstruction

# tracer parameters
export tracerName, tracerBatch, tracerVolume, tracerConcentration,
       tracerSolute, tracerInjectionTime, tracerVendor

# scanner parameters
export scannerFacility, scannerOperator, scannerManufacturer, scannerName,
       scannerTopology

# acquisition parameters
export acqStartTime, acqNumFrames, acqNumAverages,
       acqGradient, acqOffsetField, acqNumPeriodsPerFrame, acqSize

# drive-field parameters
export dfNumChannels, dfStrength, dfPhase, dfBaseFrequency, dfCustomWaveform,
       dfDivider, dfWaveform, dfCycle

# receiver parameters
export rxNumChannels, rxBandwidth, rxNumSamplingPoints,
       rxTransferFunction, rxTransferFunctionFileName, rxHasTransferFunction, rxUnit,
       rxDataConversionFactor, rxInductionFactor

# measurements
export measData, measDataTDPeriods, measIsFourierTransformed, measIsTFCorrected,
       measIsTransferFunctionCorrected,
       measIsBGCorrected, measIsBackgroundCorrected, measIsFastFrameAxis,
       measIsFramePermutation, measIsFrequencySelection,
       measIsBGFrame, measIsBackgroundFrame, measIsSpectralLeakageCorrected, measFramePermutation,
       measFrequencySelection, measIsSparsityTransformed, measIsCalibProcessed

# calibrations
export calibSNR, calibSnr, calibFov, calibFieldOfView, calibFovCenter,
       calibFieldOfViewCenter, calibSize, calibOrder, calibPositions,
       calibOffsetField, calibDeltaSampleSize,
       calibMethod, calibIsMeanderingGrid

# reconstruction results
export recoData, recoFov, recoFieldOfView, recoFovCenter, recoFieldOfViewCenter,
       recoSize, recoOrder, recoPositions

# additional functions that should be implemented by an MPIFile
export filepath, systemMatrixWithBG, systemMatrix

export selectedChannels
### Interface of an MPIFile ###

abstract type MPIFile end

# general parameters
@mustimplement version(f::MPIFile)
@mustimplement uuid(f::MPIFile)
@mustimplement time(f::MPIFile)

# study parameters
@mustimplement studyName(f::MPIFile)
@mustimplement studyNumber(f::MPIFile)
@mustimplement studyUuid(f::MPIFile)
@mustimplement studyDescription(f::MPIFile)
@mustimplement studyTime(f::MPIFile)

# experiment parameters
@mustimplement experimentName(f::MPIFile)
@mustimplement experimentNumber(f::MPIFile)
@mustimplement experimentUuid(f::MPIFile)
@mustimplement experimentDescription(f::MPIFile)
@mustimplement experimentSubject(f::MPIFile)
@mustimplement experimentIsSimulation(f::MPIFile)
@mustimplement experimentIsCalibration(f::MPIFile)
@mustimplement experimentHasReconstruction(f::MPIFile)
@mustimplement experimentHasMeasurement(f::MPIFile)

# tracer parameters
@mustimplement tracerName(f::MPIFile)
@mustimplement tracerBatch(f::MPIFile)
@mustimplement tracerVolume(f::MPIFile)
@mustimplement tracerConcentration(f::MPIFile)
@mustimplement tracerSolute(f::MPIFile)
@mustimplement tracerInjectionTime(f::MPIFile)

# scanner parameters
@mustimplement scannerBoreSize(f::MPIFile)
@mustimplement scannerFacility(f::MPIFile)
@mustimplement scannerOperator(f::MPIFile)
@mustimplement scannerManufacturer(f::MPIFile)
@mustimplement scannerName(f::MPIFile)
@mustimplement scannerTopology(f::MPIFile)

# acquisition parameters
@mustimplement acqStartTime(f::MPIFile)
@mustimplement acqNumAverages(f::MPIFile)
@mustimplement acqNumPeriodsPerFrame(f::MPIFile)
@mustimplement acqNumFrames(f::MPIFile)
@mustimplement acqGradient(f::MPIFile)
@mustimplement acqOffsetField(f::MPIFile)

# drive-field parameters
@mustimplement dfNumChannels(f::MPIFile)
@mustimplement dfStrength(f::MPIFile)
@mustimplement dfPhase(f::MPIFile)
@mustimplement dfBaseFrequency(f::MPIFile)
@mustimplement dfCustomWaveform(f::MPIFile)
@mustimplement dfDivider(f::MPIFile)
@mustimplement dfWaveform(f::MPIFile)
@mustimplement dfCycle(f::MPIFile)

# receiver properties
@mustimplement rxNumChannels(f::MPIFile)
@mustimplement rxBandwidth(f::MPIFile)
@mustimplement rxNumSamplingPoints(f::MPIFile)
@mustimplement rxTransferFunction(f::MPIFile)
@mustimplement rxTransferFunctionFileName(f::MPIFile)
@mustimplement rxHasTransferFunction(f::MPIFile)
@mustimplement rxInductionFactor(f::MPIFile)
@mustimplement rxUnit(f::MPIFile)
@mustimplement rxDataConversionFactor(f::MPIFile)

# measurements
@mustimplement measData(f::MPIFile)
@mustimplement measDataTD(f::MPIFile)
@mustimplement measDataTDPeriods(f::MPIFile, periods)
@mustimplement measFramePermutation(f::MPIFile)
@mustimplement measFrequencySelection(f::MPIFile)
@mustimplement measIsBGCorrected(f::MPIFile)
@mustimplement measIsBGFrame(f::MPIFile)
@mustimplement measIsFastFrameAxis(f::MPIFile)
@mustimplement measIsFourierTransformed(f::MPIFile)
@mustimplement measIsFramePermutation(f::MPIFile)
@mustimplement measIsFrequencySelection(f::MPIFile)
@mustimplement measIsSparsityTransformed(f::MPIFile)
@mustimplement measIsSpectralLeakageCorrected(f::MPIFile)
@mustimplement measIsTFCorrected(f::MPIFile)
@mustimplement measSparsityTransformation(f::MPIFile)
@mustimplement measSubsamplingIndices(f::MPIFile)
@mustimplement measIsCalibProcessed(b::MPIFile)

# calibrations
@mustimplement calibSNR(f::MPIFile)
@mustimplement calibFov(f::MPIFile)
@mustimplement calibFovCenter(f::MPIFile)
@mustimplement calibSize(f::MPIFile)
@mustimplement calibOrder(f::MPIFile)
@mustimplement calibPositions(f::MPIFile)
@mustimplement calibOffsetField(f::MPIFile)
@mustimplement calibDeltaSampleSize(f::MPIFile)
@mustimplement calibMethod(f::MPIFile)
@mustimplement calibIsMeanderingGrid(f::MPIFile)

# reconstruction results
@mustimplement recoData(f::MPIFile)
@mustimplement recoFov(f::MPIFile)
@mustimplement recoFovCenter(f::MPIFile)
@mustimplement recoSize(f::MPIFile)
@mustimplement recoOrder(f::MPIFile)
@mustimplement recoPositions(f::MPIFile)
@mustimplement recoIsOverscanRegion(f::MPIFile)

# additional functions that should be implemented by an MPIFile
@mustimplement filepath(f::MPIFile)

include("Derived.jl")
include("Custom.jl")
include("Utils.jl")
include("FramePermutation.jl")

### Concrete implementations ###
include("MDF.jl")
include("MDFInMemory.jl")
include("MDFCommon.jl")
include("Brukerfile.jl")
include("IMT.jl")

# This dispatches on the file extension and automatically
# generates the correct type
function MPIFile(filename::AbstractString; kargs...)
  filenamebase, ext = splitext(filename)
  if ext == ".mdf" || ext == ".hdf" || ext == ".h5"
    file = h5open(filename,"r")
    if haskey(file, "/version")
      return MDFFile(filename, file) # MDFFile currently has no kargs
    else
      return IMTFile(filename, file; kargs...)
    end
  else
    if isfile(joinpath(filename,"mdf"))
      filenameMDF = readline(joinpath(filename,"mdf"))
      return MDFFile(filenameMDF)
    else
      return BrukerFile(filename; kargs...)
    end
  end
end

function show(io::IO, f::MPIFile)
  print(io,supertype(typeof(f)))
  print(io,"\n\tStudy: ")
  show(io, studyName(f))
  print(io,", ")
  show(io,studyTime(f))
  print(io,"\n\tExperiment: ")
  show(io,experimentName(f))
  print(io,", ")
  show(io,acqStartTime(f))
  print(io,"\n")
end

# Opens a set of MPIFiles
function MPIFile(filenames::Vector)
  return map(x->MPIFile(x),filenames)
end

# For the do block
function MPIFile(h::Function, args...; kargs...)
  f = MPIFile(args...; kargs...)
  try
      h(f)
  finally
      close(f)
  end
end

Base.length(f::MPIFile) = 1
Base.close(f::MPIFile) = nothing

include("TransferFunction.jl")
include("MultiMPIFile.jl")
include("Measurements.jl")
include("SystemMatrix.jl")
include("FrequencyFilter.jl")
include("Conversion.jl")
include("RecoData.jl")
include("DatasetStore/DatasetStore.jl")
include("MixingFactors.jl")
include("Positions/Positions.jl")
include("MagneticFieldMeasurement.jl")

end # module


text:
function build_robot!(vis::Visualizer, model::WalledCartpole; name::Symbol=:WalledCartpole,
		r=0.05, α=1.0)
	r = convert(Float32, r)
	r_contact = convert(Float32, r * 1.5)

	wall_thickness = 0.15
	body_mat = MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, α))
	contact_mat = MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, α))
	wall_mat = MeshPhongMaterial(color = RGBA(0.7, 0.7, 0.7, 1.0))

	default_background!(vis)

	link = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.l), r/2)
	setobject!(vis[name][:robot]["link"], link, body_mat)
	vec_base = Vec(2r_contact, 4r, 1.5r)
	vec_shaft = Vec(3model.w + 2r_contact + 2wall_thickness, 0.5r, 0.5r)
	setobject!(vis[name][:robot]["linkbase"], Rect(-0.5 * vec_base, vec_base), body_mat)
	setobject!(vis[name][:robot]["shaft"], Rect(-0.5 * vec_shaft, vec_shaft), body_mat)
	setobject!(vis[name][:robot]["contact"], Sphere(Point3f0(0.0), r_contact), contact_mat)

	setobject!(vis[name][:env]["wall1"]["init"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness, wall_thickness, 0.75)), wall_mat)
	setobject!(vis[name][:env]["wall2"]["init"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness, wall_thickness, 0.75)), wall_mat)
	settransform!(vis[name][:env]["wall1"]["init"], Translation([-wall_thickness - model.w - 1.5 * r; -wall_thickness/2; 3r]))
	settransform!(vis[name][:env]["wall2"]["init"], Translation([model.w + 1.5 * r;                   -wall_thickness/2; 3r]))
	return nothing
end

function set_robot!(vis::Visualizer, model::WalledCartpole, q::AbstractVector;
		name::Symbol=:WalledCartpole)

	pt = [_kinematics(model, q, mode =  :tip)[1], 0.0, _kinematics(model, q, mode =  :tip)[2]]
    pb = [_kinematics(model, q, mode = :base)[1], 0.0, _kinematics(model, q, mode = :base)[2]]

	settransform!(vis[name][:robot]["link"], cable_transform(pt, pb))
	settransform!(vis[name][:robot]["linkbase"], Translation(pb))
	settransform!(vis[name][:robot]["contact"], Translation(pt))

	xw1 = q[3]
	xw2 = q[4]
	settransform!(vis[name][:env]["wall1"], Translation([xw1, 0.0, 0.0]))
	settransform!(vis[name][:env]["wall2"], Translation([xw2, 0.0, 0.0]))

	return nothing
end


text:
using Test
using Pluto
using Pluto: update_run!, ServerSession, ClientSession, Cell, Notebook

@testset "CellDepencencyVisualization" begin
    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = false

    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    notebook = Notebook([
                Cell("x = 1"), # prerequisite of test cell
                Cell("f(x) = x + y"), # depends on test cell
                Cell("f(3)"),

                Cell("""begin
                    g(a) = x
                    g(a,b) = y
                end"""), # depends on test cell
                Cell("y = x"), # test cell below
                Cell("g(6) + g(6,6)"),
                Cell("using Dates"),
            ])
    fakeclient.connected_notebook = notebook
    update_run!(🍭, notebook, notebook.cells)
    state = Pluto.notebook_to_js(notebook)

    id(i) = notebook.cells[i].cell_id

    order_of(i) = findfirst(isequal(id(i)), state["cell_execution_order"])
    @test order_of(7) < order_of(1) < order_of(5) < order_of(2) < order_of(3)

    deps = state["cell_dependencies"]

    @test deps[id(5)]["downstream_cells_map"] |> keys == Set(["y"])
    @test deps[id(5)]["downstream_cells_map"]["y"] == [id(2), id(4)]

    @test deps[id(5)]["upstream_cells_map"] |> keys == Set(["x"])
    @test deps[id(5)]["upstream_cells_map"]["x"] == [id(1)]

    # test if this also works for function definitions
    @test deps[id(2)]["downstream_cells_map"] |> keys == Set(["f"])
    @test deps[id(2)]["downstream_cells_map"]["f"] == [id(3)]

    @test deps[id(2)]["upstream_cells_map"] |> keys == Set(["y", "+"])
    @test deps[id(2)]["upstream_cells_map"]["y"] == [id(5)]
    @test deps[id(2)]["upstream_cells_map"]["+"] == [] # + function is not defined / extended in the notebook

    @test deps[id(1)]["precedence_heuristic"] > deps[id(7)]["precedence_heuristic"]
end


text:
using Documenter

push!(LOAD_PATH, "../src/")
using PointSpreadFunctions

DEPLOYDOCS = (get(ENV, "CI", nothing) == "true")

makedocs(;
         modules=[PointSpreadFunctions],
         format=Documenter.HTML(;
                                prettyurls = DEPLOYDOCS,
                                ),
         pages=[
             "Home" => "index.md",
             "Reference" => "library.md",
         ],
         repo="https://github.com/emmt/PointSpreadFunctions.jl/blob/{commit}{path}#L{line}",
         sitename = "Point spread functions for Julia",
         authors="Éric Thiébaut <eric.thiebaut@univ-lyon1.fr>",
)

if DEPLOYDOCS
    deploydocs(;
               repo="github.com/emmt/PointSpreadFunctions.jl",
               )
end


text:
import FireflyModule.Firefly, FireflyModule.copy

function generate_mutated_fireflies(firefly::Firefly, distance_matrix, r::Int, k::Int)::Array{Firefly}
    """Generates mutated fireflies using inversion_mutation function.
    Returns the new generated firefly array.
    """
    return [inversion_mutation(copy(firefly), distance_matrix, r) for _ in 1:k]
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
using Random

include("usefulFunctions.jl");

intervalExperiment = nothing;
# intervalExperiment = (1, 4000);

fileNames = datasetNames();
for fileName in fileNames
    for strategy in 1:4
        fid = open(string("script_", fileName,"_strategy_",strategy),"w");
        write(fid,"#!/bin/bash\n");
        write(fid,"#SBATCH --job-name=", fileName, "_strategy_", string(strategy), "\n");
        if intervalExperiment!=nothing
            write(fid,"#SBATCH --array=", string(intervalExperiment[1]), "-", string(intervalExperiment[2]), "\n");
        end;

        write(fid,"#SBATCH --time=4-04:00:00\n");
        write(fid,"#SBATCH -p shared\n");
        write(fid,"#SBATCH --qos shared\n");
        write(fid,"#SBATCH -C cl2s\n");

        write(fid,"\n");
        write(fid,"module load julia/1.4.2\n");
        write(fid,"julia experimentsDoME.jl ", fileName, " ", string(strategy), " \$SLURM_ARRAY_TASK_ID\n");
        close(fid);
    end;
end;

fid = open("executeAll","w");
for numFileName in 1:length(fileNames)
    fileName = fileNames[numFileName];
    for strategy in 1:4
        write(fid,"sbatch script_", fileName, "_strategy_", string(strategy), "\n");
    end;
end;
close(fid);

fid = open("cancelAll","w");
for fileName in fileNames
    for strategy in 1:4
        write(fid,"scancel -n ", fileName, "_strategy_", string(strategy), "\n");
    end;
end;
close(fid);


text:
using ForeignCallbacks
using Test
using InteractiveUtils

struct Message
    id::Int
    data::Ptr{Cvoid}
end

struct Message2
    id::Int
    data::Ref{Int}
end

@testset "Constructor" begin
    @test_throws AssertionError ForeignCallbacks.LockfreeQueue{Ref{Int}}()
    @test_throws AssertionError ForeignCallbacks.LockfreeQueue{Base.RefValue{Int}}()
    @test_throws AssertionError ForeignCallbacks.LockfreeQueue{Message2}()
    @test_throws AssertionError ForeignCallbacks.LockfreeQueue{Array}()
    @test_throws AssertionError ForeignCallbacks.LockfreeQueue{Array{Int, 1}}()
end

@testset "Queue" begin
    lfq = ForeignCallbacks.LockfreeQueue{Int}()

    @test ForeignCallbacks.dequeue!(lfq) === nothing
    ForeignCallbacks.enqueue!(lfq, 1)
    @test ForeignCallbacks.dequeue!(lfq) === Some(1)
    @test ForeignCallbacks.dequeue!(lfq) === nothing

    GC.@preserve lfq begin
        ptr = Base.pointer_from_objref(lfq)
        ForeignCallbacks.unsafe_enqueue!(ptr, 2)
    end
    # TODO: Test no load from TLS in `unsafe_enqueue!`
    @test ForeignCallbacks.dequeue!(lfq) === Some(2)
    @test ForeignCallbacks.dequeue!(lfq) === nothing
end

@testset "callback" begin
    ch = Channel{Int}(Inf)
    callback = ForeignCallbacks.ForeignCallback{Int}() do val
        put!(ch, val)
        return
    end

    GC.@preserve callback begin
        token = ForeignCallbacks.ForeignToken(callback)
        ForeignCallbacks.notify!(token, 1)
        @test fetch(ch) === 1
    end
end

@testset "IR" begin
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.enqueue!, Tuple{ForeignCallbacks.LockfreeQueue{Int}, Int}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.unsafe_enqueue!, Tuple{Ptr{Cvoid}, Int}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.notify!, Tuple{ForeignCallbacks.ForeignToken, Int}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.unsafe_enqueue!, Tuple{Ptr{Cvoid}, Message}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.notify!, Tuple{ForeignCallbacks.ForeignToken, Message}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
end

if Threads.nthreads() == 1 && Sys.CPU_THREADS > 1
    @info "relaunching with" threads = Sys.CPU_THREADS
    cmd = `$(Base.julia_cmd()) --threads=$(Sys.CPU_THREADS) $(@__FILE__)`
    @test success(pipeline(cmd, stdout=stdout, stderr=stderr))
    exit()
end

function producer!(lfq)
    for i in 1:100
        ForeignCallbacks.enqueue!(lfq, i)
        yield()
    end
end

function unsafe_producer!(lfq)
    for i in 1:100
        GC.@preserve lfq begin
            ptr = Base.pointer_from_objref(lfq)
            ForeignCallbacks.unsafe_enqueue!(ptr, i)
        end
        yield()
    end
end

function consumer!(lfq)
    acc = 0

    done = false
    while !done
        data = ForeignCallbacks.dequeue!(lfq)
        if data !== nothing
            acc += something(data)
        end
        done = acc == sum(1:100)*2*Threads.nthreads()
        yield()
    end
end

@testset "Queue threads" begin
    @test Threads.nthreads() == Sys.CPU_THREADS
    let lfq = ForeignCallbacks.LockfreeQueue{Int}()
        @sync begin
            for n in 1:2*Threads.nthreads()
                Threads.@spawn producer!(lfq)
            end
            Threads.@spawn consumer!(lfq)
        end
        @test true
    end

    let lfq = ForeignCallbacks.LockfreeQueue{Int}()
        @sync begin
            for n in 1:2*Threads.nthreads()
                Threads.@spawn unsafe_producer!(lfq)
            end
            Threads.@spawn consumer!(lfq)
        end
        @test true
    end
end

@testset "Callback threads" begin
    ch = Channel{Int}(Inf)
    callback = ForeignCallbacks.ForeignCallback{Int}() do val
        put!(ch, val)
        return
    end

    GC.@preserve callback begin
        token = ForeignCallbacks.ForeignToken(callback)
        for n in 1:2*Threads.nthreads()
            Threads.@spawn ForeignCallbacks.notify!(token, 1)
        end
    end

    acc = 0
    while acc < 2*Threads.nthreads()
        acc += fetch(ch)
    end
    @test acc == 2*Threads.nthreads()
end


text:
using ITensors
using ITensorGLMakie

include(joinpath("..", "src", "visualize_dangling.jl"))

i = Index(2, "i")
j = Index(2, "j")
k = Index(2, "k")

V = ITensor(i)
M = ITensor(i, j)
T = ITensor(i, j, k)

fig_V = visualize_dangling(V; vertex_labels="V")
fig_M = visualize_dangling(M; vertex_labels="M")
fig_T = visualize_dangling(T; vertex_labels="T")

fig_MV = visualize_dangling([M, V]; vertex_labels=["M", "V"])
fig_MM = visualize_dangling([M, M]; vertex_labels=["M", "M"])
fig_TV = visualize_dangling([T, V]; vertex_labels=["T", "V"])
fig_TM = visualize_dangling([T, M]; vertex_labels=["T", "M"])

# using Makie
# save("figure.png", fig)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    KrigingEstimator

A Kriging estimator (e.g. Simple Kriging).
"""
abstract type KrigingEstimator end

"""
    KrigingState(data, LHS, RHS, VARTYPE)

A Kriging state stores information needed
to perform estimation at any given location.
"""
mutable struct KrigingState{D<:Data,F<:Factorization,T,V}
  data::D
  LHS::F
  RHS::Vector{T}
  VARTYPE::V
end

"""
    KrigingWeights(λ, ν)

An object storing Kriging weights `λ` and Lagrange multipliers `ν`.
"""
struct KrigingWeights{T<:Real,A<:AbstractVector{T}}
  λ::A
  ν::A
end

"""
    FittedKriging(estimator, state)

An object that can be used for making predictions using the
parameters in `estimator` and the current Kriging `state`.
"""
struct FittedKriging{E<:KrigingEstimator,S<:KrigingState}
  estimator::E
  state::S
end

"""
    status(fittedkrig)

Return the status of the `fittedkrig` object, meaning
the factorization of the Kriging system was successful.
"""
status(fittedkrig::FittedKriging) = issuccess(fittedkrig.state.LHS)

#--------------
# FITTING STEP
#--------------

"""
    fit(estimator, data)

Build Kriging system from `data` and return a fitted estimator.
"""
function fit(estimator::KrigingEstimator, data)
  # variogram and domain
  γ = estimator.γ
  D = domain(data)

  # build Kriging system
  LHS = lhs(estimator, D)
  RHS = Vector{eltype(LHS)}(undef, size(LHS,1))

  # factorize LHS
  FLHS = factorize(estimator, LHS)

  # variance type
  VARTYPE = Variography.result_type(γ, first(D), first(D))

  # record Kriging state
  state = KrigingState(data, FLHS, RHS, VARTYPE)

  # return fitted estimator
  FittedKriging(estimator, state)
end

"""
    lhs(estimator, domain)

Return LHS of Kriging system for the elements in the `domain`.
"""
function lhs(estimator::KrigingEstimator, domain)
  γ = estimator.γ
  nobs = nelements(domain)
  ncon = nconstraints(estimator)

  # pre-allocate memory for LHS
  u  = first(domain)
  V² = Variography.result_type(γ, u, u)
  m  = nobs + ncon
  G  = Matrix{V²}(undef, m, m)

  # set variogram/covariance block
  pairwise!(G, γ, domain)
  if isstationary(γ)
    σ² = sill(γ)
    for j in 1:nobs, i in 1:nobs
      @inbounds G[i,j] = σ² - G[i,j]
    end
  end

  # strip units if necessary
  LHS = ustrip.(G)

  # set blocks of constraints
  set_constraints_lhs!(estimator, LHS, domain)

  LHS
end

"""
    nconstraints(estimator)

Return number of constraints for `estimator`.
"""
function nconstraints end

"""
    set_constraints_lhs!(estimator, LHS, X)

Set constraints in LHS of Kriging system.
"""
function set_constraints_lhs! end

"""
    factorize(estimator, LHS)

Factorize LHS of Kriging system with appropriate factorization method.
"""
function factorize end

#-----------------
# PREDICTION STEP
#-----------------

"""
    predict(estimator, var, uₒ)

Compute mean and variance of variable `var` using the
`estimator` at point or geometry `uₒ`.
"""
function predict(fitted::FittedKriging, var, uₒ)
  data = fitted.state.data
  combine(fitted, weights(fitted, uₒ), data[var])
end

"""
    weights(estimator, uₒ)

Compute the weights λ (and Lagrange multipliers ν) for the
`estimator` at point or geometry `uₒ`.
"""
function weights(fitted::FittedKriging, uₒ)
  nobs = nelements(fitted.state.data)

  set_rhs!(fitted, uₒ)

  # solve Kriging system
  s = fitted.state.LHS \ fitted.state.RHS

  λ = view(s, 1:nobs)
  ν = view(s, nobs+1:length(s))

  KrigingWeights(λ, ν)
end

"""
    set_rhs!(estimator, uₒ)

Set RHS of Kriging system at point or geometry `uₒ`.
"""
function set_rhs!(fitted::FittedKriging, uₒ)
  γ = fitted.estimator.γ
  dom = domain(fitted.state.data)
  nel = nelements(dom)
  RHS = fitted.state.RHS

  # RHS variogram/covariance
  g = map(u -> γ(u, uₒ), dom)
  RHS[1:nel] .= ustrip.(g)
  if isstationary(γ)
    σ² = ustrip(sill(γ))
    RHS[1:nel] .= σ² .- RHS[1:nel]
  end

  set_constraints_rhs!(fitted, uₒ)
end

"""
    set_constraints_rhs!(estimator, xₒ)

Set constraints in RHS of Kriging system.
"""
function set_constraints_rhs! end

"""
    combine(estimator, weights, z)

Combine `weights` with values `z` to produce mean and variance
using the appropriate formulas for the `estimator`.
"""
function combine(fitted::FittedKriging,
                 weights::KrigingWeights,
                 z::AbstractVector)
  γ  = fitted.estimator.γ
  b  = fitted.state.RHS
  V² = fitted.state.VARTYPE
  λ  = weights.λ
  ν  = weights.ν

  # compute b⋅[λ;ν]
  nobs = length(λ)
  c₁ = view(b, 1:nobs) ⋅ λ
  c₂ = view(b, nobs+1:length(b)) ⋅ ν
  c = c₁ + c₂

  if isstationary(γ)
    σ² = sill(γ)
    sum(λ.*z), σ² - V²(c)
  else
    sum(λ.*z), V²(c)
  end
end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("estimators/sk.jl")
include("estimators/ok.jl")
include("estimators/uk.jl")
include("estimators/edk.jl")


text:
# Class weight : w_j = n_samples / (n_classes * n_samples_j)

"""
    metrics(cm)
    metrics(cm, met...)
* `:TPR` true positive rate
* `:TNR` true negative rate
* `:NPV` negative predictive value
* `:FNR` false negative rate
* `:FPR` false positive rate
* `:FDR` false discovery rate
* `:FOR` false omission rate
* `:ACC` accuracy
* `(:Fscore, β)` Fβ score
* `:PT` prevalence threshold
"""
function metrics(cm::ConfMat{2, T}) where T
    TP = cm.counts[1, 1] |> float
    TN = cm.counts[2, 2] |> float
    FP = cm.counts[2, 1] |> float
    FN = cm.counts[1, 2] |> float
    Dict(:TP => TP, :TN => TN, :FP => FP, :FN => FN)
end

function metrics(cm::ConfMat)
    TP, TN, FP, FN = [], [], [], []
    for l ∈ cm.labels
        push!(TP, cm.counts[l, l] .|> float)
        push!(TN, sum(cm.counts[Not(l), Not(l)]) .|> float)
        push!(FP, sum(cm.counts[Not(l), l]) .|> float)
        push!(FN, sum(cm.counts[l, Not(l)]) .|> float)
    end
    Dict(:TP => TP, :TN => TN, :FP => FP, :FN => FN)
end

function metrics(cm::ConfMat, met...)
    metrics(cm) |> m -> metrics!(m, met...)
end

tosymbol(x) = x
function tosymbol(x::Tuple{Symbol, Real})
    s, v = x
    st = string(s) * "_" * string(v)
    st = replace(st, "." => "_")
    Symbol(st)
end

function metrics!(mm::Dict, met...)
    for m ∈ met
        if !haskey(mm, m)
            # F-score special case handling : tuple is transformed to a symbol
            k = tosymbol(m)
            mm[k] = metric(mm, m)
        end
    end
    mm
end

@inline nozero(x) = iszero(float(x)) ? eps(float(x)) : x

#avoidzero2(x::Int) = (avoidzero2(x) |> float)

# function safediv(n, d)
#     d = float(d)
#     float(n) / (iszero(d) ? eps(d) : d)
# end

# macro safediv(d)
#     :(iszero.($d) ? eps.(float.($d)) : float.($d))
# end

# vérif les formules

metric(b, m::Symbol) = metric(b, Val{m}())

metric(b, m::Val{:TPR}) = @. b[:TP] / nozero(b[:TP] + b[:FN])

metric(b, m::Val{:TNR}) = @. b[:TN] / nozero(b[:TN] + b[:FP])

metric(b, m::Val{:NPV}) = @. b[:TN] / nozero(b[:TN] + b[:FN])

metric(b, m::Val{:FNR}) = @. b[:FN] / nozero(b[:FN] + b[:TP])

metric(b, m::Val{:FPR}) = @. b[:FP] / nozero(b[:FP] + b[:TN])

metric(b, m::Val{:FDR}) = @. b[:FP] / nozero(b[:FP] + b[:TP])

metric(b, m::Val{:FOR}) = @. b[:FN] / nozero(b[:FN] + b[:TN])

metric(b, m::Val{:ACC}) =
    @. (b[:TP] + b[:TN]) / nozero(b[:TP] + b[:TN] + b[:FP] + b[:FN])

#_metric(b, t::Int) = println("pouet")

metric(b, m::Symbol, p) = metric(b, Val{m}(), p)

metric(b, t::Tuple{Symbol, Real}) = metric(b, t...)

metric(b, m::Val{:Fscore}, β) = begin
    a = (1 + β^2)
    @. a * b[:TP] / nozero(a * b[:TP] + β^2 * b[:FN] + b[:FP])
end

function metric(b, m::Val{:PT})
    fpr = haskey(b, :FPR) ? b[:FPR] : metric(b, :FPR)
    tpr = haskey(b, :TPR) ? b[:TPR] : metric(b, :TPR)
    a = @. sqrt(fpr)
    @. a / nozero(sqrt(tpr) + a)
end


text:
#
# This is the benefit/tax credit/IT/MinWage/NI rates from April 2022
# sources:
# https://www.gov.uk/government/publications/benefit-and-pension-rates-2022-to-2023
# https://www.gov.uk/government/publications/rates-and-allowances-tax-credits-child-benefit-and-guardians-allowance/tax-credits-child-benefit-and-guardians-allowance
# https://www.gov.uk/income-tax-rates
# https://www.gov.uk/guidance/rates-and-thresholds-for-employers-2022-to-2023
# https://www.gov.scot/publications/scottish-income-tax-2022-2023/
#
sys.it.savings_rates = [0.0, 20.0, 40.0, 45.0]
sys.it.savings_thresholds = [5_000.0, 37_700.0, 150_000.0]
sys.it.savings_basic_rate = 2 # above this counts as higher rate

sys.it.non_savings_rates = [19.0,20.0,21.0,41.0,46.0]
sys.it.non_savings_thresholds = [2_162, 13_118, 31_092, 150_000.0]
sys.it.non_savings_basic_rate = 2 # above this counts as higher rate rate FIXME 3???

sys.it.dividend_rates = [0.0, 8.75,33.75,39.35]
sys.it.dividend_thresholds = [2_000.0, 37_700.0, 150_000.0] # FIXME this gets the right answers & follows Melville, but the 2k is called 'dividend allowance in HMRC docs'
sys.it.dividend_basic_rate = 2 # above this counts as higher

sys.it.personal_allowance   = 12_570.00
sys.it.personal_allowance_income_limit = 100_000.00
sys.it.personal_allowance_withdrawal_rate = 50.0
sys.it.blind_persons_allowance  = 2_600.00

sys.it.married_couples_allowance = 9_415.00
sys.it.mca_minimum     = 3_640.00
sys.it.mca_income_maximum   = 31_400.00
sys.it.mca_credit_rate    = 10.0
sys.it.mca_withdrawal_rate  = 50.0

sys.it.marriage_allowance   = 1_260.00
sys.it.personal_savings_allowance = 1_000.00

# 29/2022 CHECK
# FIXME better to have it straight from
# the book with charges per CO2 range
# and the data being an estimate of CO2 per type
merge( sys.it.company_car_charge_by_CO2_emissions,
 Dict([
  Missing_Fuel_Type=>0.1,
  No_Fuel=>0.1,
  Other=>0.1,
  Dont_know=>0.1,
  Petrol=>0.25,
  Diesel=>0.37,
  Hybrid_use_a_combination_of_petrol_and_electricity=>0.16,
  Electric=>0.02,
  LPG=>0.02,
  Biofuel_eg_E85_fuel=>0.02 ]))

sys.it.fuel_imputation = 25_300.00 # 22/23

#
# pensions !!! CHECKUNCHANED 2022
#
sys.it.pension_contrib_basic_amount = 3_600.00
sys.it.pension_contrib_annual_allowance = 40_000.00
sys.it.pension_contrib_annual_minimum = 4_000.00
sys.it.pension_contrib_threshold_income = 240_000.00
sys.it.pension_contrib_withdrawal_rate = 50.0

# sys.it.non_savings_income = NON_SAVINGS_INCOME
# sys.it.all_taxable = ALL_TAXABLE_INCOME
# sys.it.savings_income = SAVINGS_INCOME
# sys.it.dividend_income = DIVIDEND_INCOME

# sys.it.mca_date = MCA_DATE

# NOTE THESE DO NOT INCLUDE JULY CHANGES

sys.ni.abolished = false
sys.ni.primary_class_1_rates = [0.0, 0.0, 13.25, 3.25 ]
sys.ni.primary_class_1_bands = [123.0, 190.0, 967.0, 9999999999999.9] # the '-1' here is because json can't write inf
sys.ni.secondary_class_1_rates = [0.0, 15.05, 15.05 ] # keep 2 so
sys.ni.secondary_class_1_bands = [175.0, 967.0, 99999999999999.9 ]
sys.ni.state_pension_age = 66; # fixme move
sys.ni.class_2_threshold = 6_725.0;
sys.ni.class_2_rate = 3.15;
sys.ni.class_4_rates = [0.0, 10.25, 3.25 ]
sys.ni.class_4_bands = [9_880.0, 50_270.0, 99999999999999.9 ]
# sys.ni.class_1_income = IncludedItems([WAGES],[PENSION_CONTRIBUTIONS_EMPLOYER])
# sys.ni.class_4_income = [SELF_EMPLOYMENT_INCOME]

# sys.uc.

sys.lmt.isa_jsa_esa_abolished = false
sys.lmt.pen_credit_abolished = false
## FIXME we can't turn off pension credit individually here..

sys.lmt.premia.family = 17.85
sys.lmt.premia.family_lone_parent = 22.20
sys.lmt.premia.carer_single = 38.85
sys.lmt.premia.carer_couple = 2*38.85
sys.lmt.premia.disabled_child = 68.04
sys.lmt.premia.disability_single = 36.20
sys.lmt.premia.disability_couple = 51.60
sys.lmt.premia.enhanced_disability_child = 26.67
sys.lmt.premia.enhanced_disability_single = 17.75
sys.lmt.premia.enhanced_disability_couple = 25.35
sys.lmt.premia.severe_disability_single = 69.40
sys.lmt.premia.severe_disability_couple = 138.80
sys.lmt.premia.pensioner_is = 157.65

sys.lmt.allowances.age_18_24 = 61.05
sys.lmt.allowances.age_25_and_over = 77.00
sys.lmt.allowances.age_18_and_in_work_activity = 77.00
sys.lmt.allowances.over_pension_age = 197.10 #
sys.lmt.allowances.lone_parent = 77.00
sys.lmt.allowances.lone_parent_over_pension_age = 197.10
sys.lmt.allowances.couple_both_under_18 = 61.05
sys.lmt.allowances.couple_both_over_18 = 121.05
sys.lmt.allowances.couple_over_pension_age = 294.90
sys.lmt.allowances.couple_one_over_18_high = 121.05
sys.lmt.allowances.couple_one_over_18_med = 77.00
sys.lmt.allowances.pa_couple_one_over_18_low = 61.05
sys.lmt.allowances.child = 70.80
sys.lmt.allowances.pc_mig_single = 182.60
sys.lmt.allowances.pc_mig_couple = 278.70
sys.lmt.allowances.pc_child = 56.35

# sys.lmt.income_rules.
sys.lmt.income_rules.permitted_work =152.00
sys.lmt.income_rules.lone_parent_hb = 25.0
sys.lmt.income_rules.high = 20.0
sys.lmt.income_rules.low_couple = 10.0
sys.lmt.income_rules.low_single = 5.0
sys.lmt.income_rules.hb_additional = 17.10
sys.lmt.income_rules.childcare_max_1 = 175.0
sys.lmt.income_rules.childcare_max_2 = 300.0
sys.lmt.income_rules.incomes  = LEGACY_MT_INCOME
sys.lmt.income_rules.hb_incomes = LEGACY_HB_INCOME
sys.lmt.income_rules.pc_incomes = LEGACY_PC_INCOME
sys.lmt.income_rules.sc_incomes = LEGACY_SAVINGS_CREDIT_INCOME
sys.lmt.income_rules.capital_min = 6_000.0
sys.lmt.income_rules.capital_max = 16_000.0
sys.lmt.income_rules.pc_capital_min = 10_000.0
sys.lmt.income_rules.pc_capital_max = 99999999999999.9
sys.lmt.income_rules.pensioner_capital_min = 10_000.0
sys.lmt.income_rules.pensioner_capital_max = 16_000.0

sys.lmt.income_rules.capital_tariff = 250
sys.lmt.income_rules.pensioner_tariff = 500
# FIXME why do we need a seperate copy of HoursLimits here?

sys.lmt.hours_limits.lower = 16
sys.lmt.hours_limits.med = 24
sys.lmt.hours_limits.higher = 30

sys.lmt.savings_credit.abolished = false
sys.lmt.savings_credit.withdrawal_rate = 60.0
sys.lmt.savings_credit.threshold_single = 158.47
sys.lmt.savings_credit.threshold_couple = 251.70
sys.lmt.savings_credit.max_single = 14.48
sys.lmt.savings_credit.max_couple = 16.20
sys.lmt.savings_credit.available_till = Date( 2016, 04, 06 )

sys.lmt.child_tax_credit.abolished = false
sys.lmt.child_tax_credit.family = 545.0
sys.lmt.child_tax_credit.child = 2_935.0
sys.lmt.child_tax_credit.disability = 3_545.0
sys.lmt.child_tax_credit.severe_disability = 4825
sys.lmt.child_tax_credit.threshold = 16_480.0

sys.lmt.working_tax_credit.basic = 2_070.00
sys.lmt.working_tax_credit.lone_parent = 2_125.00
sys.lmt.working_tax_credit.couple = 2_125.00
sys.lmt.working_tax_credit.hours_ge_30 = 860.00
sys.lmt.working_tax_credit.disability = 3_345.00
sys.lmt.working_tax_credit.severe_disability = 1_445.00
sys.lmt.working_tax_credit.age_50_plus = 1_365.00 # ??
sys.lmt.working_tax_credit.age_50_plus_30_hrs = 2_030.00
sys.lmt.working_tax_credit.childcare_max_2_plus_children = 300.0 # pw
sys.lmt.working_tax_credit.childcare_max_1_child = 175.0
sys.lmt.working_tax_credit.childcare_proportion = 70.0
sys.lmt.working_tax_credit.taper = 41.0
sys.lmt.working_tax_credit.threshold = 6_770.0
sys.lmt.working_tax_credit.non_earnings_minima = 300.0 # FIXME check

sys.lmt.hb.taper = 65.0
sys.lmt.hb.ndd_deductions = [16.45,37.80,51.85,84.85,96.60,106.05]
sys.lmt.hb.ndd_incomes = [154.0,224.0,292.0,389.0,484.0,99999999999999.9]

sys.lmt.ctr.taper = 20.0
sys.lmt.ctr.ndd_deductions = []
sys.lmt.ctr.ndd_incomes = []

sys.uc.abolished = false
sys.uc.threshold = 2_500.0 ## NOT USED
sys.uc.age_18_24 = 265.31
sys.uc.age_25_and_over = 334.91

sys.uc.couple_both_under_25 = 416.45
sys.uc.couple_oldest_25_plus = 525.72

sys.uc.first_child = 290.00
sys.uc.subsequent_child = 244.58
sys.uc.disabled_child_lower = 132.89
sys.uc.disabled_child_higher = 414.88
sys.uc.limited_capcacity_for_work_activity = 354.28
sys.uc.carer = 168.81

sys.uc.ndd = 77.87

sys.uc.childcare_max_2_plus_children = 1_108.04 # pm
sys.uc.childcare_max_1_child = 646.35
sys.uc.childcare_proportion = 85.0 # pct

sys.uc.minimum_income_floor_hours = 35*WEEKS_PER_MONTH

sys.uc.work_allowance_w_housing = 344.00
sys.uc.work_allowance_no_housing = 573.00
sys.uc.other_income = UC_OTHER_INCOME
# sys.uc.earned_income :: IncludedItems = UC_EARNED_INCOME
sys.uc.capital_min = 6_000.0
sys.uc.capital_max = 16_000.0
# £1 *per week* ≆ 4.35 pm FIXME make 4.35 WEEKS_PER_MONTH?
sys.uc.capital_tariff = 250.0/4.35
sys.uc.taper = 63.0
sys.uc.ctr_taper = 20.0

sys.age_limits.state_pension_ages = pension_ages()
sys.age_limits.savings_credit_to_new_state_pension = Date( 2016, 04, 06 )

sys.hours_limits.lower = 16
sys.hours_limits.med = 24
sys.hours_limits.higher = 30

sys.child_limits.max_children = 2

# THESE ARE TILL APRIL 2022
sys.minwage.ages = [16,18,21,23]
sys.minwage.wage_per_hour = [4.62, 6.56, 8.36, 8.91]
sys.minwage.apprentice_rate = 4.30

sys.hr.maximum_rooms = 4
sys.hr.rooms_rent_reduction = [14.0,25.0]
sys.hr.single_room_age = 35
#
# These are unchanged in 3 years; see:
# https://www.gov.scot/publications/local-housing-allowance-rates-2021-2022/
#
sys.hr.brmas = loadBRMAs( 4, T, DEFAULT_BRMA_2021 )

sys.nmt_bens.attendance_allowance.abolished = false
sys.nmt_bens.attendance_allowance.higher = 92.40
sys.nmt_bens.attendance_allowance.lower = 61.85

sys.nmt_bens.child_benefit.abolished = false
sys.nmt_bens.child_benefit.first_child = 21.80
sys.nmt_bens.child_benefit.other_children = 14.45
sys.nmt_bens.child_benefit.high_income_thresh = 50_000.0
sys.nmt_bens.child_benefit.withdrawal = 1/100
sys.nmt_bens.child_benefit.guardians_allowance = 18.55

sys.nmt_bens.dla.abolished = false
sys.nmt_bens.dla.care_high = 92.40
sys.nmt_bens.dla.care_middle = 61.85
sys.nmt_bens.dla.care_low = 24.45
sys.nmt_bens.dla.mob_high = 64.50
sys.nmt_bens.dla.mob_low = 24.45

sys.nmt_bens.carers.abolished = false
sys.nmt_bens.carers.allowance = 69.70

sys.nmt_bens.carers.scottish_supplement = 231.40 # FROM APRIL 2021
sys.nmt_bens.carers.hours :: Int = 35
sys.nmt_bens.carers.gainful_employment_min = 132.00

sys.nmt_bens.pip.abolished = false
sys.nmt_bens.pip.dl_standard = 61.85
sys.nmt_bens.pip.dl_enhanced = 92.40
sys.nmt_bens.pip.mobility_standard = 24.45
sys.nmt_bens.pip.mobility_enhanced = 64.50

sys.nmt_bens.esa.abolished = false
sys.nmt_bens.esa.assessment_u25 = 61.05
sys.nmt_bens.esa.assessment_25p = 77.00
sys.nmt_bens.esa.main   = 77.00
sys.nmt_bens.esa.work   = 30.60
sys.nmt_bens.esa.support  = 40.60

sys.nmt_bens.jsa.abolished = false
sys.nmt_bens.jsa.u25 = 61.05
sys.nmt_bens.jsa.o24 = 77.00

sys.nmt_bens.pensions.abolished = false
sys.nmt_bens.pensions.new_state_pension = 185.15
# pension_start_date = Date( 2016, 04, 06 )
sys.nmt_bens.pensions.cat_a  = 141.85
sys.nmt_bens.pensions.cat_b  = 141.85
sys.nmt_bens.pensions.cat_b_survivor = 85.00
sys.nmt_bens.pensions.cat_d  = 85.00

sys.nmt_bens.bereavement.abolished = false
# higher effectively just means 'with children';
sys.nmt_bens.bereavement.lump_sum_higher = 3_500
sys.nmt_bens.bereavement.lump_sum_lower = 2_500
sys.nmt_bens.bereavement.higher = 350
sys.nmt_bens.bereavement.lower = 100

sys.nmt_bens.widows_pension.abolished = false
sys.nmt_bens.widows_pension.industrial_higher = 141.85
sys.nmt_bens.widows_pension.industrial_lower = 42.56
sys.nmt_bens.widows_pension.standard_rate = 126.35
sys.nmt_bens.widows_pension.parent = 126.35
sys.nmt_bens.widows_pension.ages = collect(54:-1:45)
sys.nmt_bens.widows_pension.age_amounts = [117.51,108.66,99.82,90.97,82.13,73.28,64.44,55.59,46.75,37.91]

#
# young carer grant
sys.nmt_bens.maternity.abolished = false
sys.nmt_bens.maternity.rate = 156.66

sys.nmt_bens.smp = 156.66 ## 90% of earn cpag 21/2 812
# = XX

sys.bencap.abolished = false
sys.bencap.outside_london_single = 257.69
sys.bencap.outside_london_couple = 384.62
# not really needed, but anyway ..
sys.bencap.inside_london_single = 296.35
sys.bencap.inside_london_couple = 442.31
sys.bencap.uc_incomes_limit = 617

sys.scottish_child_payment.amount = 10.0
sys.scottish_child_payment.maximum_age = 5

sys.ubi.abolished = true
sys.ubi.adult_amount = 4_800.0
sys.ubi.child_amount= 3_000.0
sys.ubi.universal_pension = 8_780.0
sys.ubi.adult_age = 17
sys.ubi.retirement_age = 66

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
export HyperRect

"""
    HyperRect{N, Shape, BC} <: BoundedLattice{N, BC}

Hyper Rectangle lattice on `N` dimension with size `Shape`.
"""
struct HyperRect{N, BC, Shape <: Tuple} <: BoundedLattice{N, BC}
end

HyperRect(shape::Int...; boundary=Periodic) = HyperRect(shape; boundary=boundary)
HyperRect(shape::NTuple{N, Int}; boundary=Periodic) where N = HyperRect{N, boundary, Tuple{shape...}}()
nameof(::HyperRect{N}) where N = "HyperRect{$N}"
Base.@pure Base.size(::HyperRect{N, BC, Shape}) where {N, BC, Shape} = tuple(Shape.parameters...)

function sites(ltc::HyperRect)
    Base.CartesianIndices((Base.OneTo(k + 1) for k in size(ltc))...)
end

function edges(ltc::HyperRect; order=1)
    error("Not Implemented Yet")
end

function neighbors(::HyperRect{N, Periodic, Shape}, P::CartesianIndex; order=1) where {N, Shape}
    @boundscheck all(x->x>0, Tuple(P)) || error("Lattice configuration position starts from 1")
    hyper_rect_neighbors(P, Shape, Val(length))
end

# NOTE: we assume HyperRect lattice only work for low dimension case, therefore
#       we will unroll the neighbors iterator to a tuple directly.
@generated function hyper_rect_neighbors(P::CartesianIndex{D}, L::Type{Shape}, ::Val{len}) where D where {Shape <: Tuple, len}
    ex = Expr(:tuple)

    for i in 1:D
        push!(ex.args, gen_cartesian_nn_index_expr(:P, D, i, Shape.parameters[i], len))
        push!(ex.args, gen_cartesian_nn_index_expr(:P, D, i, Shape.parameters[i], -len))
    end
    ex
end

function gen_cartesian_nn_index_expr(P::Symbol, D::Int, i::Int, L::Int, k::Int)
    ind_ex = Expr(:call, :CartesianIndex)
    for index in 1:D
        if index == i
            push!(ind_ex.args, :(mod1($P[$index] + $k, $(L + 1))))
        else
            push!(ind_ex.args, :($P[$index]))
        end
    end
    ind_ex
end