text:
function (+)(x::SimpleQuaternion, y::SimpleQuaternion)
    cz = get_parts(x) .+ get_parts(y)
    return SimpleQuaternion(cz...)
end

(+)(x::SimpleQuaternion, y::Number) = x + SimpleQuaternion(y)
(+)(x::Number, y::SimpleQuaternion) = SimpleQuaternion(x) + y

function (-)(x::SimpleQuaternion)
    a,b,c,d = get_parts(x)
    return SimpleQuaternion(-a,-b,-c,-d)
end

function (-)(x::SimpleQuaternion,y::SimpleQuaternion)
    cz = get_parts(x) .- get_parts(y)
    return SimpleQuaternion(cz...)
end

(-)(x::SimpleQuaternion,y::Number) = x - SimpleQuaternion(y)
(-)(x::Number, y::SimpleQuaternion) = SimpleQuaternion(x) - y

function (*)(x::SimpleQuaternion, y::SimpleQuaternion)
    a,b,c,d = get_parts(x)
    aa,bb,cc,dd, = get_parts(y)

    w = a*aa - b*bb - c*cc - d*dd
    x = a*bb + b*aa + c*dd - d*cc
    y = a*cc - b*dd + c*aa + d*bb
    z = a*dd + b*cc - c*bb + d*aa

    return SimpleQuaternion(w,x,y,z)
end

(*)(x::SimpleQuaternion, y::Number) = x * SimpleQuaternion(y)
(*)(x::Number, y::SimpleQuaternion) = SimpleQuaternion(x) * y

function norm(x::SimpleQuaternion)
    a,b,c,d = get_parts(x)
    return sqrt(a*a + b*b + c*c + d*d)
end

function (inv)(x::SimpleQuaternion)
    a,b,c,d = get_parts(x)
    n = a*a + b*b + c*c + d*d
    return SimpleQuaternion(a/n, -b/n, -c/n, -d/n)
end

function (_inv)(x::SimpleQuaternion)
    a,b,c,d = get_parts(x)
    n = a*a + b*b + c*c + d*d
    return SimpleQuaternion(a//n, -b//n, -c//n, -d//n)
end

(/)(x::SimpleQuaternion, y::SimpleQuaternion) =  x * inv(y)
(/)(x::SimpleQuaternion,y::Number) =  x * inv(y)
(/)(x::Number,y::SimpleQuaternion) =  x * inv(y)

(//)(x::SimpleQuaternion, y::SimpleQuaternion) = x * _inv(y)
(//)(x::SimpleQuaternion, y::Number) = x * (1//y)
(//)(x::Number,y::SimpleQuaternion) = SimpleQuaternion(x) // y

"""
`adjoint(x::SimpleQuaternion)` or `x'` returns the adjoint of `x`.

If `x` is `a + b*im + c*jm + d*km` then `x'` is `a - b*im - c*im - d*im`.
"""
adjoint(x::SimpleQuaternion) = conj(x)

"""
`conj(x::SimpleQuaternion)` returns the conjugate of `x`.

If `x` is `a + b*im + c*jm + d*km` then `conj(x)` is `a - b*im - c*im - d*im`.
"""
function conj(x::SimpleQuaternion)
    a,b,c,d = get_parts(x)
    return SimpleQuaternion(a,-b,-c,-d)
end


text:
const Entry = Dict{String, Any}

allkeys(entry::Entry) = all(k -> haskey(entry, k), ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"])

validbirth(year::Int) = 1920 ≤ year ≤ 2002
validissue(year::Int) = 2010 ≤ year ≤ 2020
validexpiration(year::Int) = 2020 ≤ year ≤ 2030

validheight((height, unit)) = if unit == "cm"
    150 ≤ height ≤ 193
elseif unit == "in"
    59 ≤ height ≤ 76
else
    false
end

validhair(color::AbstractString) = !isnothing(match(r"#[0-9a-f]{6}$", color))
valideye(color::AbstractString) = color in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
validpid(pid::AbstractString) = !isnothing(match(r"^\d{9}$", pid))

const Validator = Dict(
    "byr" => validbirth,
    "iyr" => validissue,
    "eyr" => validexpiration,
    "hgt" => validheight,
    "hcl" => validhair,
    "ecl" => valideye,
    "pid" => validpid
)

const Parser = Dict(
    "byr" => v -> parse(Int, v),
    "iyr" => v -> parse(Int, v),
    "eyr" => v -> parse(Int, v),
    "hgt" => v -> (parse(Int, v[1:end-2]), v[end-1:end]),
)

validate(entry::Entry) = all(entry) do (key, value)
    try
        !haskey(Validator, key) || Validator[key](value)
    catch
        false
    end
end

@enum Part Part1 Part2

isvalid(part::Part, entry) = if part == Part1
    allkeys(entry)
else
    isvalid(Part1, entry) && validate(entry)
end

let part1 = 0, part2 = 0
    entry = Entry()

    for line in readlines("input.txt")
        if isempty(line)
            part1 += isvalid(Part1, entry)
            part2 += isvalid(Part2, entry)
            empty!(entry)
        else
            fields = split(line, " ")
            for field in fields
                key, value = field[1:3], field[5:end]
                if haskey(Parser, key) &&
                    try
                        entry[key] = Parser[key](value)
                    catch
                    end
                else
                    entry[key] = value
                end
            end
        end
    end

    part1 += isvalid(Part1, entry)
    part2 += isvalid(Part2, entry)

    @info "Solutions" part1 part2
end

ext: jl

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
module projecttwo

include("mini.jl")

function line_optimize(f, a, b; eps=1e-3, method::SVOptMethod=mini())
	optimizer = method
	optimizer(f, a, b)
end

export line_optimize

end # module

max_stars_repo_name: madzias/Julia

max_issues_repo_name: madzias/Julia

max_forks_repo_name: madzias/Julia


text:
#=
dtft.jl
discrete-time Fourier transform
The main purpose of this code is for testing the NUFFT code.
2019-06-05, Jeff Fessler, University of Michigan
Note: the multi-dimensional version uses collect(Tuple(idx))
that might be inefficient, but for efficiency one should just use NUFFT.

see also MIRT/time/dtft.jl
=#

export dtft_init, dtft, dtft_adj

#using LinearAlgebra: norm
using Distributed: @sync, @distributed, pmap
using SharedArrays: SharedArray, sdata
using LinearMapsAA: LinearMapAA, LinearMapAO

"""
    X = dtft(w, x ; n_shift=?)
1D DTFT

``X[m] = \\sum_{n=0}^{N-1} x[n] \\exp(-i w[m] (n - n_{shift})), m=1,…,M``

in
- `w::AbstractVector{<:Real}`	`[M]` frequency locations ("units" radians/sample)
- `x::AbstractVector{<:Number}`	`[N]` 1D signal

option
- `n_shift::Real` often is N/2; default 0

out
- `X::AbstractVector{ComplexF64}` `[M]` DTFT
"""
function dtft(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	return dtft_loop_n(w, x ; n_shift=n_shift)
end

"""
    X = dtft(w, x ; n_shift=?)
multi-dimensional DTFT (DSFT)

``X[m] = \\sum_{n=0}^{N-1} x[n] \\exp(-i w[m,:] (n - n_{shift})), m=1,…,M``
where here `n` is a `CartesianIndex`

in
- `w::AbstractMatrix{<:Real}`	`[M,D]` frequency locations ("units" radians/sample)
- `x::AbstractArray{<:Number}`	`[(Nd)]` D-dimensional signal

option
- `n_shift::AbstractVector{<:Real}` often is N/2; default zeros(D)

out
- `X::AbstractVector{ComplexF64}` `[M]` DTFT
"""
function dtft(w::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Number}
	; n_shift::AbstractVector{<:Real} = zeros(Int, ndims(x)),
)
	return dtft_loop_n(w, x ; n_shift=n_shift)
end

"""
    d = dtft_init(w, N ; n_shift=?)
for 1D DTFT

in
- `w::AbstractVector{<:Real}`	`[M]` frequency locations ("units" radians/sample)
- `N::Int` size of signal `x`

option
- `n_shift::Real` often is N/2; default 0
- `T::DataType` default `ComplexF64` for testing NUFFT accuracy

out
- `d::NamedTuple` with fields
	`dtft = x -> dtft(x), adjoint = y -> dtft_adj(y), A=LinearMapAO`
"""
function dtft_init(
	w::AbstractVector{<:Real},
	N::Int ;
	n_shift::Real = 0,
	T::DataType = ComplexF64,
)

	M = length(w)
	forw = x -> dtft(w, x ; n_shift=n_shift)
	back = y -> dtft_adj(w, y, N ; n_shift=n_shift)
	A = LinearMapAA(forw, back, (M, N),
		(name="dtft1", N=(N,)) ; T = ComplexF64,
		operator = true,
	)
	return (dtft=forw, adjoint=back, A=A)
end

"""
    d = dtft_init(w, N ; n_shift=?)
for multi-dimensional DTFT (DSFT)

in
- `w::AbstractMatrix{<:Real}`	`[M,D]` frequency locations ("units" radians/sample)
- `N::Dims` `[D]` dimensions of signal `x`

option
- `n_shift::AbstractVector{<:Real}` often is N/2; default zeros(D)
- `T::DataType` default `ComplexF64` for testing NUFFT accuracy

out
- `d::NamedTuple` with fields
	`dtft = x -> dtft(x), adjoint = y -> dtft_adj(y), A=LinearMapAO`
"""
function dtft_init(
	w::AbstractMatrix{<:Real},
	N::Dims ;
	n_shift::AbstractVector{<:Real} = zeros(Int, length(N)),
)

	D = size(w,2)
	length(N) != D && throw(DimensionMismatch("length(N) vs D=$D"))
	length(n_shift) != D && throw(DimensionMismatch("length(n_shift) vs D=$D"))
	M = size(w,1)
	forw = x -> dtft(w, x ; n_shift=n_shift)
	back = y -> dtft_adj(w, y, N ; n_shift=n_shift)
	A = LinearMapAA(forw, back, (M, prod(N)),
		(name="dtft$(length(N))", N=N) ; T = ComplexF64,
		operator = true, idim = N,
	)
	return (dtft=forw, adjoint=back, A=A)
end

"""
    x = dtft_adj(w, X, N ; n_shift=?)
adjoint for 1D DTFT

``x[n] = \\sum_{m=1}^M X[m] \\exp(i w[m] (n - n_{shift})), n=0,…,N-1``

This is the *adjoint* (transpose) of `dtft`, not an *inverse* DTFT.

in
- `w::AbstractVector{<:Real}`	`[M]` frequency locations ("units" radians/sample)
- `X::AbstractVector{ComplexF64}`	`[M]` spectrum values
- `N::Int` size of signal `x`

option
- `n_shift::Real` often is N/2; default 0

out
- `x::AbstractVector{<:Number}`	signal [N]
"""
function dtft_adj(w::AbstractVector{<:Real}, X::AbstractVector{<:Number},
	N::Int ; n_shift::Real = 0,
)
	M = length(w)
	out = Array{ComplexF64}(undef, N)
	nshift1 = n_shift + 1
	X = ComplexF64.(X)
	for n=1:N # loop over output signal samples
		out[n] = sum(cis.((n - nshift1) * w) .* X)
	end
	return out
end

"""
    x = dtft_adj(w, X, N ; n_shift=?)
adjoint for multi-dimensional DTFT (DSFT)

``x[n] = \\sum_{m=1}^M X[m] \\exp(i w[m,:] (n - n_{shift})), n=0,…,N-1``
where here `n` is a `CartesianIndex`

in
- `X::AbstractVector{ComplexF64}` `[M]` DTFT
- `w::AbstractMatrix{<:Real}` `[M,D]` frequency locations ("units" radians/sample)
- `N::Dims`	`[D]` dimensions of signal `x`

option
- `n_shift::AbstractVector{<:Real}` often is `N/2`; default `zeros(D)`

out
- `x::AbstractArray{<:Number}` `[(N)]` `D`-dimensional signal
"""
function dtft_adj(
	w::AbstractMatrix{<:Real},
	X::AbstractVector{<:Number},
	N::Dims ;
	n_shift::AbstractVector{<:Real} = zeros(Int, ndims(x)),
)

	M,D = size(w)
	nshift1 = n_shift .+ 1

	T = (eltype(w) == Float64) ? ComplexF64 : ComplexF32
	out = Array{T}(undef, N)

	for n=1:prod(N)
		idx = CartesianIndices(N)[n]
		tmp = cis.(w * (collect(Tuple(idx)) - nshift1))
		tmp = transpose(cis.(w * (collect(Tuple(idx)) - nshift1))) * X
		out[idx] = cis.(-w * (collect(Tuple(idx)) - nshift1))' * X
	end
	return out
end

# 1D loop over signal values
function dtft_loop_n(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
	; n_shift::Real = 0,
)
	N = length(x)
	out = fill(ComplexF64(0), size(w))
	nshift1 = n_shift + 1
	for n=1:N # tried @simd but not faster
		out .+= x[n] * cis.(-(n-nshift1) * w)
	end
	return out
end

# dD loop over signal values
# * `w [M,D]`
# * `x [(Nd)]`
function dtft_loop_n(
	w::AbstractMatrix{<:Real},
	x::AbstractMatrix{<:Number} ;
	n_shift::AbstractVector{<:Real} = zeros(Int, size(w,2)),
)
	M,D = size(w)
	N = size(x)
	length(N) != D && throw(DimensionMismatch("length(N) vs D=$D"))
	length(n_shift) != D && throw(DimensionMismatch("length(n_shift) vs D=$D"))
	ndims(x) != D && throw(DimensionMismatch("ndims(x) vs D=$D"))
	nshift1 = n_shift .+ 1
	T = (eltype(w) == Float64) ? ComplexF64 : ComplexF32
	out = fill(T(0), M)
	for n=1:prod(N)
		idx = CartesianIndices(N)[n]
		out .+= x[idx] * cis.(-w * (collect(Tuple(idx)) - nshift1))
	end
	return out
end

# 1D loop over frequencies
function dtft_loop_m(
	w::AbstractVector{<:Real},
	x::AbstractVector{<:Number} ;
	n_shift::Real = 0,
)
	M = length(w)
	N = length(x)
	out = Array{ComplexF64}(undef, M)
	nn = -((0:(N-1)) .- n_shift)
	x = ComplexF64.(x)
	for m=1:M
		out[m] = sum(cis.(nn * w[m]) .* x)
	end
	return out
end

# 1D matrix * vector
function dtft_matvec(
	w::AbstractVector{<:Real},
	x::AbstractVector{<:Number} ;
	n_shift::Real = 0,
)
	N = length(x)
	return cis.(w * (-((0:(N-1)) .- n_shift))') * x
end

# 1D single frequency
function dtft_one_w(
	w::Real,
	x::AbstractVector{<:Number} ;
	n_shift::Real = 0,
)
	N = length(x)
	return sum(cis.(-((0:(N-1)) .- n_shift) * w) .* x)
end

# 1D pmap
function dtft_pmap_m(
	w::AbstractVector{<:Real},
	x::AbstractVector{<:Number}
	; n_shift::Real = 0,
)
	tmp = w -> dtft_one_w(w, x ; n_shift=n_shift)
	return pmap(tmp, w)
end

# 1D distributed loop over frequencies
function dtft_dist_m(
	w::AbstractVector{<:Real},
	x::AbstractVector{<:Number} ;
	n_shift::Real = 0,
)
	M = length(w)
	N = length(x)
	out = SharedArray{ComplexF64}(M)
	nn = -((0:(N-1)) .- n_shift)
	x = ComplexF64.(x)
	@sync @distributed for m=1:M
		out[m] = sum(cis.(nn * w[m]) .* x)
	end
	return sdata(out)
end


text:
struct ZIInverseGaussianExpert{T<:Real} <: ZIContinuousExpert
    p::T
    μ::T
    λ::T
    ZIInverseGaussianExpert{T}(p::T, µ::T, λ::T) where {T<:Real} = new{T}(p, µ, λ)
end

function ZIInverseGaussianExpert(p::T, μ::T, λ::T; check_args=true) where {T <: Real}
    check_args && @check_args(ZIInverseGaussianExpert, 0 <= p <= 1 && μ >= zero(μ) && λ > zero(λ))
    return ZIInverseGaussianExpert{T}(p, μ, λ)
end

#### Outer constructors
ZIInverseGaussianExpert(p::Real, μ::Real, λ::Real) = ZIInverseGaussianExpert(promote(p, μ, λ)...)
ZIInverseGaussianExpert(p::Integer, μ::Integer, λ::Integer) = ZIInverseGaussianExpert(float(p), float(μ), float(λ))
ZIInverseGaussianExpert() = ZIInverseGaussianExpert(0.5, 1.0, 1.0)

## Conversion
function convert(::Type{ZIInverseGaussianExpert{T}}, p::S, μ::S, λ::S) where {T <: Real, S <: Real}
    λ(T(p), T(μ), T(λ))
end
function convert(::Type{ZIInverseGaussianExpert{T}}, d::ZIInverseGaussianExpert{S}) where {T <: Real, S <: Real}
    ZIInverseGaussianExpert(T(d.p), T(d.μ), T(d.λ), check_args=false)
end
copy(d::ZIInverseGaussianExpert) = ZIInverseGaussianExpert(d.p, d.μ, d.λ, check_args=false)

## Loglikelihood of Expoert
logpdf(d::ZIInverseGaussianExpert, x...) = isinf(x...) ? -Inf : Distributions.logpdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)
pdf(d::ZIInverseGaussianExpert, x...) = isinf(x...) ? 0.0 : Distributions.pdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)
logcdf(d::ZIInverseGaussianExpert, x...) = isinf(x...) ? 0.0 : Distributions.logcdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)
cdf(d::ZIInverseGaussianExpert, x...) = isinf(x...) ? 1.0 : Distributions.cdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)

expert_ll_exact(d::ZIInverseGaussianExpert, x::Real) = (x == 0.) ? log(p_zero(d)) :  log(1-p_zero(d)) + LRMoE.logpdf(d, x)
function expert_ll(d::ZIInverseGaussianExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll_pos = LRMoE.expert_ll(LRMoE.InverseGaussianExpert(d.μ, d.λ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_ll = (yl == 0.) ? log.(p0 + (1-p0)*exp.(expert_ll_pos)) : log.(0.0 + (1-p0)*exp.(expert_ll_pos))
    expert_ll = (tu == 0.) ? log.(p0) : expert_ll
    return expert_ll
end
function expert_tn(d::ZIInverseGaussianExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_pos = LRMoE.expert_tn(LRMoE.InverseGaussianExpert(d.μ, d.λ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_tn = (tl == 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_pos))
    expert_tn = (tu == 0.) ? log.(p0) : expert_tn
    return expert_tn
end
function expert_tn_bar(d::ZIInverseGaussianExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar_pos = LRMoE.expert_tn_bar(LRMoE.InverseGaussianExpert(d.μ, d.λ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_tn_bar = (tl > 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_bar_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_bar_pos))
    return expert_tn_bar
end

exposurize_expert(d::ZIInverseGaussianExpert; exposure = 1) = d

## Parameters
params(d::ZIInverseGaussianExpert) = (d.p, d.μ, d.λ)
p_zero(d::ZIInverseGaussianExpert) = d.p
function params_init(y, d::ZIInverseGaussianExpert)
    p_init = sum(y .== 0.0) / sum(y .>= 0.0)
    pos_idx = (y .> 0.0)
    μ, σ2 = mean(y[pos_idx]), var(y[pos_idx])
    μ_init = μ
    λ_init = μ^3 / σ2
    if isnan(μ_init) || isnan(λ_init)
        return ZIInverseGaussianExpert()
    else
        return ZIInverseGaussianExpert(p_init, μ_init, λ_init)
    end
end

## KS stats for parameter initialization
function ks_distance(y, d::ZIInverseGaussianExpert)
    p_zero = sum(y .== 0.0) / sum(y .>= 0.0)
    return max(abs(p_zero-d.p), (1-d.p)*HypothesisTests.ksstats(y[y .> 0.0], Distributions.InverseGaussian(d.μ, d.λ))[2])
end

## Simululation
sim_expert(d::ZIInverseGaussianExpert) = (1 .- Distributions.rand(Distributions.Bernoulli(d.p), 1)[1]) .* Distributions.rand(Distributions.InverseGaussian(d.μ, d.λ), 1)[1]

## penalty
penalty_init(d::ZIInverseGaussianExpert) = [1.0 Inf 1.0 Inf]
no_penalty_init(d::ZIInverseGaussianExpert) = [1.0 Inf 1.0 Inf]
penalize(d::ZIInverseGaussianExpert, p) = (p[1]-1)*log(d.μ) - d.μ/p[2] + (p[3]-1)*log(d.λ) - d.λ/p[4]

## statistics
mean(d::ZIInverseGaussianExpert) = (1-d.p)*mean(Distributions.InverseGaussian(d.μ, d.λ))
var(d::ZIInverseGaussianExpert) = (1-d.p)*var(Distributions.InverseGaussian(d.μ, d.λ)) + d.p*(1-d.p)*(mean(Distributions.InverseGaussian(d.μ, d.λ)))^2
quantile(d::ZIInverseGaussianExpert, p) = p <= d.p ? 0.0 : quantile(Distributions.InverseGaussian(d.μ, d.λ), p-d.p)
lev(d::ZIInverseGaussianExpert, u) = (1-d.p)*lev(InverseGaussianExpert(d.μ, d.λ), u)
excess(d::ZIInverseGaussianExpert, u) = mean(d) - lev(d, u)

## EM: M-Step
function EM_M_expert(d::ZIInverseGaussianExpert,
                     tl, yl, yu, tu,
                     exposure,
                     z_e_obs, z_e_lat, k_e;
                     penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Update zero probability
    expert_ll_pos = expert_ll.(LRMoE.InverseGaussianExpert(d.μ, d.λ), tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(LRMoE.InverseGaussianExpert(d.μ, d.λ), tl, yl, yu, tu)

    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(yl, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    z_zero_e_lat = z_e_lat .* EM_E_z_zero_lat(tl, p_old, expert_tn_bar_pos)
    z_pos_e_lat = z_e_lat .- z_zero_e_lat
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, z_zero_e_lat, z_pos_e_lat, k_e)

    # Update parameters: call its positive part
    tmp_exp = InverseGaussianExpert(d.μ, d.λ)
    tmp_update = EM_M_expert(tmp_exp,
                            tl, yl, yu, tu,
                            exposure,
                            z_pos_e_obs, z_pos_e_lat, k_e,
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIInverseGaussianExpert(p_new, tmp_update.μ, tmp_update.λ)
end

## EM: M-Step, exact observations
function EM_M_expert_exact(d::ZIInverseGaussianExpert,
                     ye, exposure,
                     z_e_obs;
                     penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Update zero probability
    expert_ll_pos = expert_ll_exact.(LRMoE.InverseGaussianExpert(d.μ, d.λ), ye)

    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(ye, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, 0.0, 0.0, 0.0)

    # Update parameters: call its positive part
    tmp_exp = InverseGaussianExpert(d.μ, d.λ)
    tmp_update = EM_M_expert_exact(tmp_exp,
                            ye, exposure,
                            z_pos_e_obs;
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIInverseGaussianExpert(p_new, tmp_update.μ, tmp_update.λ)
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
# NOTE: Since  the internals have changed in julia v1.3.0, we need VERSION checks inside
# these functions.

# Helper for `pmatch`: Mirrors Base.PCRE.exec
function _pexec(re, subject, offset, options, match_data)
    rc = ccall((:pcre2_match_8, Base.PCRE.PCRE_LIB), Cint,
               (Ptr{Cvoid}, Ptr{UInt8}, Csize_t, Csize_t, Cuint, Ptr{Cvoid}, Ptr{Cvoid}),
               re, subject, sizeof(subject), offset, options, match_data,
               @static if VERSION >= v"1.3-"
                   Base.PCRE.get_local_match_context()
               else
                   Base.PCRE.MATCH_CONTEXT[]
               end)
    # rc == -1 means no match, -2 means partial match.
    rc < -2 && error("PCRE.exec error: $(err_message(rc))")
    (rc >= 0 ||
        # Allow partial matches if they were requested in the options.
        ((options & Base.PCRE.PARTIAL_HARD != 0 || options & Base.PCRE.PARTIAL_SOFT != 0) && rc == -2))
end

# Helper to call `_pexec` w/ thread-safe match data: Mirrors Base.PCRE.exec_r_data
# Only used for Julia Version >= v"1.3-"
function _pexec_r_data(re, subject, offset, options)
    match_data = Base.PCRE.create_match_data(re)
    ans = _pexec(re, subject, offset, options, match_data)
    return ans, match_data
end

"""
   pmatch(r::Regex, s::AbstractString[, idx::Integer[, addopts]])

Variant of `Base.match` that supports partial matches (when `Base.PCRE.PARTIAL_HARD`
is set in `re.match_options`).
"""
function pmatch(re::Regex, str::Union{SubString{String}, String}, idx::Integer, add_opts::UInt32=UInt32(0))
    Base.compile(re)
    opts = re.match_options | add_opts
    @static if VERSION >= v"1.3-"
        # rc == -1 means no match, -2 means partial match.
        matched, data = _pexec_r_data(re.regex, str, idx-1, opts)
        if !matched
            Base.PCRE.free_match_data(data)
            return nothing
        end
        n = div(Base.PCRE.ovec_length(data), 2) - 1
        p = Base.PCRE.ovec_ptr(data)
        mat = SubString(str, unsafe_load(p, 1)+1, prevind(str, unsafe_load(p, 2)+1))
        cap = Option{SubString{String}}[unsafe_load(p,2i+1) == Base.PCRE.UNSET ? nothing :
                                            SubString(str, unsafe_load(p,2i+1)+1,
                                                      prevind(str, unsafe_load(p,2i+2)+1)) for i=1:n]
        off = Int[ unsafe_load(p,2i+1)+1 for i=1:n ]
        result = RegexMatch(mat, cap, unsafe_load(p,1)+1, off, re)
        Base.PCRE.free_match_data(data)
        return result
    else  #  Julia VERSION < v"1.3"
        # rc == -1 means no match, -2 means partial match.
        matched = _pexec(re.regex, str, idx-1, opts, re.match_data)
        if !matched
            return nothing
        end
        ovec = re.ovec
        n = div(length(ovec),2) - 1
        mat = SubString(str, ovec[1]+1, prevind(str, ovec[2]+1))
        cap = Option{SubString{String}}[ovec[2i+1] == PCRE.UNSET ? nothing :
                                            SubString(str, ovec[2i+1]+1,
                                                      prevind(str, ovec[2i+2]+1)) for i=1:n]
        off = Int[ ovec[2i+1]+1 for i=1:n ]
        RegexMatch(mat, cap, ovec[1]+1, off, re)
    end
end

pmatch(r::Regex, s::AbstractString) = pmatch(r, s, firstindex(s))
pmatch(r::Regex, s::AbstractString, i::Integer) = throw(ArgumentError(
    "regex matching is only available for the String type; use String(s) to convert"
))

"""
    partial(str::AbstractString)

Constructs a regular expression to perform partial matching.
"""
partial(str::AbstractString) = Regex(str, Base.DEFAULT_COMPILER_OPTS,
    Base.DEFAULT_MATCH_OPTS | Base.PCRE.PARTIAL_HARD)

"""
    exact(str::AbstractString)

Constructs a regular expression to perform exact maching.
"""
exact(str::AbstractString) = Regex(str)


text:
using TheNaNTrap
using Documenter

makedocs(;
    modules=[TheNaNTrap],
    authors="Julia Computing, Matt Bauman <mbauman@juliacomputing.com>, and contributors",
    repo="https://github.com/mbauman/TheNaNTrap.jl/blob/{commit}{path}#L{line}",
    sitename="TheNaNTrap.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mbauman.github.io/TheNaNTrap.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mbauman/TheNaNTrap.jl",
)


text:
"""
FDI with simple time delay.
Λ_func: function of `t`.
"""
struct DelayFDI <: AbstractFDI
    τ
    function DelayFDI(τ)
        τ == 0.0 ? @warn("no delay") : nothing
        @assert τ >= 0
        new(τ)
    end
end

function Estimate(fdi::DelayFDI, Λ_func)
    @unpack τ = fdi
    return function (t)
        Λ̂ = Λ_func(t - τ)
    end
end

# function State(fdi::DelayFDI)
#     return function (dim_input::Int)
#         ComponentArray(Λ̂=ones(dim_input) |> Diagonal |> Matrix)
#     end
# end

# """
# Λ: effectiveness matrix
# # Notes
# For delayed FDI (DelayFDI),
# the function of effectiveness matrix should be known.
# For the explanation of the below implementation,
# see https://discourse.julialang.org/t/differentialequations-jl-discrete-delay-problem/62829/10?u=ihany.
# # NOTICE
# For some technical issue,
# the matrix-form state is dealt with by ComponentArray (default setting).
# """
# function Dynamics!(fdi::DelayFDI, Λ_func)
#     @unpack τ = fdi
#     dim_input = Λ_func(0.0) |> length  # auxiliary evaluation
#     dynamics! = function (dX, X, p, t)
#         dX.Λ̂ .= X.Λ̂ - Λ_func(t - τ)
#         nothing
#     end
#     ODEFunction(dynamics!, mass_matrix=zeros(dim_input) |> Diagonal |> Matrix)
# end


text:
include("estimated-cases.jl")
include("confirmed-cases.jl")


text:
#=
Project:      MIRTdev
File:         image_geom_mri.jl
Description:  Image geom function
Author:       Rajas Gupta
Date Created: October 1st 2019
Date Update:  October 7th 2019

2007-02-22, Jeff Fessler, University of Michigan (image_geom_mri.m)
=#

export image_geom_mri

using Test: @test, @inferred
#using MIRT: image_geom, MIRT_image_geom

"""
    ig = image_geom_mri(varargin)

Same as `image_geom()` but default offsets are 0.5 so that
image pixel indices go from -N/2 to N/2-1.

"""
function image_geom_mri( ; kwargs...)
    return image_geom( ; offset_x = .5, offset_y = .5, offset_z = .5, kwargs...)
end

"""
    image_geom_mri(:test)
self test
"""
function image_geom_mri(test::Symbol)
    test != :test && throw(DomainError(test, "Not valid"))
    @test (@inferred image_geom_mri( ; nx = 64, dx = 2)) isa MIRT_image_geom
    true
end


text:
__precompile__()

module EmpiricalOrthogonalFunctions

using Statistics
using LinearAlgebra
using Logging
import Base.truncate

include("eoftypes.jl")
include("eoffuncs.jl")

export EmpiricalOrthogonalFunction,
    pcs,
    eofs,
    correlationmap,
    covariancemap,
    eigenvalues,
    variancefraction,
    totalanomalyvar,
    northtest,
    projectfield,
    reconstruct,
    orthorotation,
    orthorotation!,
    truncate,
    truncate!

end # module


text:
abstract type SamplingAlgorithm end

"""
GridSample{T}

T is the step dx for lb:dx:ub
"""
struct GridSample{T} <: SamplingAlgorithm
    dx::T
end

struct UniformSample <: SamplingAlgorithm end
struct SobolSample <: SamplingAlgorithm end
struct LatinHypercubeSample <: SamplingAlgorithm end

"""
LowDiscrepancySample{T}

T is the base for the sequence
"""
struct LowDiscrepancySample{T} <: SamplingAlgorithm
    base::T
end

struct RandomSample <: SamplingAlgorithm end

struct KroneckerSample{A,B} <: SamplingAlgorithm
    alpha::A
    s0::B
end

struct GoldenSample <: SamplingAlgorithm end

"""
sample(n,lb,ub,S::GridSample)

Returns a tuple containing numbers in a grid.

"""
function sample(n,lb,ub,S::GridSample)
    dx = S.dx
    if lb isa Number
        return vec(rand(lb:S.dx:ub,n))
    else
        d = length(lb)
        x = [[rand(lb[j]:dx[j]:ub[j]) for j = 1:d] for i in 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,lb,ub,::UniformRandom)

Returns a Tuple containing uniform random numbers.
"""
function sample(n,lb,ub,::UniformSample)
    if lb isa Number
        return rand(Uniform(lb,ub),n)
    else
        d = length(lb)
        x = [[rand(Uniform(lb[j],ub[j])) for j in 1:d] for i in 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,lb,ub,::SobolSampling)

Returns a Tuple containing Sobol sequences.
"""
function sample(n,lb,ub,::SobolSample)
    s = SobolSeq(lb,ub)
    skip(s,n)
    if lb isa Number
        return [next!(s)[1] for i = 1:n]
    else
        return Tuple.([next!(s) for i = 1:n])
    end
end

"""
sample(n,lb,ub,::LatinHypercube)

Returns a Tuple containing LatinHypercube sequences.
"""
function sample(n,lb,ub,::LatinHypercubeSample)
    d = length(lb)
    if lb isa Number
        x = vec(LHCoptim(n,d,1)[1])
        # x∈[0,n], so affine transform
        return @. (ub-lb) * x/(n) + lb
    else
        lib_out = float(LHCoptim(n,d,1)[1])
        # x∈[0,n], so affine transform column-wise
        @inbounds for c = 1:d
            lib_out[:,c] = (ub[c]-lb[c])*lib_out[:,c]/n .+ lb[c]
        end
        x = [lib_out[i,:] for i = 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,lb,ub,S::LowDiscrepancySample)

Low discrepancy sample:
- Dimension 1: Van der Corput sequence
- Dimension > 1: Halton sequence
If dimension d > 1, all bases must be coprime with each other.
"""
function sample(n,lb,ub,S::LowDiscrepancySample)
    d = length(lb)
    if d == 1
        #Van der Corput
        b = S.base
        x = zeros(Float32,n)
        for i = 1:n
            expansion = digits(i,base = b)
            L = length(expansion)
            val = zero(Float32)
            for k = 1:L
                val += expansion[k]*float(b)^(-(k-1)-1)
            end
            x[i] = val
        end
        # It is always defined on the unit interval, resizing:
        return @. (ub-lb) * x + lb
    else
        #Halton sequence
        x = zeros(Float32,n,d)
        for j = 1:d
            b = S.base[j]
            for i = 1:n
                val = zero(Float32)
                expansion = digits(i, base = b)
                L = length(expansion)
                val = zero(Float32)
                for k = 1:L
                    val += expansion[k]*float(b)^(-(k-1)-1)
                end
                x[i,j] = val
            end
        end
        #Resizing
        # x∈[0,1], so affine transform column-wise
        @inbounds for c = 1:d
            x[:,c] = (ub[c]-lb[c])*x[:,c] .+ lb[c]
        end

        y = [x[i,:] for i = 1:n]
        return Tuple.(y)
    end
end

"""
sample(n,d,D::Distribution)

Returns a Tuple containing numbers distributed as D
"""
function sample(n,d,D::Distribution)
    if d == 1
        return rand(D,n)
    else
        x = [[rand(D) for j in 1:d] for i in 1:n]
        return Tuple.(x)
    end
end

"""
sample(n,d,K::KroneckerSample)

Returns a Tuple containing numbers following the Kronecker sample
"""
function sample(n,lb,ub,K::KroneckerSample)
    d = length(lb)
    alpha = K.alpha
    s0 = K.s0
    if d == 1
        x = zeros(n)
        @inbounds for i = 1:n
            x[i] = (s0+i*alpha)%1
        end
        return @. (ub-lb) * x + lb
    else
        x = zeros(n,d)
        @inbounds for j = 1:d
            for i = 1:n
                x[i,j] = (s0[j] + i*alpha[j])%i
            end
        end
        #Resizing
        # x∈[0,1], so affine transform column-wise
        @inbounds for c = 1:d
            x[:,c] = (ub[c]-lb[c])*x[:,c] .+ lb[c]
        end

        y = [x[i,:] for i = 1:n]
        return Tuple.(y)
    end
end

function sample(n,lb,ub,G::GoldenSample)
    d = length(lb)
    if d == 1
        x = zeros(n)
        g = (sqrt(5)+1)/2
        a = 1.0/g
        for i = 1:n
            x[i] = (0.5+a*i)%1
        end
        return @. (ub-lb) * x + lb
    else
        x = zeros(n,d)
        for j = 1:d
            #Approximate solution of x^(d+1) = x + 1, a simple newton is good enough
            y = 2.0
            for s = 1:10
                g = (1+y)^(1/(j+1))
            end
            a = 1.0/g
            for i = 1:n
                x[i,j] = (0.5+a*i)%1
            end
        end
        @inbounds for c = 1:d
            x[:,c] = (ub[c]-lb[c])*x[:,c] .+ lb[c]
        end
        y = [x[i,:] for i = 1:n]
        return Tuple.(y)
    end
end


text:
function test_poly_mi_cvx(s)
    global solver = s

    test_dir = dirname(@__FILE__)

    @testset "Convex MIPolynomial Models" begin end
end


text:
####################################################################################
# Two-dimensional shapes with parametric boundary
####################################################################################

# FIXME: this should probably be called a disk
struct Circle <: AbstractEntity
    # dim = 2
    tag::Int
    center::SVector{2,Float64}
    radius::Float64
    boundary::Vector{ParametricEntity}
    function Circle(t,c,r,bnd)
        ent = new(t,c,r,bnd)
        global_add_entity!(ent)
        return ent
    end
end

function Circle(;center=(0, 0),radius=1)
    f          = (s) -> center .+ radius .* SVector(cospi(2 * s[1]), sinpi(2 * s[1]))
    domain     = HyperRectangle(0,1)
    ent        = ParametricEntity(f, domain)
    tag        = new_tag(2) # generate a unique tag for entities of dimension 2
    return Circle(tag,center, radius, [ent])
end
Base.in(pt,circ::Circle) = norm(pt .- circ.center) < circ.radius
key(ent::Circle) = (2,ent.tag)
geometric_dimension(ent::Circle) = 2
ambient_dimension(ent::Circle) = 2

struct Kite <: AbstractEntity
    # dim = 2
    tag::Int
    center::SVector{2,Float64}
    radius::Float64
    boundary::Vector{ParametricEntity}
    function Kite(t,c,r,bnd)
        ent = new(t,c,r,bnd)
        global_add_entity!(ent)
        return ent
    end
end

function Kite(;radius=1,center=(0, 0))
    f = (s) -> center .+ radius .* SVector(cospi(2 * s[1]) + 0.65 * cospi(4 * s[1]) - 0.65,
                              1.5 * sinpi(2 * s[1]))
    domain = HyperRectangle(0,1)
    surf   = ParametricEntity(f, domain)
    t      = new_tag(2)
    return Kite(t,center, radius, [surf])
end
key(ent::Kite) = (2,ent.tag)
geometric_dimension(ent::Kite) = 2

struct Droplet <: AbstractEntity
    # dim = 2
    tag::Int
    center::SVector{2,Float64}
    radius::Float64
    boundary::Vector{ParametricEntity}
    function Droplet(t,c,r,bnd)
        ent = new(t,c,r,bnd)
        global_add_entity!(ent)
        return ent
    end
end

function Droplet(;radius=1,center=(0,0))
    f = (s) -> center .+ radius .* SVector(2 * sinpi(s[1]),-sinpi(2*s[1]))
    domain = HyperRectangle(0,1)
    surf   = ParametricEntity(f, domain)
    t      = new_tag(2)
    return Droplet(t,center, radius, [surf])
end
key(ent::Droplet) = (2,ent.tag)
geometric_dimension(ent::Droplet) = 2

struct Boomerang <: AbstractEntity
    # dim = 2
    tag::Int
    center::SVector{2,Float64}
    radius::Float64
    boundary::Vector{ParametricEntity}
    function Boomerang(t,c,r,bnd)
        ent = new(t,c,r,bnd)
        global_add_entity!(ent)
        return ent
    end
end

function Boomerang(;radius=1,center=(0,0))
    f = (s) -> center .+ radius .* SVector(-2/3 * sinpi(3*s[1]),-sinpi(2*s[1]))
    domain = HyperRectangle(0,1)
    surf   = ParametricEntity(f, domain)
    t      = new_tag(2)
    return Droplet(t,center, radius, [surf])
end
geometric_dimension(ent::Boomerang) = 2

####################################################################################
# Three-dimensional shapes with parametric boundary
####################################################################################

# TODO: rename to Ball?
struct Sphere <: AbstractEntity
    # dim = 3
    tag::Int
    center::SVector{3,Float64}
    radius::Float64
    boundary::Vector{ParametricEntity}
    function Sphere(t,c,r,bnd)
        ent = new(t,c,r,bnd)
        global_add_entity!(ent)
        return ent
    end
end

function Sphere(;center=(0,0,0),radius=0.5)
    nparts = 6
    domain = HyperRectangle((-1,-1),(1,1))
    parts  = Vector{ParametricEntity}(undef,nparts)
    for id=1:nparts
        param     = (x) -> _sphere_parametrization(x[1],x[2],id,radius,center)
        parts[id] = ParametricEntity(param,domain)
    end
    tag = new_tag(3)
    return Sphere(tag,center,radius,parts)
end
Base.in(pt,sph::Sphere) = norm(pt .- sph.center) < sph.radius
geometric_dimension(ent::Sphere) = 3
ambient_dimension(ent::Sphere) = 3

struct Ellipsoid <: AbstractEntity
    tag::Int
    center::SVector{3,Float64}
    paxis::SVector{3,Float64}
    boundary::Vector{ParametricEntity}
    function Ellipsoid(t,c,paxis,bnd)
        ent = new(t,c,paxis,bnd)
        global_add_entity!(ent)
        return ent
    end
end

function Ellipsoid(;center=zeros(3),paxis=ones(3))
    nparts = 6
    domain = HyperRectangle((-1.,-1.),(1.,1.))
    parts  = Vector{ParametricEntity}(undef,nparts)
    for id=1:nparts
        param     = (x) -> _ellipsoid_parametrization(x[1],x[2],id,paxis,center)
        parts[id] = ParametricEntity(param,domain,[domain])
    end
    return Ellipsoid(center,paxis,parts)
end

struct Bean <: AbstractEntity
    # dim = 3
    tag::Int
    center::SVector{3,Float64}
    paxis::SVector{3,Float64}
    boundary::Vector{ParametricEntity}
    function Bean(t,c,paxis,bnd)
        ent = new(t,c,paxis,bnd)
        global_add_entity!(ent)
        return ent
    end
end
geometric_dimension(ent::Bean) = 3
ambient_dimension(ent::Bean)   = 3

function Bean(;center=(0,0,0),paxis=(1,1,1))
    nparts = 6
    domain = HyperRectangle((-1,-1),(1,1))
    parts  = Vector{ParametricEntity}(undef,nparts)
    for id=1:nparts
        param     = (x) -> _bean_parametrization(x[1],x[2],id,paxis,center)
        parts[id] = ParametricEntity(param,domain)
    end
    tag = new_tag(3)
    return Bean(tag,center,paxis,parts)
end

struct Cube <: AbstractEntity
    # dim = 3
    tag::Int
    low_corner::SVector{3,Float64}
    widths::SVector{3,Float64}
    boundary::Vector{ParametricEntity}
    function Cube(t,c,p,bnd)
        ent = new(t,c,p,bnd)
        global_add_entity!(ent)
        return ent
    end
end

function Cube(;center=SVector(0,0,0),widths=(2,2,2))
    nparts = 6
    domain = HyperRectangle((-1.,-1.),(1.,1.))
    parts  = Vector{ParametricEntity}(undef,nparts)
    for id=1:nparts
        param     = (x) -> _cube_parametrization(x[1],x[2],id,widths,center)
        parts[id] = ParametricEntity(param,domain)
    end
    tag = new_tag(3)
    return Cube(tag,center,widths,parts)
end
geometric_dimension(ent::Cube) = 3
ambient_dimension(ent::Cube) = 3

function _cube_parametrization(u,v,id,widths,center)
    @assert -1 ≤ u ≤ 1
    @assert -1 ≤ v ≤ 1
    if id==1
        x = SVector(1.,u,v)
    elseif id==2
        x = SVector(-u,1.,v)
    elseif id==3
        x = SVector(u,v,1.)
    elseif id==4
        x = SVector(-1.,-u,v)
    elseif id==5
        x = SVector(u,-1.,v)
    elseif id==6
        x = SVector(-u,v,-1.)
    end
    return @. center + widths/2 * x
end

function _sphere_parametrization(u,v,id,rad=1,center=SVector(0,0,0))
    @assert -1 ≤ u ≤ 1
    @assert -1 ≤ v ≤ 1
    # parametrization of 6 patches. First gets a point on the cube [-1,1] ×
    # [-1,1] × [-1,1], the maps it onto the sphere
    if id==1
        x = SVector(1.,u,v)
    elseif id==2
        x = SVector(-u,1.,v)
    elseif id==3
        x = SVector(u,v,1.)
    elseif id==4
        x = SVector(-1.,-u,v)
    elseif id==5
        x = SVector(u,-1.,v)
    elseif id==6
        x = SVector(-u,v,-1.)
    end
    return center .+ rad.*x./sqrt(u^2+v^2+1)
end

function _ellipsoid_parametrization(u,v,id,paxis,center)
    x = _sphere_parametrization(u,v,id)
    return x .* paxis .+ center
end

function _bean_parametrization(u,v,id,paxis,center)
    x̂ = _sphere_parametrization(u,v,id)
    a = 0.8; b = 0.8; alpha1 = 0.3; alpha2 = 0.4; alpha3=0.1
    x = SVector(
        a*sqrt(1.0-alpha3*cospi(x̂[3])).*x̂[1],
        -alpha1*cospi(x̂[3])+b*sqrt(1.0-alpha2*cospi(x̂[3])).*x̂[2],
        x̂[3])
    return x .* paxis .+ center
end

function _acorn_parametrization(u,v,id,radius,center,rot)
    Rx = @SMatrix [1 0 0;0 cos(rot[1]) sin(rot[1]);0 -sin(rot[1]) cos(rot[1])]
    Ry = @SMatrix [cos(rot[2]) 0 -sin(rot[2]);0 1 0;sin(rot[2]) 0 cos(rot[2])]
    Rz = @SMatrix [cos(rot[3]) sin(rot[3]) 0;-sin(rot[3]) cos(rot[3]) 0;0 0 1]
    R  = Rz*Ry*Rx;
    x = _sphere_parametrization(u,v,id)
    th,phi,_ = cart2sph(x...)
    r=0.6+sqrt(4.25+2*cos(3*(phi+pi/2)))
    x[1] = r.*cos(th).*cos(phi)
    x[2] = r.*sin(th).*cos(phi)
    x[3] = r.*sin(phi)
    x    = R*x
    return radius.*x .+ center
end

function _cushion_parametrization(u,v,id,radius,center,rot)
    Rx = @SMatrix [1 0 0;0 cos(rot[1]) sin(rot[1]);0 -sin(rot[1]) cos(rot[1])]
    Ry = @SMatrix [cos(rot[2]) 0 -sin(rot[2]);0 1 0;sin(rot[2]) 0 cos(rot[2])]
    Rz = @SMatrix [cos(rot[3]) sin(rot[3]) 0;-sin(rot[3]) cos(rot[3]) 0;0 0 1]
    R  = Rz*Ry*Rx;
    x = _sphere_parametrization(u,v,id)
    th,phi,_ = cart2sph(x...)
    r = sqrt(0.8+0.5*(cos(2*th)-1).*(cos(4*phi)-1));
    x[1] = r.*cos(th).*cos(phi)
    x[2] = r.*sin(th).*cos(phi)
    x[3] = r.*sin(phi)
    x    = R*x
    return radius.*x .+ center
end

max_stars_repo_name: WaveProp/ParametricSurfaces

max_issues_repo_name: WaveProp/ParametricSurfaces

max_forks_repo_name: WaveProp/ParametricSurfaces


text:
using JuLIP
using SKTB

# @testset "Compare SKTB.jl against QUIP" begin

beta = 10.0

TB = SKTB
at = bulk("Si", cubic = true) * 2
rattle!(at, 0.02)
tbm = TB.NRLTB.NRLTBModel(:Si, FermiDiracSmearing(beta), orbitals = :sp, cutoff=:original)

H, M = hamiltonian(tbm, at)
H = Array(H)
M = Array(M)

using PyCall
@pyimport quippy

function make_quip_tb(at, Nk)
    k, weight = TB.monkhorstpackgrid(cell(at), (0, 0, Nk))
    kpoint_xml =  "<KPoints N=\"$(length(weight))\">\n" *
        join(["<point weight=\"$w\">$kx $ky $kz</point>\n" for (w, (kx, ky, kz)) in zip(weight, k)]) *
        "</KPoints>"
    println(kpoint_xml)

    xml_str = "<params>" *
        readstring(open("/Users/ortner/gits/QUIP/share/Parameters/tightbind.parms.NRL_TB.Si.xml", "r")) *
        kpoint_xml * "</params>"

    pot = quippy.Potential("TB NRL-TB", param_str=xml_str)
    quip_tb = ASECalculator(pot)
    return quip_tb
end

function tb_matrices(calc, at, tbm)
   norb = 4
   H = zeros(norb * length(at), norb  * length(at))
   M = zeros(norb * length(at), norb  * length(at))
   qsi = quippy.Atoms(at.po)
   calc.po[:calc_tb_matrices](qsi, hd=H, sd=M)

   # permute
   for n = 1:length(at)
      In = TB.indexblock(n, tbm.H)[2:4]
      Inp = [In[3], In[1], In[2]]
      H[:, In] = H[:, Inp]
      M[:, In] = M[:, Inp]
      H[In, :] = H[Inp, :]
      M[In, :] = M[Inp, :]
   end

   return H, M
end

quip_tb = make_quip_tb(at, 0)
Hq, Mq = tb_matrices(quip_tb, at, tbm)

M, H, Mq, Hq = real(M), real(H), real(Mq), real(Hq)

@show norm(H - Hq, Inf)
@test norm(H - Hq, Inf) < 1e-5
@show norm(M - Mq, Inf)
@test norm(M - Mq, Inf) < 1e-5

E = sort(real(eigvals(H, M)))
Eq = sort(real(eigvals(Hq, Mq)))
@show norm(E - Eq, Inf)
@test norm(E - Eq, Inf) < 1e-5

# println("JuLIP H: " )
# display(round(H[1:12, 1:16], 3)); println()
# println("QUIP H: ")
# display(round(Hq[1:12,1:16], 3)); println()
#
# println("JuLIP M: " )
# display(round(M[1:12, 1:16], 3)); println()
# println("QUIP M: ")
# display(round(Mq[1:12,1:16], 3)); println()

# end


text:
resolve_keys(key::Symbol, keys_children) = Symbol.((string(key) * "_") .* string.(keys_children))

Flux.Parallel(connection, layers::Vector) = Flux.Parallel(connection, layers...)

function fmap_custom(f, x, k=nothing; exclude=isleaf, cache=IdDict())
    # Because of https://github.com/FluxML/Flux.jl/blob/878b39cfd07b2366cd86f6a2b5b820990968c0ff/src/layers/basic.jl#L44
    k = (x isa Chain || k isa FChain) ? (k === nothing ? Symbol("_layers") : Symbol(string(k) * "_layers")) :
        (k === nothing ? Symbol() : k)

    haskey(cache, x) && return cache[x]
    y = if exclude(x)
        f(x, k)
    else
        func, re = functor(x)
        ks = resolve_keys(k, keys(func))
        res = re(map((x, y) -> first(fmap_custom(f, x, y; exclude=exclude, cache=cache)), values(func), ks))
        res = func isa NTuple ? Tuple(res) : res
    end
    x !== nothing && (cache[x] = (y, k))
    return y, k
end

"""
    destructure_parameters_states(model)
    destructure_parameters_states(model, ps)

# Arguments

Destructures the model and returns 3 values:

    - `parameters`: Trainable parameters of the `model` (by default all arrays in the @functor model are trainable. This can be changed by defining the `Flux.trainable` function for the `model`)
    - `states`: Arrays accessible by `fmap` but not trainable parameters
    - `restructure closure`: Can be used to restore the model given parameter and state vectors

Pass an `IdSet` to `ps` to only include certain `parameters`. Rest will be treated as `states`.

See also: [`destructure_parameters`](@ref)
"""
function destructure_parameters_states(m, ps::Zygote.IdSet=Flux.params(m).params)
    parameters = Zygote.Buffer([])
    states = Zygote.Buffer([])
    parameters_keys = Dict{Symbol,Function}()
    states_keys = Dict{Symbol,Function}()

    fmap_custom(m) do x, k
        if x isa AbstractArray{<:Number}
            v = vec(x)
            if x ∈ ps
                parameters_keys[k] = () -> zero(v)
                push!(parameters, v)
            else
                states_keys[k] = () -> zero(v)
                push!(states, v)
            end
        end
        return x
    end

    return (vcat(copy(parameters)...), vcat(copy(states)...),
            (p, s) -> _restructure_parameters_states(m, parameters_keys, states_keys, p, s))
end

function destructure_parameters_states(m, parameters_keys, states_keys)
    parameters = Zygote.Buffer([])
    states = Zygote.Buffer([])

    fmap_custom(m) do x, k
        if x === nothing
            if k ∈ keys(parameters_keys)
                push!(parameters, parameters_keys[k]())
            elseif k ∈ keys(states_keys)
                push!(states, states_keys[k]())
            end
        elseif x isa AbstractArray{<:Number}
            if k ∈ keys(parameters_keys)
                push!(parameters, vec(x))
            elseif k ∈ keys(states_keys)
                push!(states, vec(x))
            else
                error("Unknown Key Encountered: " * k)
            end
        end
        return x
    end

    return (vcat(copy(parameters)...), vcat(copy(states)...),
            (p, s) -> _restructure_parameters_states(m, parameters_keys, states_keys, p, s))
end

function _restructure_parameters_states(m, parameters_keys, states_keys, parameters, states)
    i, j = 0, 0

    m̄ = fmap_custom(m) do x, k
        !(x isa AbstractArray{<:Number}) && return x
        if k ∈ keys(parameters_keys)
            x = reshape(parameters[i .+ (1:length(x))], size(x))
            i += length(x)
        elseif k ∈ keys(states_keys)
            x = reshape(states[j .+ (1:length(x))], size(x))
            j += length(x)
        else
            error("Unknown Key Encountered: " * k)
        end
        return x
    end[1]

    exp_vector_size = i + j
    recv_vector_size = length(parameters) + length(states)
    recv_vector_size == exp_vector_size || @warn "Expected $(exp_vector_size) params, got $(recv_vector_size)"
    return m̄
end

Zygote.@adjoint function _restructure_parameters_states(m, parameters_keys, states_keys, parameters, states)
    m̄ = _restructure_parameters_states(m, parameters_keys, states_keys, parameters, states)
    numel = length(parameters) + length(states)
    function _restructure_parameters_states_pullback(dm)
        ∇parameters, ∇states, _ = destructure_parameters_states(dm, parameters_keys, states_keys)
        recv_numel = length(∇parameters) + length(∇states)
        numel != recv_numel && @warn "Expected $(numel) params, got $(recv_numel)"
        return (nothing, nothing, nothing, ∇parameters, ∇states)
    end
    return m̄, _restructure_parameters_states_pullback
end

"""
    destructure_parameters(model)
    destructure_parameters(model, ps)

# Arguments

Destructures the model and returns 2 values:

    - `parameters`: Trainable parameters of the `model` (by default all arrays in the @functor model are trainable. This can be changed by defining the `Flux.trainable` function for the `model`)
    - `restructure closure`: Can be used to restore the model given parameter and state vectors

Pass an `IdSet` to `ps` to only include certain `parameters`.

See also: [`destructure_parameters_states`](@ref)
"""
function destructure_parameters(m, ps::Zygote.IdSet=Flux.params(m).params)
    parameters = Zygote.Buffer([])
    parameters_keys = Dict{Symbol,Function}()

    fmap_custom(m) do x, k
        if x isa AbstractArray{<:Number}
            if x ∈ ps
                v = vec(x)
                parameters_keys[k] = () -> zero(v)
                push!(parameters, v)
            end
        end
        return x
    end

    return (vcat(copy(parameters)...), p -> _restructure_parameters(m, parameters_keys, p))
end

function destructure_parameters(m, parameters_keys)
    parameters = Zygote.Buffer([])

    fmap_custom(m) do x, k
        if x === nothing
            k ∈ keys(parameters_keys) && push!(parameters, parameters_keys[k]())
        elseif x isa AbstractArray{<:Number}
            k ∈ keys(parameters_keys) && push!(parameters, vec(x))
        end
        return x
    end

    return (vcat(copy(parameters)...), p -> _restructure_parameters(m, parameters_keys, p))
end

function _restructure_parameters(m, parameters_keys, parameters)
    i = 0

    m̄ = fmap_custom(m) do x, k
        !(x isa AbstractArray{<:Number}) && return x
        if k ∈ keys(parameters_keys)
            x = reshape(parameters[i .+ (1:length(x))], size(x))
            i += length(x)
        end
        return x
    end[1]

    exp_vector_size = i
    recv_vector_size = length(parameters)
    recv_vector_size == exp_vector_size || @warn "Expected $(exp_vector_size) params, got $(recv_vector_size)"
    return m̄
end

Zygote.@adjoint function _restructure_parameters(m, parameters_keys, parameters)
    m̄ = _restructure_parameters(m, parameters_keys, parameters)
    numel = length(parameters)
    function _restructure_parameters_pullback(dm)
        ∇parameters, _ = destructure_parameters(dm, parameters_keys)
        recv_numel = length(∇parameters)
        numel != recv_numel && @warn "Expected $(numel) params, got $(recv_numel)"
        return (nothing, nothing, ∇parameters)
    end
    return m̄, _restructure_parameters_pullback
end


text:
export example_domain, init_nnfem, example_global_data, example_static_domain1, find_boundary

@doc raw"""
init_nnfem(domain::Domain)

Prepares `domain` for use in custom operators.
"""
function init_nnfem(domain::Domain)
    LIB = joinpath(@__DIR__, "../../deps/CustomOp/build/libDataLib")

    if Sys.iswindows()
        LIB = joinpath(@__DIR__, "../../deps/CustomOp/build/DataLib")
    end

    LIB = abspath(LIB)
    tf.load_op_library(LIB*"."*dlext) # ensure that the tensorflow library is first loaded

    @eval begin
        ccall((:init_mesh, $LIB), Cvoid, ())

        for (iele,e) in enumerate($domain.elements)

            el_eqns = getEqns($domain, iele)
            el_eqns_active = findall(el_eqns .>= 1)
            ccall((:create_mesh, $LIB), Cvoid,
            (
            Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint,
            Ptr{Cdouble}, Cint, Ptr{Cint}
            ), Int32.(e.elnodes[:]), e.coords[:], vcat([x[:] for x in e.dhdx]...), e.weights, hcat([x[:] for x in e.hs]...), Int32(length(e.elnodes)), Int32(length(e.weights)),
            Int32.(el_eqns_active), Int32(length(el_eqns_active)), Int32.(el_eqns))
        end

        ccall((:create_domain, $LIB), Cvoid,
        (
        Ptr{Cdouble}, Cint, Cint, Cint
        ), $domain.nodes[:], Int32.($domain.nnodes), Int32.($domain.neqs), Int32.($domain.neles))
        ccall((:print_mesh, $LIB), Cvoid, ())

    end
end

function example_domain(m::Int64 = 10, n::Int64 = 10, h::Float64 = 1.0)
    # Create a very simple mesh
    elements = []
    prop = Dict("name"=> "PlaneStrain", "rho"=> 1.0, "E"=> 2.0, "nu"=> 0.35)
    coords = zeros((m+1)*(n+1), 2)
    for j = 1:n
        for i = 1:m
            idx = (m+1)*(j-1)+i
            elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
            ngp = 2
            nodes = [
            (i-1)*h (j-1)*h
            i*h (j-1)*h
            i*h j*h
            (i-1)*h j*h
            ]
            coords[elnodes, :] = nodes
            push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
        end
    end

    # fixed on the bottom, push on the right
    EBC = zeros(Int64, (m+1)*(n+1), 2)
    FBC = zeros(Int64, (m+1)*(n+1), 2)
    g = zeros((m+1)*(n+1), 2)
    f = zeros((m+1)*(n+1), 2)

    for i = 2:m
        idx = n*(m+1)+i
        FBC[idx,:] .= -1
        FBC[i,:] .= -1
    end

    for j = 1:n+1
        idx = (j-1)*(m+1)+1
        FBC[idx,:] .= -1
        idx = (j-1)*(m+1)+m+1
        FBC[idx,:] .= -1
    end

    ndims = 2
    domain = Domain(coords, elements, ndims, EBC, g, FBC, f)
end

function example_static_domain1(m::Int64 = 10, n::Int64 = 10, h::Float64 = 1.0)
    # Create a very simple mesh
    elements = []
    prop = Dict("name"=> "Scalar1D", "kappa"=>1.0 )
    coords = zeros((m+1)*(n+1), 2)
    for j = 1:n
        for i = 1:m
            idx = (m+1)*(j-1)+i
            elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
            ngp = 2
            nodes = [
                (i-1)*h (j-1)*h
                i*h (j-1)*h
                i*h j*h
                (i-1)*h j*h
            ]
            coords[elnodes, :] = nodes
            push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
        end
    end

    # free boundary on all sides
    EBC = zeros(Int64, (m+1)*(n+1))
    FBC = zeros(Int64, (m+1)*(n+1))
    g = zeros((m+1)*(n+1))
    f = zeros((m+1)*(n+1))

    for j = 1:n+1
        FBC[(j-1)*(m+1) + 1] = -1
        FBC[(j-1)*(m+1) + m + 1] = -1
    end

    for i = 1:m+1
        FBC[  i ] = -1
        FBC[ n*(m+1) + i ] = -1
    end

    ndims = 1
    domain = StaticDomain1(coords, elements, EBC, g, FBC, f)
end

"""
    example_global_data(domain::Domain)

Returns an examplary global data structure. Initial displacement, velocity, and acceleration are zero. Time-dependent displacement and force function are set to `missing`.
"""
function example_global_data(domain::Domain)
    Dstate = zeros(domain.neqs)
    state = zeros(domain.neqs)
    velo = zeros(domain.neqs)
    acce = zeros(domain.neqs)
    gt = nothing
    ft = nothing
    globdat = GlobalData(state, Dstate, velo, acce, domain.neqs, gt, ft)
end

"""
    find_boundary(nodes::Array{Float64, 2}, elements::Array{Int64, 2})

Given nodes and elements arrays, this function returns the boundary indices.
"""
function find_boundary(nodes::Array{Float64, 2}, elements::Array{Int64, 2})
    ef = (i, j)->(min(i, j), max(i, j))
    push_pop! = x -> (x in eset ? pop!(eset, x) : push!(eset, x))
    eset = Set([])
    for k = 1:size(elements, 1)
        push_pop!(ef(elements[k,1], elements[k,2]))
        push_pop!(ef(elements[k,2], elements[k,3]))
        push_pop!(ef(elements[k,3], elements[k,4]))
        push_pop!(ef(elements[k,4], elements[k,1]))
    end
    nset = Set(Int64[])
    for e in collect(eset)
        push!(nset, e[1])
        push!(nset, e[2])
    end
    return collect(nset)
end


text:
# # Using External Libraries

# It is straightforward to make use of external libraries
# written in other language. For example, the following
# codes import modules from the [OpenFermion](https://github.com/quantumlib/OpenFermion)
# and construct a qubit Hamiltonian for molecules as quantum blocks in Yao.

using Yao
using LinearAlgebra
using PyCall

# Install OpenFermion and PySCF if you don't have them
pip = pyimport("pip._internal.main")
pip.main(["install", "pyscf", "openfermion", "openfermionpyscf"])

# First we import hamiltonians from OpenFermion and PySCF

of_hamil = pyimport("openfermion.hamiltonians")
of_trsfm = pyimport("openfermion.transforms")
of_pyscf = pyimport("openfermionpyscf")

# following the OpenFermion and PySCF tutorial, we define molecules from data.

diatomic_bond_length = 1.0
geometry = [("H", (0., 0., 0.)), ("H", (0., 0., diatomic_bond_length))]
basis = "sto-3g"
multiplicity = 1
charge = 0
description = string(diatomic_bond_length)

molecule = of_hamil.MolecularData(geometry, basis, multiplicity, charge, description)
molecule = of_pyscf.run_pyscf(molecule,run_scf=1,run_fci=1)

# Then obtain its Hamiltonian

m_h = molecule.get_molecular_hamiltonian()
nbits = m_h.n_qubits

jw_h = of_trsfm.jordan_wigner(of_trsfm.get_fermion_operator(m_h))

# now, let's construct the corresponding Hamiltonian in Yao

function yao_hamiltonian(nbits, jw_h)
    gates = Dict("X"=>X, "Y"=>Y, "Z"=>Z)
    h = Add{nbits}()
    for (k, v) in jw_h.terms
        op = v*put(nbits, 1=>I2)
        for t in k
            site, opname = t
            op = op*put(nbits, site+1=>gates[opname])
        end
        push!(h, op)
    end
    h
end

# Yah!
yao_h = yao_hamiltonian(nbits, jw_h)

# now let's try to calculate its eigen values
w, v = eigen(Matrix(yao_h))

max_stars_repo_name: logankilpatrick/tutorials

max_issues_repo_name: logankilpatrick/tutorials

max_forks_repo_name: logankilpatrick/tutorials


text:
module ConstraintSolver

using MatrixNetworks
using JSON

CS = ConstraintSolver

mutable struct Variable
    idx         :: Int
    from        :: Int
    to          :: Int
    first_ptr   :: Int
    last_ptr    :: Int
    values      :: Vector{Int}
    indices     :: Vector{Int}
    offset      :: Int
    min         :: Int
    max         :: Int
    changes     :: Vector{Vector{Tuple{Symbol,Int64,Int64,Int64}}}
end

mutable struct CSInfo
    pre_backtrack_calls :: Int
    backtracked         :: Bool
    backtrack_fixes     :: Int
    in_backtrack_calls  :: Int
    backtrack_reverses  :: Int
end

abstract type Constraint
end

abstract type ObjectiveFunction
end

mutable struct MinMaxObjective <: ObjectiveFunction
    fct                 :: Function
    indices             :: Vector{Int}
end

mutable struct BasicConstraint <: Constraint
    idx                 :: Int
    fct                 :: Function
    indices             :: Vector{Int}
    pvals               :: Vector{Int}
    BasicConstraint() = new()
end

mutable struct LinearVariables
    indices             :: Vector{Int}
    coeffs              :: Vector{Int}
end

mutable struct LinearConstraint <: Constraint
    idx                 :: Int
    fct                 :: Function
    indices             :: Vector{Int}
    pvals               :: Vector{Int}
    coeffs              :: Vector{Int}
    operator            :: Symbol
    rhs                 :: Int
    in_all_different    :: Bool
    LinearConstraint() = new()
end

mutable struct BacktrackObj
    idx                 :: Int
    parent_idx          :: Int
    depth               :: Int
    status              :: Symbol
    variable_idx        :: Int
    pval                :: Int
    constraint_idx      :: Vector{Int}
    best_bound          :: Int

    BacktrackObj() = new()
end

mutable struct TreeLogNode
    id              :: Int
    status          :: Symbol
    best_bound      :: Int
    step_nr         :: Int
    var_idx         :: Int
    set_val         :: Int
    var_states      :: Dict{Int64,Vector{Int64}}
    var_changes     :: Dict{Int64,Vector{Tuple{Symbol, Int64, Int64, Int64}}}
    children        :: Vector{TreeLogNode}
    TreeLogNode() = new()
end

mutable struct CoM
    init_search_space   :: Vector{Variable}
    search_space        :: Vector{Variable}
    subscription        :: Vector{Vector{Int}}
    constraints         :: Vector{Constraint}
    bt_infeasible       :: Vector{Int}
    c_backtrack_idx     :: Int
    backtrack_vec       :: Vector{BacktrackObj}
    sense               :: Symbol
    objective           :: ObjectiveFunction
    best_sol            :: Int # Objective of the best solution
    best_bound          :: Int # Overall best bound
    solutions           :: Vector{Int} # saves only the id to the BacktrackObj
    info                :: CSInfo
    input               :: Dict{Symbol,Any}
    logs                :: Vector{TreeLogNode}

    CoM() = new()
end

include("printing.jl")
include("logs.jl")
include("Variable.jl")
include("objective.jl")
include("linearcombination.jl")
include("all_different.jl")
include("eq_sum.jl")
include("equal.jl")
include("not_equal.jl")

function init()
    com = CoM()
    com.constraints         = Vector{Constraint}()
    com.subscription        = Vector{Vector{Int}}()
    com.init_search_space   = Vector{Variable}()
    com.search_space        = Vector{Variable}()
    com.bt_infeasible       = Vector{Int}()
    com.c_backtrack_idx     = 1
    com.sense               = :None # means no objective
    com.backtrack_vec       = Vector{BacktrackObj}()
    com.solutions           = Vector{Int}()
    com.info                = CSInfo(0, false, 0, 0, 0)
    com.input               = Dict{Symbol, Bool}()
    com.logs                = Vector{TreeLogNode}()
    return com
end

function add_var!(com::CS.CoM, from::Int, to::Int; fix=nothing)
    ind = length(com.search_space)+1
    changes = Vector{Vector{Tuple{Symbol,Int64,Int64,Int64}}}()
    push!(changes, Vector{Tuple{Symbol,Int64,Int64,Int64}}())
    var = Variable(ind, from, to, 1, to-from+1, from:to, 1:to-from+1, 1-from, from, to, changes)
    if fix !== nothing
        fix!(com, var, fix)
    end
    push!(com.search_space, var)
    push!(com.subscription, Int[])
    push!(com.bt_infeasible, 0)
    return var
end

function fulfills_constraints(com::CS.CoM, index, value)
    # variable doesn't have any constraint
    if index > length(com.subscription)
        return true
    end
    constraints = com.constraints[com.subscription[index]]
    feasible = true
    for constraint in constraints
        feasible = constraint.fct(com, constraint, value, index)
        if !feasible
            break
        end
    end
    return feasible
end

"""
    fixed_vs_unfixed(search_space, indices)

Returns the fixed_vals as well as the unfixed_indices
"""
function fixed_vs_unfixed(search_space, indices)
    # get all values which are fixed
    fixed_vals = Int[]
    unfixed_indices = Int[]
    for (i,ind) in enumerate(indices)
        if isfixed(search_space[ind])
            push!(fixed_vals, value(search_space[ind]))
        else
            push!(unfixed_indices, i)
        end
    end
    return fixed_vals, unfixed_indices
end

"""
    set_pvals!(com::CS.CoM, constraint::Constraint)

Computes the possible values inside this constraint and sets it as constraint.pvals
"""
function set_pvals!(com::CS.CoM, constraint::Constraint)
    indices = constraint.indices
    variables = Variable[v for v in com.search_space[indices]]
    pvals_intervals = Vector{NamedTuple}()
    push!(pvals_intervals, (from = variables[1].from, to = variables[1].to))
    for (i,ind) in enumerate(indices)
        extra_from = variables[i].from
        extra_to   = variables[i].to
        comp_inside = false
        for cpvals in pvals_intervals
            if extra_from >= cpvals.from && extra_to <= cpvals.to
                # completely inside the interval already
                comp_inside = true
                break
            elseif extra_from >= cpvals.from && extra_from <= cpvals.to
                extra_from = cpvals.to+1
            elseif extra_to <= cpvals.to && extra_to >= cpvals.from
                extra_to = cpvals.from-1
            end
        end
        if !comp_inside && extra_to >= extra_from
            push!(pvals_intervals, (from = extra_from, to = extra_to))
        end
    end
    pvals = collect(pvals_intervals[1].from:pvals_intervals[1].to)
    for interval in pvals_intervals[2:end]
        pvals = vcat(pvals, collect(interval.from:interval.to))
    end
    constraint.pvals = pvals
end

"""
    add_constraint!(com::CS.CoM, constraint::Constraint)

Add a constraint to the model
"""
function add_constraint!(com::CS.CoM, constraint::Constraint)
    constraint.idx = length(com.constraints)+1
    push!(com.constraints, constraint)
    set_pvals!(com, constraint)
    for (i,ind) in enumerate(constraint.indices)
        push!(com.subscription[ind], constraint.idx)
    end
end

function set_objective!(com::CS.CoM, sense::Symbol, objective::ObjectiveFunction)
    if sense != :Min && sense != :Max
        throw(ErrorException("The objective sense must be :Min or :Max"))
    end
    com.sense = sense
    com.objective = objective
end

function get_weak_ind(com::CS.CoM)
    lowest_num_pvals = typemax(Int)
    biggest_inf = -1
    best_ind = -1
    biggest_dependent = typemax(Int)
    found = false

    for ind in 1:length(com.search_space)
        if !isfixed(com.search_space[ind])
            num_pvals = nvalues(com.search_space[ind])
            inf = com.bt_infeasible[ind]
            if inf >= biggest_inf
                if inf > biggest_inf || num_pvals < lowest_num_pvals
                    lowest_num_pvals = num_pvals
                    biggest_inf = inf
                    best_ind = ind
                    found = true
                end
            end
        end
    end
    return found, best_ind
end

function prune!(com::CS.CoM, prune_steps::Vector{Int})
    search_space = com.search_space
    for backtrack_idx in prune_steps
        for var in search_space
            for change in var.changes[backtrack_idx]
                fct_symbol = change[1]
                val = change[2]
                if fct_symbol == :fix
                    fix!(com, var, val; changes=false)
                elseif fct_symbol == :rm
                    rm!(com, var, val; changes=false)
                elseif fct_symbol == :remove_above
                    remove_above!(com, var, val; changes=false)
                elseif fct_symbol == :remove_below
                    remove_below!(com, var, val; changes=false)
                else
                    throw(ErrorException("There is no pruning function for $fct_symbol"))
                end
            end
        end
    end
end

"""
    prune!(com::CS.CoM; pre_backtrack=false)

Prune based on changes by initial solve or backtracking. The information for it is stored in each variable
Returns whether it's still feasible
"""
function prune!(com::CS.CoM; pre_backtrack=false)
    feasible = true
    prev_var_length = zeros(Int, length(com.search_space))
    constraint_idxs_vec = zeros(Int, length(com.constraints))
    # get all constraints which need to be called (only once)
    current_backtrack_id = com.c_backtrack_idx
    for var in com.search_space
        if length(var.changes[current_backtrack_id]) > 0
            inner_constraints = com.constraints[com.subscription[var.idx]]
            for constraint in inner_constraints
                constraint_idxs_vec[constraint.idx] = 1
            end
        end
    end

    current_level = 1
    level_increased = false
    # while we haven't called every constraint
    while true
        b_open_constraint = false
        # will be changed or b_open_constraint => false
        constraint = com.constraints[1]
        for ci in 1:length(com.constraints)
            if constraint_idxs_vec[ci] == current_level
                constraint = com.constraints[ci]
                constraint_idxs_vec[ci] = 0
                b_open_constraint = true
                break
            end
        end
        if !b_open_constraint && !level_increased
            level_increased = true
            current_level += 1
            continue
        end
        if !b_open_constraint && level_increased
            break
        end
        level_increased = false
        if all(v->isfixed(v), com.search_space[constraint.indices])
            continue
        end
        feasible = constraint.fct(com, constraint; logs = false)
        if !pre_backtrack
            com.info.in_backtrack_calls += 1
        else
            com.info.pre_backtrack_calls += 1
        end

        if !feasible
            break
        end

        # if we changed another variable increase the level of the constraints to call them later
        for var in com.search_space
            new_var_length = length(var.changes[current_backtrack_id])
            if new_var_length > prev_var_length[var.idx]
                prev_var_length[var.idx] = new_var_length
                inner_constraints = com.constraints[com.subscription[var.idx]]
                for constraint in inner_constraints
                    if constraint_idxs_vec[constraint.idx] < current_level
                        constraint_idxs_vec[constraint.idx] = current_level + 1
                    end
                end
            end
        end
    end
    return feasible
end

function single_reverse_pruning!(search_space, index::Int, prune_int::Int, prune_fix::Int)
    if prune_int > 0
        var = search_space[index]
        l_ptr = max(1,var.last_ptr)

        new_l_ptr = var.last_ptr + prune_int
        @views min_val = minimum(var.values[l_ptr:new_l_ptr])
        @views max_val = maximum(var.values[l_ptr:new_l_ptr])
        if min_val < var.min
            var.min = min_val
        end
        if max_val > var.max
            var.max = max_val
        end
        var.last_ptr = new_l_ptr
    end
    if prune_fix > 0
        var = search_space[index]
        f_ptr = max(1,var.first_ptr)
        var.last_ptr = prune_fix

        @views min_val = minimum(var.values[1:prune_fix])
        @views max_val = maximum(var.values[1:prune_fix])
        if min_val < var.min
            var.min = min_val
        end
        if max_val > var.max
            var.max = max_val
        end
        var.first_ptr = 1
    end
end

"""
    reverse_pruning!(com, backtrack_idx)

Reverse the changes made by a specific backtrack object
"""
function reverse_pruning!(com::CS.CoM, backtrack_idx::Int)
    search_space = com.search_space
    for var in search_space
        v_idx = var.idx

        for change in Iterators.reverse(var.changes[backtrack_idx])
            single_reverse_pruning!(search_space, v_idx, change[4], change[3])
        end
    end
end

function get_best_bound(com::CS.CoM; var_idx=0, val=0)
    if com.sense == :None
        return 0
    end
    return com.objective.fct(com, var_idx, val)
end

function checkout_from_to!(com::CS.CoM, from_idx::Int, to_idx::Int)
    backtrack_vec = com.backtrack_vec
    from = backtrack_vec[from_idx]
    to = backtrack_vec[to_idx]
    if to.parent_idx == from.idx
        return
    end
    reverse_pruning!(com, from.idx)

    prune_steps = Vector{Int}()
    # first go to same level if new is higher in the tree
    if to.depth < from.depth
        depth = from.depth
        parent_idx = from.parent_idx
        parent = backtrack_vec[parent_idx]
        while to.depth < depth
            reverse_pruning!(com, parent_idx)
            parent = backtrack_vec[parent_idx]
            parent_idx = parent.parent_idx
            depth -= 1
        end
        if parent_idx == to.parent_idx
            return
        else
            from = parent
        end
    elseif from.depth < to.depth
        depth = to.depth
        parent_idx = to.parent_idx
        parent = backtrack_vec[parent_idx]
        while from.depth < depth
            pushfirst!(prune_steps, parent_idx)
            parent = backtrack_vec[parent_idx]
            parent_idx = parent.parent_idx
            depth -= 1
        end

        to = parent
        if backtrack_vec[prune_steps[1]].parent_idx == from.parent_idx
            prune!(com, prune_steps)
            return
        end
    end
    @assert from.depth == to.depth
    # same diff but different parent
    # => level up until same parent
     while from.parent_idx != to.parent_idx
        reverse_pruning!(com, from.parent_idx)
        from = backtrack_vec[from.parent_idx]

        pushfirst!(prune_steps, to.parent_idx)
        to = backtrack_vec[to.parent_idx]
    end

    prune!(com, prune_steps)
end

function backtrack!(com::CS.CoM, max_bt_steps)
    found, ind = get_weak_ind(com)
    com.info.backtrack_fixes   = 1

    pvals = reverse!(values(com.search_space[ind]))
    dummy_backtrack_obj = BacktrackObj()
    dummy_backtrack_obj.status = :Close
    dummy_backtrack_obj.idx = 1
    dummy_backtrack_obj.variable_idx = 1

    backtrack_vec = com.backtrack_vec
    push!(backtrack_vec, dummy_backtrack_obj)

    # the first solve (before backtrack) has idx 1
    num_backtrack_objs = 1
    step_nr = 1
    for pval in pvals
        backtrack_obj = BacktrackObj()
        num_backtrack_objs += 1
        backtrack_obj.idx = num_backtrack_objs
        backtrack_obj.parent_idx = 1
        backtrack_obj.depth = 1
        backtrack_obj.status = :Open
        backtrack_obj.best_bound = get_best_bound(com; var_idx=ind, val=pval)
        backtrack_obj.variable_idx = ind
        backtrack_obj.pval = pval
        push!(backtrack_vec, backtrack_obj)
        for v in com.search_space
            push!(v.changes, Vector{Tuple{Symbol,Int64,Int64,Int64}}())
        end
        if com.input[:logs]
            push!(com.logs, log_one_node(com, length(com.search_space), num_backtrack_objs, -1))
        end
    end
    last_backtrack_id = 0

    started = true
    just_increased_depth = true

    while length(backtrack_vec) > 0
        step_nr += 1
        # get next open backtrack object
        l = 1
        if com.sense == :None
            backtrack_obj = backtrack_vec[l]
            while backtrack_obj.status != :Open
                l += 1
                if l > length(backtrack_vec)
                    break
                end
                backtrack_obj = backtrack_vec[l]
            end
        else
            sorted_backtrack_idxs_nvals = sortperm(backtrack_vec; by=o->o.depth, rev=true)
            sorted_backtrack_idxs_bb = sortperm(backtrack_vec[sorted_backtrack_idxs_nvals]; by=o->o.best_bound, rev=com.sense == :Max )
            sorted_backtrack_idxs = sorted_backtrack_idxs_nvals[sorted_backtrack_idxs_bb]

            backtrack_obj = backtrack_vec[sorted_backtrack_idxs[l]]
            while backtrack_obj.status != :Open
                l += 1
                if l > length(backtrack_vec)
                    break
                end
                backtrack_obj = backtrack_vec[sorted_backtrack_idxs[l]]
            end
        end

        # no open node => Infeasible
        if l > length(backtrack_vec)
            break
        end
        com.best_bound = backtrack_obj.best_bound
        ind = backtrack_obj.variable_idx

        com.c_backtrack_idx = backtrack_obj.idx

        if !started
            com.c_backtrack_idx = 0
            checkout_from_to!(com, last_backtrack_id, backtrack_obj.idx)
            com.c_backtrack_idx = backtrack_obj.idx
        end

        if com.input[:logs]
            com.logs[backtrack_obj.idx].step_nr = step_nr
        end

        started = false
        last_backtrack_id = backtrack_obj.idx

        just_increased_depth = false
        backtrack_obj.status = :Closed

        pval = backtrack_obj.pval

        # check if this value is still possible
        constraints = com.constraints[com.subscription[ind]]
        feasible = true
        for constraint in constraints
            feasible = constraint.fct(com, constraint, pval, ind)
            if !feasible
                break
            end
        end
        if !feasible
            continue
        end
        # value is still possible => set it
        fix!(com, com.search_space[ind], pval)

        for constraint in constraints
            feasible = constraint.fct(com, constraint; logs = false)
            if !feasible
                feasible = false
                break
            end
        end
        if !feasible
            com.info.backtrack_reverses += 1
            continue
        end

        # prune on changed values
        feasible = prune!(com)

        if !feasible
            com.info.backtrack_reverses += 1
            continue
        end

        found, ind = get_weak_ind(com)
        if !found
            com.best_sol = get_best_bound(com)
            if com.best_sol == com.best_bound
                push!(com.solutions, backtrack_obj.idx)
                return :Solved
            else
                continue
            end
        end

        if com.info.backtrack_fixes + 1 > max_bt_steps
            return :NotSolved
        end

        com.info.backtrack_fixes   += 1

        if com.input[:logs]
            com.logs[backtrack_obj.idx] = log_one_node(com, length(com.search_space), backtrack_obj.idx, step_nr)
        end

        pvals = reverse!(values(com.search_space[ind]))
        last_backtrack_obj = backtrack_vec[last_backtrack_id]
        for pval in pvals
            backtrack_obj = BacktrackObj()
            num_backtrack_objs += 1
            backtrack_obj.parent_idx = last_backtrack_obj.idx
            backtrack_obj.depth = last_backtrack_obj.depth + 1
            backtrack_obj.idx = num_backtrack_objs
            backtrack_obj.status = :Open
            backtrack_obj.best_bound = get_best_bound(com; var_idx=ind, val=pval)
            backtrack_obj.variable_idx = ind
            backtrack_obj.pval = pval
            push!(backtrack_vec, backtrack_obj)
            for v in com.search_space
                push!(v.changes, Vector{Tuple{Symbol,Int64,Int64,Int64}}())
            end
            if com.input[:logs]
                push!(com.logs, log_one_node(com, length(com.search_space), num_backtrack_objs, -1))
            end
        end
        just_increased_depth = true
    end
    if length(com.solutions) > 0
        println("Actually found one")
    end
    return :Infeasible
end

function arr2dict(arr)
    d = Dict{Int,Bool}()
    for v in arr
        d[v] = true
    end
    return d
end

function simplify!(com)
    # check if we have all_different and sum constraints
    # (all different where every value is used)
    b_all_different = false
    b_all_different_sum = false
    b_sum = false
    for constraint in com.constraints
        if nameof(constraint.fct) == :all_different
            b_all_different = true
            if length(constraint.indices) == length(constraint.pvals)
                b_all_different_sum = true
            end
        elseif nameof(constraint.fct) == :eq_sum
            b_sum = true
        end
    end
    # l = 0
    if b_all_different_sum && b_sum
        # for each all_different constraint
        # which can be formulated as a sum constraint
        # check which sum constraints are completely inside all different
        # which are partially inside
        # compute inside sum and total sum
        # println("lc :", length(com.constraints))
        n_constraints_before = length(com.constraints)
        for constraint_idx in 1:length(com.constraints)
            constraint = com.constraints[constraint_idx]

            if nameof(constraint.fct) == :all_different
                add_sum_constraint = true
                if length(constraint.indices) == length(constraint.pvals)
                    all_diff_sum = sum(constraint.pvals)
                    in_sum = 0
                    total_sum = 0
                    outside_indices = Int[]
                    cons_indices_dict = arr2dict(constraint.indices)
                    for variable_idx in keys(cons_indices_dict)
                        found_sum_constraint = false
                        for sub_constraint_idx in com.subscription[variable_idx]
                            # don't mess with constraints added later on
                            if sub_constraint_idx > n_constraints_before
                                continue
                            end
                            sub_constraint = com.constraints[sub_constraint_idx]
                            # it must be an equal constraint and all coefficients must be 1 otherwise we can't add a constraint
                            if nameof(sub_constraint.fct) == :eq_sum && all(c->c==1, sub_constraint.coeffs)
                                found_sum_constraint = true
                                total_sum += sub_constraint.rhs
                                all_inside = true
                                for sub_variable_idx in sub_constraint.indices
                                    if !haskey(cons_indices_dict, sub_variable_idx)
                                        all_inside = false
                                        push!(outside_indices, sub_variable_idx)
                                    else
                                        delete!(cons_indices_dict, sub_variable_idx)
                                    end
                                end
                                if all_inside
                                    in_sum += sub_constraint.rhs
                                end
                                break
                            end
                        end
                        if !found_sum_constraint
                            add_sum_constraint = false
                            break
                        end
                    end

                    # make sure that there are not too many outside indices
                    if add_sum_constraint && length(outside_indices) < 3
                        add_constraint!(com, sum(com.search_space[outside_indices]) == total_sum - all_diff_sum)
                    end
                end
            end
        end
    end
end

function set_in_all_different!(com::CS.CoM)
    for constraint in com.constraints
        if :in_all_different in fieldnames(typeof(constraint))
            if !constraint.in_all_different
                subscriptions_idxs = [[i for i in com.subscription[v]] for v in constraint.indices]
                intersects = intersect(subscriptions_idxs...)

                for i in intersects
                    if nameof(com.constraints[i].fct) == :all_different
                        constraint.in_all_different = true
                        break
                    end
                end
            end
        end
    end
end

function solve!(com::CS.CoM; backtrack=true, max_bt_steps=typemax(Int64), keep_logs=false)
    com.input[:logs] = keep_logs
    if keep_logs
        com.init_search_space = deepcopy(com.search_space)
    end

    set_in_all_different!(com)

    # check for better constraints
    simplify!(com)

    # check if all feasible even if for example everything is fixed
    feasible = true
    for constraint in com.constraints
        com.info.pre_backtrack_calls += 1
        feasible = constraint.fct(com, constraint)
        if !feasible
            feasible = false
            break
        end
    end

    if !feasible
        return :Infeasible
    end

    if all(v->isfixed(v), com.search_space)
        com.best_bound = get_best_bound(com)
        com.best_sol = com.best_bound
        return :Solved
    end
    feasible = prune!(com; pre_backtrack=true)

    com.best_bound = get_best_bound(com)
    if keep_logs
        push!(com.logs, log_one_node(com, length(com.search_space), 1, 1))
    end

    if !feasible
        return :Infeasible
    end

    if all(v->isfixed(v), com.search_space)
        com.best_sol = com.best_bound
        return :Solved
    end
    if backtrack
        com.info.backtracked = true
        return backtrack!(com, max_bt_steps)
    else
        @info "Backtracking is turned off."
        return :NotSolved
    end
end

export add_var!, add_constraint!, solve!, value, set_objective!

end # module


text:
mutable struct PQueue{T}
    min::Int
    nbins::Int
    content::Vector{Vector{T}}
end

# initialize new queue
PQueue{T}(nbins) where T = PQueue(nbins + 1, nbins, [Vector{T}() for _ in 1:nbins])
PQueue(nbins, item::T, weight=1) where T = enqueue!(PQueue{T}(nbins), item, weight)

Base.isempty(q::PQueue) = q.min > q.nbins

function enqueue!(q::PQueue, item, weight)
    push!(q.content[weight], item)
    q.min = min(q.min, weight)
    return q
end

function dequeue!(q::PQueue)
    elem = pop!(q.content[q.min])
    # increase smallestbin, if elem was last in bin
    while q.min ≤ q.nbins && isempty(q.content[q.min])
        q.min += 1
    end
    return elem
end


text:
LCR6xxx_Horizontal_Dict = Dict()

# Nr of acquisition points
function set_nr_pts(obj::T;npts=1e6) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.Horizontal.Maximize=0")
	write(obj,"VBS app.Acquisition.Horizontal.MaxSamples=$(Int64(npts))")
end
# Sampling rate
function set_srate(obj::T;sr=1e6) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.Horizontal.Maximize=1")
	write(obj,"VBS app.Acquisition.Horizontal.SampleRate=$(Int64(sr))")
end
# Acquisition duration
function set_hduration(obj::T;d=1e-3) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.Horizontal.HorScale=$(d/10)")
end
# Horizontal offset origin
function set_hoffs_div(obj::T;href::Int=1) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.Horizontal.HorOffsetOrigin=$href")
end


text:
module TeleMap

	abstract type TelemetricMap end;
	abstract type TelemetricEvent end;
	abstract type TelemetricEvents end;
	abstract type SpeedType end;

	import Dates, KernelDensity, Rmath, DataFrames, Query, RecipesBase

	export TelemetricMap, TelemetricEvent, TelemetricEvents, SpeedType, KinematicEvent, KinematicEvents, KinematicMap, DeviationEvents, MPH, KMH

	include("Kinematic.jl")

end # module


text:
using GeoData, Test, Statistics, Dates, Plots, DimensionalData, RasterDataSources
import ArchGDAL, NCDatasets
using GeoData: window, mode, span, sampling, name, bounds

include(joinpath(dirname(pathof(GeoData)), "../test/test_utils.jl"))

path = maybedownload("https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif")

@testset "GDALarray" begin
    gdalarray = geoarray(path; mappedcrs=EPSG(4326), name=:test)

    @testset "open" begin
        @test open(A -> A[Y=1], gdalarray) == gdalarray[:, 1, :]
    end

    @testset "read" begin
        A = read(gdalarray)
        @test A isa GeoArray
        @test parent(A) isa Array
    end

    @testset "array properties" begin
        @test size(gdalarray) == (514, 515, 1)
        @test gdalarray isa GDALarray{UInt8,3}
    end

    @testset "dimensions" begin
        @test length(dims(gdalarray, X)) == 514
        @test ndims(gdalarray) == 3
        @test dims(gdalarray) isa Tuple{<:X,<:Y,<:Band}
        @test mode(gdalarray, Band) == DimensionalData.Categorical(Ordered())
        @test span(gdalarray, (Y, X)) ==
            (Regular(-60.02213698319351), Regular(60.02213698319374))
        @test sampling(gdalarray, (Y, X)) ==
            (Intervals(Start()), Intervals(Start()))
        @test refdims(gdalarray) == ()
        # Bounds calculated in python using rasterio
        @test all(bounds(gdalarray, Y) .≈ (4224973.143255847, 4255884.5438021915))
        @test all(bounds(gdalarray, X) .≈ (-28493.166784412522, 2358.211624949061))
    end

    @testset "other fields" begin
        # This file has an inorrect missing value
        @test missingval(gdalarray) == nothing
        @test metadata(gdalarray) isa Metadata{:GDAL}
        @test basename(metadata(gdalarray).val[:filepath]) == "cea.tif"
        @test name(gdalarray) == :test
        @test label(gdalarray) == "test"
        @test units(gdalarray) == nothing
        @test mappedcrs(dims(gdalarray, Y)) == EPSG(4326)
        @test mappedcrs(dims(gdalarray, X)) == EPSG(4326)
        @test mappedcrs(gdalarray) == EPSG(4326)
        @test mappedcrs(gdalarray[Y(1)]) == EPSG(4326)
        @test_throws ErrorException mappedcrs(gdalarray[Y(1), X(1)])
        wkt = WellKnownText(GeoFormatTypes.CRS(),
          "PROJCS[\"unnamed\",GEOGCS[\"NAD27\",DATUM[\"North_American_Datum_1927\",SPHEROID[\"Clarke 1866\",6378206.4,294.978698213898,AUTHORITY[\"EPSG\",\"7008\"]],AUTHORITY[\"EPSG\",\"6267\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4267\"]],PROJECTION[\"Cylindrical_Equal_Area\"],PARAMETER[\"standard_parallel_1\",33.75],PARAMETER[\"central_meridian\",-117.333333333333],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")
        @test crs(dims(gdalarray, Y)) == wkt
        @test crs(dims(gdalarray, X)) == wkt
        @test crs(gdalarray) == wkt
        @test crs(gdalarray[Y(1)]) == wkt
        @test_throws ErrorException crs(gdalarray[Y(1), X(1)])
    end

    @testset "indexing" begin
        @test gdalarray[Band(1)] isa GeoArray{UInt8,2}
        @test gdalarray[Y(1), Band(1)] isa GeoArray{UInt8,1}
        @test gdalarray[X(1), Band(1)] isa GeoArray{UInt8,1}
        @test gdalarray[X(1), Y(1), Band(1)] isa UInt8
        @test gdalarray[1, 1, 1] isa UInt8
        # Value indexed in python with rasterio
        @test gdalarray[Y(21), X(21), Band(1)] == 82
    end

    @testset "methods" begin
        mean(gdalarray; dims=Y) == mean(data(gdalarray); dims=2)
    end

    @testset "selectors" begin
        # TODO verify the value with R/gdal etc
        @test gdalarray[Y(Contains(33.8)), X(Contains(-117.5)), Band(1)] isa UInt8
        @test gdalarray[Y(Between(33.7, 33.9)), Band(1)] isa GeoArray
    end

    @testset "conversion to GeoArray" begin
        geoA = gdalarray[X(1:50), Y(1:1), Band(1)]
        @test size(geoA) == (50, 1)
        @test eltype(geoA) <: UInt8
        @time geoA isa GeoArray{UInt8,1}
        @test dims(geoA) isa Tuple{<:X,Y}
        @test refdims(geoA) isa Tuple{<:Band}
        @test metadata(geoA) == metadata(gdalarray)
        @test missingval(geoA) == nothing
        @test name(geoA) == :test
    end

    @testset "save" begin
        gdalarray = GDALarray(path; mappedcrs=EPSG(4326), name=:test);

        @testset "2d" begin
            geoA = gdalarray[Band(1)]
            filename = tempname() * ".asc"
            write(filename, geoA)
            saved1 = GDALarray(filename; mappedcrs=EPSG(4326))[Band(1)];
            @test all(saved1 .== saved1 ≈ geoA)
            @test typeof(saved1) !== typeof(geoA)
            @test val(dims(saved1, X)) ≈ val(dims(geoA, X))
            @test val(dims(saved1, Y)) ≈ val(dims(geoA, Y))
            @test all(metadata.(dims(saved1)) .== metadata.(dims(geoA)))
            @test metadata(dims(saved1)[1]) == metadata(dims(geoA)[1])
            @test missingval(saved1) === missingval(geoA)
            @test refdims(saved1) == refdims(geoA)
        end

        @testset "3d, with subsetting" begin
            geoA2 = gdalarray[Y(Between(33.7, 33.9)),
                                  X(Between(-117.6, -117.4))]
            filename2 = tempname() * ".tif"
            write(filename2, geoA2)
            saved2 = GeoArray(GDALarray(filename2; name=:test, mappedcrs=EPSG(4326)))
            @test size(saved2) == size(geoA2) == length.(dims(saved2)) == length.(dims(geoA2))
            @test refdims(saved2) == refdims(geoA2)
            #TODO test a file with more metadata
            @test val(metadata(saved2))[:filepath] == filename2
            @test missingval(saved2) === missingval(geoA2)
            @test GeoData.name(saved2) == GeoData.name(geoA2)
            @test step(mode(dims(saved2, Y))) ≈ step(mode(dims(geoA2, Y)))
            @test step(mode(dims(saved2, X))) ≈ step(mode(dims(geoA2, X)))
            @test typeof(dims(saved2)) == typeof(dims(geoA2))
            @test all(val(dims(saved2, Band)) .≈ val(dims(geoA2, Band)))
            @test all(val(dims(saved2, X)) .≈ val(dims(geoA2, X)))
            @test all(val(dims(saved2, Y)) .≈ val(dims(geoA2, Y)))
            @test all(metadata.(dims(saved2)) .== metadata.(dims(geoA2)))
            @test data(saved2) == data(geoA2)
            @test typeof(saved2) == typeof(geoA2)
            filename3 = tempname() * ".tif"
            geoA3 = cat(gdalarray[Band(1)], gdalarray[Band(1)], gdalarray[Band(1)]; dims=Band(1:3))
            write(filename3, geoA3)
            saved3 = GeoArray(GDALarray(filename3; mappedcrs=EPSG(4326)))
            @test all(saved3 .== geoA3)
            @test val(dims(saved3, Band)) == 1:3
        end

        @testset "resave current" begin
            filename = tempname() * ".rst"
            write(filename, gdalarray)
            gdalarray2 = GDALarray(filename)
            write(gdalarray2)
            @test convert(GeoArray, GDALarray(filename)) == convert(GeoArray, gdalarray2)
        end

        @testset "to grd" begin
            write("testgrd.gri", gdalarray)
            grdarray = GRDarray("testgrd.gri")
            @test crs(grdarray) == convert(ProjString, crs(gdalarray))
            @test bounds(grdarray) == (bounds(gdalarray))
            @test val(dims(grdarray, Y)) == reverse(val(dims(gdalarray, Y)))
            @test val(dims(grdarray, X)) ≈ val(dims(gdalarray, X))
            @test all(GeoArray(grdarray) .== GeoArray(gdalarray))
        end

        @testset "from GeoArray" begin
            filename = tempname() * ".tiff"
            ga = GeoArray(rand(100, 200), (X, Y))
            reorder(ga, (X(ForwardIndex()), Y(ReverseIndex())))
            write(filename, ga)
            saved = GDALarray(filename)
            @test parent(saved[Band(1)]) == parent(ga)
            @test saved[1, 1, 1] == saved[At(1), At(1), At(1)]
            @test saved[100, 200, 1] == saved[At(100), At(200), At(1)]
            filename2 = tempname() * ".tif"
            ga2 = GeoArray(rand(100, 200), (X(101:200; mode=Sampled()), Y(1:200; mode=Sampled())))
            write(filename2, ga2)
            @test parent(reorder(GDALarray(filename2)[Band(1)], ForwardArray)) == ga2
        end

        # This needs netcdf bounds variables to work
        @testset "to netcdf" begin
            filename2 = tempname() * ".nc"
            write(filename2, gdalarray[Band(1)])
            saved = GeoArray(NCDarray(filename2; crs=crs(gdalarray)))
            @test size(saved) == size(gdalarray[Band(1)])
            @test saved ≈ reverse(gdalarray[Band(1)]; dims=Lat)
            clat, clon = DimensionalData.shiftlocus.(Ref(Center()), dims(gdalarray, (Lat, Lon)))
            @test mappedindex(clat) ≈ reverse(mappedindex(saved, Lat))
            @test mappedindex(clon) ≈ mappedindex(saved, Lon)
            @test all(mappedbounds(saved, Lon) .≈ mappedbounds(clon))
            @test all(mappedbounds(saved, Lat) .≈ mappedbounds(clat))
            @test projectedindex(clon) ≈ projectedindex(saved, Lon)
            @test all(projectedbounds(clon) .≈ projectedbounds(saved, Lon))
            # reason lat crs conversion is less accrurate than lon TODO investigate further
            @test all(map((a, b) -> isapprox(a, b; rtol=1e-6),
                projectedindex(gdalarray, Lat),
                reverse(projectedindex(DimensionalData.shiftlocus(Start(), dims(saved, Lat))))
            ))
            @test all(map((a, b) -> isapprox(a, b; rtol=1e-6), projectedbounds(saved, Lat),  projectedbounds(gdalarray, Lat)))
        end

    end

    @testset "show" begin
        sh = sprint(show, MIME("text/plain"), gdalarray)
        # Test but don't lock this down too much
        @test occursin("GDALarray", sh)
        @test occursin("Y", sh)
        @test occursin("X", sh)
        @test occursin("Band", sh)
    end

    @testset "plot" begin # TODO write some tests for this
        gdalarray |> plot
        gdalarray[Y(1)] |> plot
    end

    @testset "nodatavalue type matches the array type" begin
        A = geoarray(WorldClim{Climate}, :tavg; res="10m", month=1)
        @test typeof(missingval(A)) === eltype(A)
        @test missingval(A) === -3.4f38
    end

end

@testset "GDAL stack" begin
    gdalstack = stack((a=path, b=path))

    @testset "read" begin
        st = read(gdalstack)
        @test st isa GeoStack
        @test st.data isa NamedTuple
        @test first(st.data) isa GeoArray
        @test parent(first(st.data)) isa Array
    end

    @testset "child array properties" begin
        @test size(gdalstack[:a]) == (514, 515, 1)
        @test gdalstack[:a] isa GeoArray{UInt8,3}
    end

    @testset "indexing" begin
        @test gdalstack[:a][Y(2:3), X(1), Band(1)] == [0x00, 0x6b]
        @test gdalstack[:a][Y(1), X(1), Band(1)] == 0x00
        @test gdalstack[:b, Band(1)] == gdalstack[:b][Band(1)]
        @test typeof(gdalstack[:b, Band(1)]) == typeof(gdalstack[:b][Band(1)])
    end

    @testset "window" begin
        windowedstack = GDALstack((a=path, b=path); window=(Y(1:5), X(1:5), Band(1)))
        @test window(windowedstack) == (Y(1:5), X(1:5), Band(1))
        windowedarray = GeoArray(windowedstack[:a])
        @test windowedarray isa GeoArray{UInt8,2}
        @test length.(dims(windowedarray)) == (5, 5)
        @test size(windowedarray) == (5, 5)
        @test windowedarray[1:3, 2:2] == reshape([0x00, 0x00, 0x00], 3, 1)
        @test windowedarray[1:3, 2] == [0x00, 0x00, 0x00]
        @test windowedarray[1, 2] == 0x00
        windowedstack = GDALstack((a=path, b=path); window=(Y(1:5), X(1:5), Band(1:1)))
        windowedarray = windowedstack[:b]
        @test windowedarray[1:3, 2:2, 1:1] == reshape([0x00, 0x00, 0x00], 3, 1, 1)
        @test windowedarray[1:3, 2:2, 1] == reshape([0x00, 0x00, 0x00], 3, 1)
        @test windowedarray[1:3, 2, 1] == [0x00, 0x00, 0x00]
        @test windowedarray[1, 2, 1] == 0x00
        windowedstack = GDALstack((a=path, b=path); window=(Band(1),))
        windowedarray = GeoArray(windowedstack[:b])
        @test windowedarray[1:3, 2:2] == reshape([0x00, 0x00, 0x00], 3, 1)
        @test windowedarray[1:3, 2] == [0x00, 0x00, 0x00]
        @test windowedarray[1, 2] == 0x00
    end

    # Stack Constructors
    @testset "conversion to GeoStack" begin
        geostack = GeoStack(gdalstack)
        @test Symbol.(Tuple(keys(gdalstack))) == keys(geostack)
        smallstack = GeoStack(gdalstack; keys=(:a,))
        @test keys(smallstack) == (:a,)
    end

    if VERSION > v"1.1-"
        @testset "copy" begin
            geoA = zero(GeoArray(gdalstack[:a]))
            copy!(geoA, gdalstack, :a)
            # First wrap with GeoArray() here or == loads from disk for each cell.
            # we need a general way of avoiding this in all disk-based sources
            @test all(geoA .== GeoArray(gdalstack[:a]))
        end
    end

    @testset "save" begin
        geoA = GeoArray(gdalstack[:a])
        filename = tempname() * ".tif"
        write(filename, gdalstack)
        base, ext = splitext(filename)
        filename_b = string(base, "_b", ext)
        saved = GeoArray(GDALarray(filename_b))
        @test all(saved .== geoA)
    end

    @testset "show" begin
        sh = sprint(show, MIME("text/plain"), gdalstack)
        # Test but don't lock this down too much
        @test occursin("DiskStack", sh)
        @test occursin("GDALarray", sh)
        @test occursin("Y", sh)
        @test occursin("X", sh)
        @test occursin("Band", sh)
        @test occursin(":a", sh)
        @test occursin(":b", sh)
        @test occursin("cea.tif", sh)
    end

end

@testset "GDAL series" begin
    gdalser = series([path, path], (Ti(),); childkwargs=(mappedcrs=EPSG(4326), name=:test))
    @test GeoArray(gdalser[Ti(1)]) == GeoArray(GDALarray(path; mappedcrs=EPSG(4326), name=:test))

    gdalstack = GDALstack((a=path, b=path); childtype=GDALarray, childkwargs=(mappedcrs=EPSG(4326),))
    gdalser = GeoSeries([gdalstack, gdalstack], (Ti,))
    @test gdalser[1].childkwargs == gdalstack.childkwargs
    # Rebuild the ser by wrapping the GDALarray data in Array.
    # `modify` forces `rebuild` on all containers as in-Memory variants
    modified_ser = modify(Array, gdalser)
    @test typeof(modified_ser) <: GeoSeries{<:GeoStack{<:NamedTuple{(:a,:b),<:Tuple{<:GeoArray{UInt8,3,<:Tuple,<:Tuple,<:Array{UInt8,3}},Vararg}}}}

    @testset "read" begin
        geoseries = read(gdalser)
        @test geoseries isa GeoSeries{<:GeoStack}
        @test geoseries.data isa Vector{<:GeoStack}
        @test first(geoseries.data[1].data) isa GeoArray
    end
end

nothing