text:
import ClimateMachine.BalanceLaws:
# declaration
    vars_state,
# initialization
    init_state_prognostic!,
    init_state_auxiliary!,
    nodal_init_state_auxiliary!,
# rhs computation
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
# boundary conditions
    boundary_conditions,
    boundary_state!

import ClimateMachine.DGMethods.NumericalFluxes:
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder,
    NumericalFluxFirstOrder,
    NumericalFluxSecondOrder,
    RusanovNumericalFlux,
    numerical_boundary_flux_second_order!,
    numerical_flux_second_order!,
    numerical_boundary_flux_first_order!,
    numerical_flux_first_order!,
    normal_boundary_flux_second_order!

using Unitful

##
#     Declaration of variables
##
"""
    vars_state
    - returns a NamedTuple of data types.
"""
function vars_state(model::SlabOceanModelSetup, aux::Auxiliary, T)

    @vars begin
        x::T
        y::T
        z::T
        F_ρe_prescribed::T # stores prescribed flux for secondary (land import)
    end
end

function vars_state(::SlabOceanModelSetup, ::Prognostic, T)
    @vars begin
        T_sfc::T
    end
end

vars_state(::SlabOceanModelSetup, ::Gradient, T) = @vars()

vars_state(::SlabOceanModelSetup, ::GradientFlux, T) = @vars()

##
#     Initialization of variables
##
"""
    init_state_xyz!
    - sets up the initial fields within our state variables
    (e.g., prognostic, auxiliary, etc.), however it seems to not initialized
    the gradient flux variables by default.
"""

function nodal_init_state_auxiliary!(m::SlabOceanModelSetup, state_auxiliary, tmp, geom)
    init_state_auxiliary!(m, m.physics.orientation, state_auxiliary, geom)

    state_auxiliary.F_ρe_prescribed = 0
end

function init_state_auxiliary!(::SlabOceanModelSetup, ::SphericalOrientation, state_auxiliary, geom)
    FT = eltype(state_auxiliary)
    r = norm(geom.coord)
    state_auxiliary.x = geom.coord[1]
    state_auxiliary.y = geom.coord[2]
    state_auxiliary.z = geom.coord[3]

end

function init_state_auxiliary!(::SlabOceanModelSetup, ::Union{NoOrientation, FlatOrientation}, state_auxiliary, geom)
    FT = eltype(state_auxiliary)
    @inbounds r = geom.coord[3]
    state_auxiliary.x = geom.coord[1]
    state_auxiliary.y = geom.coord[2]
    state_auxiliary.z = geom.coord[3]
end

function init_state_prognostic!(model::SlabOceanModelSetup, state::Vars, aux::Vars, localgeo, t)
    x = aux.x
    y = aux.y
    z = aux.z

    parameters = model.physics.parameters
    ic = model.initial_conditions

    FT = eltype(state)
    state.T_sfc = FT(0)
    state.T_sfc = ic.T_sfc(parameters, x, y, z)

    return nothing
end

"""
    LHS computations
"""
flux_first_order!(model::SlabOceanModelSetup, _...) = nothing
compute_gradient_argument!(model::SlabOceanModelSetup, _...) = nothing
compute_gradient_flux!(model::SlabOceanModelSetup, _...) = nothing
flux_second_order!(model::SlabOceanModelSetup, _...) = nothing

"""
    RHS computations
"""
@inline function source!(
    model::SlabOceanModelSetup,
    source::Vars,
    state::Vars,
    gradflux::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    x = aux.x
    y = aux.y
    z = aux.z
    p = model.physics.parameters
    FT = eltype(aux)

    Q = FT(0)
    if boundary_mask(p, x, y, z)
        Q = idealized_qflux(p, x, y, z)
    end

    source.T_sfc = (aux.F_ρe_prescribed + Q) / (p.ρ_o * p.c_o * p.h_o)

    return nothing
end

@inline function idealized_qflux(p, x...)
    """
    Calculates additional tropical fluxes for a more realistic representation of surface tropical heating in idealized experiments
    """
    FT = eltype(p.Q_0)
    ϕ = lat(x...)
    Q = p.Q_0 * (FT(1) - 2 * ϕ^2 / p.L_w^2) * exp(-ϕ^2 / p.L_w^2) / cos(ϕ)

    ϕ_Qcutoff = FT(85 / 180 * pi) # this is not in the GFLD model but here needed to avoid polar discontinuities
    Q = abs(ϕ) > ϕ_Qcutoff ? FT(0) : Q
end

"""
    Numerics and boundary conditions
    - no interior / exterior interface fluxes needed for the slab
"""
@inline boundary_conditions(model::SlabOceanModelSetup) = model.boundary_conditions

numerical_flux_first_order!(::Nothing, model::Union{SlabOceanModelSetup}, _...) = nothing
numerical_flux_gradient!(nf, bc, model::Union{SlabOceanModelSetup}, _...) = nothing
function numerical_flux_second_order!(::Nothing, model::Union{SlabOceanModelSetup}, fluxᵀn, _...)
    fluxᵀn.T_sfc = fluxᵀn.T_sfc .* Float64(0) # this overwrites the internal numerical flux contributions to the normal flux
end

numerical_boundary_flux_first_order!(::Nothing, bc, model::Union{SlabOceanModelSetup}, _...) = nothing
numerical_boundary_flux_gradient!(::Nothing, bc, model::Union{SlabOceanModelSetup}, _...) = nothing
function numerical_boundary_flux_second_order!(::Nothing, bc, model::Union{SlabOceanModelSetup}, fluxᵀn, _...)
    fluxᵀn.T_sfc = fluxᵀn.T_sfc .* Float64(0) # this overwrites the internal numerical flux contributions to the normal flux
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine

max_forks_repo_name: CliMA/CouplerMachine


text:
Test.@test( PredictMD.remove_all_full_stops(".1.2.3.4.") == "1234" )
Test.@test( PredictMD.remove_all_full_stops("1234") == "1234" )

PredictMD.parse_icd_icd9_ccs_appendixasingledx_file!()
PredictMD.parse_icd_icd9_ccs_appendixasingledx_file!()
PredictMD.parse_icd_icd9_ccs_appendixasingledx_file!()

Test.@test(PredictMD.single_level_dx_ccs_number_to_name(1) == "Tuberculosis")
Test.@test(PredictMD.single_level_dx_ccs_number_to_name(2) == "Septicemia (except in labor)")

Test.@test(
    all(sort(unique(
        PredictMD.single_level_dx_ccs_to_list_of_icd9_codes(107))) .==
            sort(unique(["42741", "42742", "4275"])))
    )

Test.@test(PredictMD.icd9_code_to_single_level_dx_ccs("42741") == 107)

ext: jl


text:
using Clustering
using Distances: SqEuclidean, pairwise
using Base.Test

srand(34568)

@assert RandSeedAlg <: SeedingAlgorithm
@assert KmppAlg <: SeedingAlgorithm
@assert KmCentralityAlg <: SeedingAlgorithm

alldistinct(x::Vector{Int}) = (length(Set(x)) == length(x))

function min_interdist(X::Matrix)
    dists = pairwise(SqEuclidean(), X)
    n = size(X, 2)
    r = Inf
    for i = 1:n, j = 1:n
        if i != j && dists[i,j] < r
            r = dists[i,j]
        end
    end
    return r
end

d = 3
n = 100
k = 5
X = rand(d, n)
C = pairwise(SqEuclidean(), X)

md0 = min_interdist(X)

## RandSeed

iseeds = initseeds(RandSeedAlg(), X, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

iseeds = initseeds_by_costs(RandSeedAlg(), C, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

R = copyseeds(X, iseeds)
@test isa(R, Matrix{Float64})
@test R == X[:, iseeds]

## Kmpp

iseeds = initseeds(KmppAlg(), X, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

iseeds = initseeds_by_costs(KmppAlg(), C, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

@test min_interdist(X[:, iseeds]) > 20 * md0

iseeds = kmpp(X, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

iseeds = kmpp_by_costs(C, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

## Kmcentrality

iseeds = initseeds(KmCentralityAlg(), X, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

iseeds = initseeds_by_costs(KmCentralityAlg(), C, k)
@test length(iseeds) == k
@test alldistinct(iseeds)

@test min_interdist(X[:, iseeds]) > 2 * md0


text:
## CONSTANTS

const Dist = Distributions

## EQUALITY OF DISTRIBUTIONS (minor type piracy)

# TODO: We should get rid of this. I think it is used only in
# MLJModels/test.

function ==(d1::D, d2::D) where D<:Dist.Sampleable
    ret = true
    for fld in fieldnames(D)
        ret = ret && getfield(d1, fld) == getfield(d2, fld)
    end
    return ret
end

## DISTRIBUTION AS TRAIT (needed?)

# fallback:
isdistribution(d) = isdistribution(typeof(d))
isdistribution(::Type{<:Any}) = false

# for anything sampleable in Distributions.jl:
isdistribution(::Type{<:Dist.Sampleable}) = true

## ADD TO distributions.jl TYPE HIERARCHY TO ACCOUNT FOR NON-EUCLIDEAN
## SUPPORTS

abstract type NonEuclidean <: Distributions.ValueSupport end

## A NEW TRAIT FOR DISTRIBUTIONS

support_scitype(::Type) = Unknown
support_scitype(d) = support_scitype(typeof(d))

## UNIVARIATE NOMINAL PROBABILITY DISTRIBUTION

"""
    UnivariateFinite(classes, p)

A discrete univariate distribution whose finite support is the
elements of the vector `classes`, and whose corresponding
probabilities are elements of the vector `p`, which must sum to
one. Here `classes` must have type
`AbstractVector{<:CategoricalElement}` where

    CategoricalElement = Union{CategoricalValue,CategoricalString}

and all classes are assumed to share the same categorical pool.

    UnivariateFinite(prob_given_class)

A discrete univariate distribution whose finite support is the set of
keys of the provided dictionary, `prob_given_class`. The dictionary
keys must be of type `CategoricalElement` (see above) and the dictionary
values specify the corresponding probabilities.

    classes(d::UnivariateFinite)

A list of categorial elements in the common pool of classes used to
construct `d`.

    levels(d::UnivariateFinite)

A list of the raw levels in the common pool of classes used to
construct `d`, equal to `get.(classes(d))`.

    Distributions.support(d::UnivariateFinite)

Those classes associated with non-zero probabilities.

````julia
v = categorical(["yes", "maybe", "no", "yes"])
d = UnivariateFinite(v[1:2], [0.3, 0.7])
pdf(d, "yes")     # 0.3
pdf(d, v[1])      # 0.3
pdf(d, "no")      # 0.0
pdf(d, "house")   # throws error
classes(d) # CategoricalArray{String,1,UInt32}["maybe", "no", "yes"]
levels(d) # Array{String, 1}["maybe", "no", "yes"]
support(d) # CategoricalArray{String,1,UInt32}["maybe", "no"]
mode(d)    # CategoricalString{UInt32} "maybe"
rand(d, 5) # CategoricalArray{String,1,UInt32}["maybe", "no", "maybe", "maybe", "no"] or similar
d = fit(UnivariateFinite, v)
pdf(d, "maybe") # 0.25
```

One can also do weighted fits:

```julia
w = [1, 4, 5, 1] # some weights
d = fit(UnivariateFinite, v, w)
pdf(d, "maybe") ≈ 4/11 # true
````

*Warning:* The `pdf` function will give wrong answers if the order of
 levels of any categorical element passed to the UnivariateFinite
 constructor is changed.

See also `classes`, `support`.

"""
struct UnivariateFinite{L,U,T<:Real} <: Dist.Distribution{Dist.Univariate,NonEuclidean}
    decoder::CategoricalDecoder{L,U}
    prob_given_class::LittleDict{U,T}
end

UnivariateFinite(prob_given_class::AbstractDict) =
    throw(ArgumentError("The support of a UnivariateFinite "*
                        "can consist only of `CategoricalString` "*
                        "or `CategoricalValue` elements, and "*
                        "probabilities must be `AbstractFloat`. "))

function UnivariateFinite(prob_given_class::AbstractDict{L,T}) where {U<:Unsigned,L<:CategoricalElement{U},T<:Real}

    an_element = first(keys(prob_given_class))
    decoder_ = decoder(an_element)

    p = values(prob_given_class) |> collect
    Dist.@check_args(UnivariateFinite, Dist.isprobvec(p))

    d = LittleDict{U,T}()
    for key in classes(an_element)
        haskey(prob_given_class, key) && (d[int(key)] = prob_given_class[key] )
    end
    return UnivariateFinite(decoder_, d)
end

function UnivariateFinite(classes::AbstractVector{<:CategoricalElement},
                          p::AbstractVector{<:Real})
    Dist.@check_args(UnivariateFinite, length(classes)==length(p))
    prob_given_class = LittleDict([classes[i]=>p[i] for i in eachindex(p)])
    return  UnivariateFinite(prob_given_class)
end
UnivariateFinite(classes::AbstractVector, p) =
    throw(ArgumentError("`classes` must have type `AbstractVector{T}` where "*
                        "`T <: Union{CategoricalValue,CategoricalString}. "*
                        "Perhaps you have `T=Any`? "))

classes(d::UnivariateFinite) = classes(d.decoder.pool)
levels(d::UnivariateFinite)  = d.decoder.pool.levels

# get the internal integer representations of the support
raw_support(d::UnivariateFinite) =
    sort!(collect(keys(d.prob_given_class)))

Distributions.support(d::UnivariateFinite) =
    map(d.decoder, raw_support(d))

function Base.show(stream::IO, d::UnivariateFinite)
    raw = sort!(raw_support(d)) # reflects order of pool at
                                # instantiation of d
    x1 = d.decoder(first(raw))
    p1 = d.prob_given_class[first(raw)]
    str = "UnivariateFinite($x1=>$(round(p1, sigdigits=3))"
    pairs = (d.decoder(r)=>d.prob_given_class[r] for r in raw[2:end])
    for pair in pairs
        str *= ", $(pair[1])=>$(round(pair[2], sigdigits=3))"
    end
    str *= ")"
    print(stream, str)
end

support_scitype(::Type{<:UnivariateFinite}) = Finite

"""
    isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)

Returns `true` if and only if `Set(classes(d1) == Set(classes(d2))`
and the corresponding probabilities are approximately equal. The
key-word arguments `kwargs` are passed through to each call of
`isapprox` on probabiliity pairs. Returns `false` otherwise.

"""
function Base.isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)

    classes1 = classes(d1)
    classes2 = classes(d2)

    for c in classes1
        c in classes2 || return false
        isapprox(pdf(d1, c), pdf(d2, c); kwargs...) ||
            return false # pdf defined below
    end
    return true
end

function average(dvec::AbstractVector{UnivariateFinite{L,U,T}};
                 weights=nothing) where {L,U,T}

    n = length(dvec)

    Dist.@check_args(UnivariateFinite, weights == nothing || n==length(weights))

    # check all distributions have consistent pool:
    first_index = first(dvec).decoder.pool.index
    for d in dvec
        d.decoder.pool.index == first_index ||
            error("Averaging UnivariateFinite distributions with incompatible"*
                  " pools. ")
    end

    # get all refs:
    refs = Tuple(reduce(union, [keys(d.prob_given_class) for d in dvec]))

    # initialize the prob dictionary for the distribution sum:
    prob_given_class = LittleDict{U,T}(refs, zeros(T, length(refs)))

    # make vector of all the distributions dicts padded to have same common keys:
    prob_given_class_vec = map(dvec) do d
        merge(prob_given_class, d.prob_given_class)
    end

    # sum up:
    if weights == nothing
        scale = 1/n
        for x in refs
            for k in 1:n
                prob_given_class[x] += scale*prob_given_class_vec[k][x]
            end
        end
    else
        scale = 1/sum(weights)
        for x in refs
            for k in 1:n
                prob_given_class[x] +=
                    weights[k]*prob_given_class_vec[k][x]*scale
            end
        end
    end

    return UnivariateFinite(first(dvec).decoder, prob_given_class)
end

function _pdf(d::UnivariateFinite{L,U,T}, ref) where {L,U,T}
    return get(d.prob_given_class, ref, zero(T))
end

Distributions.pdf(d::UnivariateFinite{L,U,T},
                  x::CategoricalElement) where {L,U,T} = _pdf(d, int(x))

function Distributions.pdf(d::UnivariateFinite{L,U,T},
                           level::L) where {L,U,T}
    p = d.decoder.pool
    if level in p.levels
        return _pdf(d, int(p, level))
    else
        throw(ArgumentError(""))
    end
end

function Distributions.mode(d::UnivariateFinite)
    dic = d.prob_given_class
    p = values(dic)
    max_prob = maximum(p)
    m = first(first(dic)) # mode, just some ref for now
    for (x, prob) in dic
        if prob == max_prob
            m = x
            break
        end
    end
    return d.decoder(m)
end

"""
    _cummulative(d::UnivariateFinite)

Return the cummulative probability vector `[0, ..., 1]` for the
distribution `d`, using whatever ordering is used in the dictionary
`d.prob_given_class`. Used only for to implement random sampling from
`d`.

"""
function _cummulative(d::UnivariateFinite{L,U,T}) where {L,U,T<:Real}
    p = collect(values(d.prob_given_class))
    K = length(p)
    p_cummulative = Array{T}(undef, K + 1)
    p_cummulative[1] = zero(T)
    p_cummulative[K + 1] = one(T)
    for i in 2:K
        p_cummulative[i] = p_cummulative[i-1] + p[i-1]
    end
    return p_cummulative
end

"""
_rand(p_cummulative)

Randomly sample the distribution with discrete support `1:n` which has
cummulative probability vector `p_cummulative=[0, ..., 1]` (of length
`n+1`). Does not check the first and last elements of `p_cummulative`
but does not use them either.

"""
function _rand(p_cummulative)
    real_sample = rand()
    K = length(p_cummulative)
    index = K
    for i in 2:K
        if real_sample < p_cummulative[i]
            index = i - 1
            break
        end
    end
    return index
end

function Base.rand(d::UnivariateFinite)
    p_cummulative = _cummulative(d)
    classes = d.decoder.pool.valindex[collect(keys(d.prob_given_class))]
    return classes[_rand(p_cummulative)]
end

function Base.rand(d::UnivariateFinite, n::Int)
    p_cummulative = _cummulative(d)
    classes = d.decoder.pool.valindex[collect(keys(d.prob_given_class))]
    return [classes[_rand(p_cummulative)] for i in 1:n]
end

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{L}) where L
    L <: CategoricalElement ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "CategoricalValue or CategoricalString type. ")
    y = skipmissing(v) |> collect
    isempty(y) && error("No non-missing data to fit. ")
    N = length(y)
    count_given_class = Dist.countmap(y)
    classes = Tuple(keys(count_given_class))
    probs = values(count_given_class)./N
    prob_given_class = LittleDict(classes, probs)
    return UnivariateFinite(prob_given_class)
end

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{L},
                           weights::AbstractVector{<:Real}) where L
    L <: CategoricalElement ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "CategoricalValue or CategoricalString type. ")
    y = broadcast(identity, skipmissing(v))
    isempty(y) && error("No non-missing data to fit. ")
    classes_seen = filter(in(unique(y)), classes(y[1]))

    # instantiate and initialize prob dictionary:
    prob_given_class = LittleDict{L,Float64}()
    for c in classes_seen
        prob_given_class[c] = 0
    end

    # compute unnormalized  probablilities:
    for i in eachindex(y)
        prob_given_class[y[i]] += weights[i]
    end

    # normalize the probabilities:
    S = sum(values(prob_given_class))
    for c in keys(prob_given_class)
        prob_given_class[c] /=S
    end

    return UnivariateFinite(prob_given_class)

end

## NORMAL DISTRIBUTION ARITHMETIC

# *(lambda::Number, d::Dist.Normal) = Dist.Normal(lambda*d.μ, lambda*d.σ)
# function +(ds::Dist.Normal...)
#     μ = sum([d.μ for d in ds])
#     σ = sum([d.σ^2 for d in ds]) |> sqrt
#     return Dist.Normal(μ, σ)
# end

# hack for issue #809 in Dist.jl:
Distributions.fit(::Type{Dist.Normal{T}}, args...) where T =
    Dist.fit(Dist.Normal, args...)


text:
# test_threading.jl

@testset "Test threading" begin
    con = DBInterface.connect(DuckDB.DB)

    DBInterface.execute(con, "CREATE TABLE integers AS SELECT * FROM range(100000000) t(i)")
    results = DBInterface.execute(con, "SELECT SUM(i) sum FROM integers")
    df = DataFrame(results)
    @test df.sum == [4999999950000000]

    DBInterface.close!(con)
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
function rrule(::typeof(dot), t1::AbstractTensorMap{ComplexSpace}, t2::AbstractTensorMap{ComplexSpace})
    fwd = dot(t1, t2)

    function dot_pushback(f̄wd)
        t̄1 = f̄wd' * t2
        t̄2 = f̄wd * t1
        return NoTangent(), t̄1, t̄2
    end
    return fwd, dot_pushback
end

function toarray(t::TensorMap{ComplexSpace})
    @inline tuplejoin(x, y) = (x..., y...)
    t_dims = tuplejoin(dims(codomain(t)), dims(domain(t)))
    return reshape(t.data, t_dims)
end

@inline get_chi(psi::TensorMap{ComplexSpace, 2, 1}) = dim(domain(psi))
@inline get_d(psi::TensorMap{ComplexSpace, 2, 1}) = dim(codomain(psi)) ÷ dim(domain(psi))

@non_differentiable get_chi(psi::TensorMap{ComplexSpace, 2, 1})
@non_differentiable get_d(psi::TensorMap{ComplexSpace, 2, 1})

function arr_to_TensorMap(arr::Array{ComplexF64, 3})
    chi, d, _ = size(arr)
    return TensorMap(arr, ℂ^chi*ℂ^d, ℂ^chi)
end
function rrule(::typeof(arr_to_TensorMap), arr::Array{ComplexF64, 3})
    chi, d, _ = size(arr)
    fwd = TensorMap(arr, ℂ^chi*ℂ^d, ℂ^chi)
    function arr_to_TensorMap_pushback(f̄wd)
        return NoTangent(), reshape(f̄wd.data, (chi, d, chi))
    end
    return fwd, arr_to_TensorMap_pushback
end

*(x::ComplexSpace) = x

"""
    convert_to_tensormap(arr::Array{<:Complex}, d_dom::Integer) -> TensorMap

    Given the dimension of the domain, convert an array to the tensormap.
    FIXME: fix vector cases.
"""
function convert_to_tensormap(arr::Array{<:Complex}, d_dom::Integer)
    size_arr = size(arr)
    d_arr = length(size_arr)
    dom = ifelse(d_dom > 0, prod(map(x->ℂ^x, size_arr[1:d_dom])), ProductSpace{ComplexSpace, 0})
    codom = ifelse(d_arr > d_dom, prod(map(x->ℂ^x, size_arr[d_dom+1:end])), ProductSpace{ComplexSpace, 0})
    return TensorMap(arr, dom, codom)
end

"""
    rrule(typeof(convert_to_tensormap), arr::Array{<:Complex}, d_dom::Integer)

    The rrule for function `convert_to_tensormap`.
"""
function rrule(::typeof(convert_to_tensormap), arr::Array{<:Complex}, d_dom::Integer)
    fwd = convert_to_tensormap(arr, d_dom)

    function convert_to_tensormap_pushback(f̄wd)
        return NoTangent(), convert(Array, f̄wd), NoTangent()
    end
    return fwd, convert_to_tensormap_pushback
end

"""
    quicksave(name::String, psi::TensorMap{ComplexSpace, 2, 1})

    Save the mps local tensor `psi` to a `.jld2` file.
"""
function quicksave(name::String, psi::TensorMap{ComplexSpace, 2, 1})
    # TODO. get convert(Dict, t::AbstractTensorMap) in the manual to work with JLD2?
    chi, d = get_chi(psi), get_d(psi)
    psi_dict = Dict("chi" => chi,
                    "d" => d,
                    "data" => psi.data
    )

    save("$(name).jld2", psi_dict)
end

"""
    quickload(name::String)

    load the mps local tensor `psi` from a `.jld2` file.
"""
function quickload(name::String)
    psi_dict = load("$(name).jld2")
    chi, d = psi_dict["chi"], psi_dict["d"]
    return TensorMap(psi_dict["data"], ℂ^chi*ℂ^d, ℂ^chi)
end

"""
    logsumexp(w)

    Calculation of logsumexp.
"""
function logsumexp(w::Array{<:Number, 1})
    u = maximum(real.(w))
    return u + log(sum(exp.(w .- u)))
end

"""
    rrule(::typeof(logsumexp), w::Array{<:Number, 1})

    Backwards rule for logsumexp.
    I did this because it seems Zygote's automatic backwards will fail at `maximum(norm.(w))` in the case of complex numbers.
"""
function rrule(::typeof(logsumexp), w::Array{<:Number, 1})
    u = maximum(real.(w))
    expw = exp.(w .- u)
    fwd = u + log(sum(expw))

    function logsumexp_pushback(f̄wd)
        w̄ = f̄wd * conj(expw / sum(expw))
        return NoTangent(), w̄
    end
    return fwd, logsumexp_pushback
end

# FIXME. seems buggy
function tmp_eig(t::TensorMap{ComplexSpace})
    tarr = convert_to_array(t)
    dom = domain(t)
    Wdim = dom |> dims |> prod
    tarr = reshape(tarr, (Wdim, Wdim))

    Warr, Varr = eigen(tarr)
    Warr = diagm(Warr)

    W = convert_to_tensormap(Warr, 1)
    V = TensorMap(Varr, dom, ℂ^Wdim)

    return W, V
end

# FIXME. seems buggy
function tmp_inv(t::TensorMap{ComplexSpace})
    tarr = convert_to_array(t)
    dom = domain(t)
    codom = codomain(t)
    tdim = dom |> dims |> prod
    tarr = reshape(tarr, (tdim, tdim))

    tinv_arr = inv(tarr)
    tinv = TensorMap(tinv_arr, dom, codom)
    return tinv
end


text:
@testset "Distances" begin

    dgm1 = diagram(0.5=>1.0)
    dgm2 = diagram(0.5=>1.1)
    dgm3 = diagram(0.5=>1.0, 0.6=>1.1, 0.3 => Inf)

    m = convert(Matrix, dgm1)
    @test size(m) == (2,1)
    m = convert(Matrix, dgm3)
    @test size(m) == (2,2)
    m = convert(Matrix, dgm3, skipinf=false)
    @test size(m) == (2,3)

    @test wasserstein(dgm1, dgm1) ≈ 0.0
    @test wasserstein(dgm1, dgm2) ≈ 0.1
    @test wasserstein(dgm2, dgm2) ≈ 0.0
    @test wasserstein(dgm3, dgm2) ≈ 0.2
    @test wasserstein(dgm2, dgm3) ≈ 0.2
    @test wasserstein(diagram(1.0=>2.0), Interval{Float64}[]) ≈ √2/2
    @test wasserstein(diagram(1.0=>2.0), diagram(0.0=>Inf)) ≈ √2/2
    @test wasserstein(diagram(0.0=>Inf), diagram(1.0=>2.0)) ≈ √2/2

end


text:
for T in subtypes(PolyElem)
  (f::T)(a) = subst(f, a)

  function (f::T)(a::T)
     if parent(f) != parent(a)
        return subst(f, a)
     end
     return compose(f, a)
  end

  (f::T)(a::Integer) = evaluate(f, a)
  (f::T)(a::fmpz) = evaluate(f, a)

  function (f::T)(a::RingElem)
     if parent(a) != base_ring(f)
        return subst(f, a)
     end
     return evaluate(f, a)
  end

end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.QML

using LinearAlgebra: I
using QuantumCircuits.QCircuits.Gates: Rx, Ry, Rz
using QuantumCircuits.QCircuits.Math

################################################################################
#  Get save, parameters                                                        #
################################################################################
qc = QCircuit(1)
qc.u3(0, ParameterT(1.0), ParameterT(2.0), ParameterT(3.0))
qc.u3(0, ParameterT(4.0), ParameterT(5.0), ParameterT(6.0))
@test getparameters(qc) == [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]

setparameters!(qc, [1.5, 2.5, 3.5, 4.5, 5.5, 6.5])
@test getparameters(qc) == [1.5, 2.5, 3.5, 4.5, 5.5, 6.5]

################################################################################
# Rotations
gates = [Rx, Ry, Rz]
θs = [0, π/4, π/2, π]

for g in gates
    for θ in θs
        circ = QCircuit(1)
        add!(circ, g, 0, ParameterT(θ))

        @test getparameters(circ) == [θ]
        setparameters!(circ, [1.0])
        @test getparameters(circ) == [1.0]
    end
end

################################################################################
#  Compare                                                                     #
################################################################################
qc1 = QCircuit(2)
qc1.x(0)
qc1.h(1)
qc1.cx(0, 1)

qc2 = QCircuit(2)
qc2.h(1)
qc2.x(0)
qc2.cx(0, 1)

@test qc1 == qc2

qr11 = QuantumRegister(1)
qr12 = QuantumRegister(1)
qc1 = QCircuit([qr11, qr12])
qc1.x(qr11[0])
qc1.h(qr12[0])

qc2 = QCircuit(2)
qc2.h(1)
qc2.x(0)

@test qc1 == qc2

################################################################################
function check_inv(qc)
    err = matrix_norm(tomatrix(inv(qc)) * tomatrix(qc) - I)
    @test abs(err) < 1e-16
end

# Inverse
qc = QCircuit(1)
qc.x(0)
qc.h(0)
check_inv(qc)

qc = QCircuit(1)
qc.y(0)
qc.z(0)
qc.x(0)
qc.h(0)
check_inv(qc)

qc = QCircuit(1)
qc.rx(0, π/2)
check_inv(qc)

qc = QCircuit(1)
qc.ry(0, π/2)
check_inv(qc)

qc = QCircuit(1)
qc.rz(0, π/2)
check_inv(qc)

qc = QCircuit(2)
qc.ry(0, π)
qc.cx(0, 1)
qc.rz(1, π)
qc.rx(0, π/2)
check_inv(qc)

qc = QCircuit(2)
qc.s(0)
qc.cx(1, 0)
qc.t(1)
check_inv(qc)

qc = QCircuit(2)
qc.sdg(0)
qc.cx(1, 0)
qc.tdg(1)
check_inv(qc)

qc = QCircuit(1)
qc.u3(0, π/4, π/2, π)
check_inv(qc)

qc = QCircuit(2)
qc.u3(0, π/4, π/2, π)
qc.cx(1, 0)
qc.u3(1, π/2, π/4, π/8)
check_inv(qc)

################################################################################
#  Append                                                                     #
################################################################################
qc1 = QCircuit(2)
qc1.x(0)
qc1.h(1)
qc1.cx(0, 1)
qc1.u3(1, 1, 2, 3)

qc2 = QCircuit(2)
qc2.cx(0, 1)
qc2.rx(0, 4.0)
qc2.ry(1, 5.0)
qc2.rz(1, 6.0)

testqc = QCircuit(2)
testqc.x(0)
testqc.h(1)
testqc.cx(0, 1)
testqc.u3(1, 1, 2, 3)
testqc.cx(0, 1)
testqc.rx(0, 4.0)
testqc.ry(1, 5.0)
testqc.rz(1, 6.0)

append!(qc1, qc2)

@test qc1 == testqc


text:
module Util

include("similarity_FM.jl")

end # Util

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
#=
REST API Specification
----------------------

    • Malformed GET, POST requests return a HTTP 400 i.e. Bad Request
    • Correct resutls return a HTTP 200 i.e. OK

    • GET:  /*                    returns a HTTP 501 message i.e. Not Implemented
    • GET:  /api/kill             kills the server
    • GET:  /api/read-configs     returns the searcher configurations

    • POST: /*                    returns a HTTP 501 message i.e. Not Implemented
    • POST: /api/envop            loads (deserializes)/ saves (serializes)/ re-indexes the search environment
    • POST: /api/search           triggers a search using parameters from the HTTP message body
    • POST: /api/recommed         triggers a recommendation using parameters from the HTTP message body
    • POST: /api/rank             triggers a ranking i.e. boost using parameters from the HTTP message body

HTTP Message body specification for the search, recommend, ranking and environment operations
---------------------------------------------------------------------------------------
    • /api/envop body format (JSON):
         {
          "cmd" : <the environment operation to be performed, a string; available: 'load', 'save' and 'reindex'>,
          "cmd_argument": <argument for the operation, a string; for 'load' and 'save' it should be a filepath,
                           for 'reindex' it should be either a searcher id or '*' i.e. all searchers>
         }
    • /api/search body format (JSON):
         {
          "query" : <the query to be performed, a string>,
          "input_parser": <the input parser to use; available: 'noop_input_parser', 'base_input_parser'>
          "return_fields" : <a list of string names for the fields to be returned>,
          "sort_fields" : <OPTIONAL, a list of string names for the fields to sort by when filtering. Sort precedence is given by list order>,
          "sort_reverse" : <OPTIONAL, whether to reverse the sorting i.e. largest number/letter first>,
          "search_method" : <OPTIONAL, a string defining the type of classic search method>,
          "searchable_filters" : <OPTIONAL, a list of fields whose values will also be part of search if used for filtering>
          "max_matches" : <OPTIONAL, an integer defining the maximum number of results for search>,
          "response_size" : <OPTIONAL, an integer defining the maximum number of results to be actually returned>,
          "response_page" : <OPTIONAL, an integer that specifies which page of 'response_size' results to return>,
          "max_suggestions" : <OPTIONAL, an integer defining the maximum number of suggestions / mismatches keyword>,
          "custom_weights" : <OPTIONAL, a dictionary where the keys are strings with searcher ids and values
                              the weights of the searchers in the result aggregation>
          "ranker": <OPTIONAL, the ranked to use; available: 'noop_ranker'>
         }

    • /api/recommend body format (JSON):
         {
          "recommender": <a string with the name of the recommender to use; available: 'noop_recommender', 'search_recommender'>
          "recommend_id" : <the id of the entity for which recommendations are sought>
          "recommend_id_key": <the db name of the column holding the id value>
          "input_parser": <the input parser to use; available: 'noop_input_parser', 'base_input_parser'>
          "filter_fields" : <a list of string name fields containing the fields that will be used by the recommender>,
          "return_fields" : <a list of string names for the fields to be returned>,
          "sort_fields" : <OPTIONAL, a list of string names for the fields to sort by when filtering. Sort precedence is given by list order>,
          "sort_reverse" : <OPTIONAL, whether to reverse the sorting i.e. largest number/letter first>,
          "search_method" : <OPTIONAL, a string defining the type of classic search method>,
          "searchable_filters" : <OPTIONAL, a list of fields whose values will form a search query if used in filter_fields>,
          "max_matches" : <OPTIONAL, an integer defining the maximum number of results for recommendations>,
          "response_size" : <OPTIONAL, an integer defining the maximum number of results to be actually returned>,
          "response_page" : <OPTIONAL, an integer that specifies which page of 'response_size' results to return>,
          "max_suggestions" : <OPTIONAL, an integer defining the maximum number of suggestions / mismatches keyword>,
          "custom_weights" : <OPTIONAL, a dictionary where the keys are strings with searcher ids and values
                              the weights of the searchers in the result aggregation>
          "ranker": <OPTIONAL, the ranked to use; available: 'noop_ranker'>
         }

    • /api/rank body format (JSON):
         {
            "ranker": <a string with the name of the ranker to use; available: 'noop_ranker'>
            "rank_ids": <list of ids to rank>,
            "rank_id_key": <the db name of the column holding the id values>,
            "return_fields" : <a list of string names for the fields to be returned>,
            "response_size" : <OPTIONAL, an integer defining the maximum number of results to be actually returned>,
            "response_page" : <OPTIONAL, an integer that specifies which page of 'response_size' results to return>
         }
=#

"""
    rest_server(port::Integer, io_port::Integer, search_server_ready::Condition [;ipaddr::String])

Starts a bi-directional HTTP REST server at address `ipaddr::String`
(defaults to `"0.0.0.0"` i.e. all ip's) that uses the TCP port `port`
and communicates with the search server through the TCP port `io_port`.
The server is started once the condition `search_server_ready`
is triggered.
"""
function rest_server(port::Integer,
                     io_port::Integer,
                     search_server_ready::Condition;
                     ipaddr::String="0.0.0.0")
    _HEADERS = ["Access-Control-Allow-Origin"=>"*"]
    #Checks
    if port <= 0
        @error "Please specify a HTTP REST port of positive integer value."
    end

    # Small closure that prints debug information for handlers
    __debug_wrapper(handle) = (req::HTTP.Request)->begin
        @debug "HTTP request $(req.target) received."
        handle(req)
    end

    # Define REST endpoints to dispatch to "service" functions
    GARAMOND_REST_ROUTER = HTTP.Router()
    HTTP.@register(GARAMOND_REST_ROUTER, "GET", "/*", __debug_wrapper(noop_req_handler))
    HTTP.@register(GARAMOND_REST_ROUTER, "GET", "/api/kill", __debug_wrapper(kill_req_handler))
    HTTP.@register(GARAMOND_REST_ROUTER, "GET", "/api/read-configs", __debug_wrapper(read_configs_req_handler))
    HTTP.@register(GARAMOND_REST_ROUTER, "POST", "/*", __debug_wrapper(noop_req_handler))
    HTTP.@register(GARAMOND_REST_ROUTER, "POST", "/api/envop", __debug_wrapper(envop_req_handler))
    HTTP.@register(GARAMOND_REST_ROUTER, "POST", "/api/search", __debug_wrapper(search_req_handler))
    HTTP.@register(GARAMOND_REST_ROUTER, "POST", "/api/recommend", __debug_wrapper(recommend_req_handler))
    HTTP.@register(GARAMOND_REST_ROUTER, "POST", "/api/rank", __debug_wrapper(rank_req_handler))

    # Wait for search server to be ready
    wait(search_server_ready)
    @info "REST server online @$ipaddr:$port..."

    # Start serving requests
    @async HTTP.serve(Sockets.IPv4(ipaddr), port, readtimeout=0) do http_req::HTTP.Request
        srchsrv_req = try
            HTTP.Handlers.handle(GARAMOND_REST_ROUTER, http_req)
        catch e
            @debug "Error handling HTTP request.\n$e"
            ERRORED_REQUEST
        end

        if srchsrv_req === UNINITIALIZED_REQUEST
            # An unsupported endpoint was called
            return HTTP.Response(501)
        elseif srchsrv_req === ERRORED_REQUEST
            # Something went wrong during handling
            return HTTP.Response(400)
        elseif srchsrv_req isa InternalRequest
            # All OK, send request to search server and get response
            ssconn = TCPSocket()
            try
                Sockets.connect!(ssconn, Sockets.localhost, io_port)
            catch
                return HTTP.Response(503)
            end
            println(ssconn, request2json(srchsrv_req))                   # writes a "\n" as well
            response = ifelse(isopen(ssconn), readline(ssconn), "")      # expects a "\n" terminator
            close(ssconn)
            return HTTP.Response(200, _HEADERS, body=response)
        else
            return HTTP.Response(400)
        end
    end
end

# Handlers (receive a HTTP request and return a search server request)
noop_req_handler(req::HTTP.Request) = UNINITIALIZED_REQUEST

kill_req_handler(req::HTTP.Request) = KILL_REQUEST

read_configs_req_handler(req::HTTP.Request) = READCONFIGS_REQUEST

envop_req_handler(req::HTTP.Request) = begin
    return InternalRequest(operation=:envop, query=httpbody2string(req))
end

search_req_handler(req::HTTP.Request) = begin
    parameters = httpbody2dict(req)
    return InternalRequest(
                operation = :search,
                query = parameters["query"],  # if missing, throws
                input_parser = Symbol(parameters["input_parser"]),  # if missing, throws
                return_fields = Symbol.(parameters["return_fields"]),  # if missing, throws
                sort_fields = Symbol.(get(parameters, "sort_fields", DEFAULT_SORT_FIELDS)),
                sort_reverse = get(parameters, "sort_reverse", DEFAULT_SORT_REVERSE),
                search_method = Symbol(get(parameters, "search_method", DEFAULT_SEARCH_METHOD)),
                searchable_filters = Symbol.(get(parameters, "searchable_filters", String[])),
                max_matches = get(parameters, "max_matches", DEFAULT_MAX_MATCHES),
                response_size = get(parameters, "response_size", DEFAULT_RESPONSE_SIZE),
                response_page = get(parameters, "response_page", DEFAULT_RESPONSE_PAGE),
                max_suggestions = get(parameters, "max_suggestions", DEFAULT_MAX_SUGGESTIONS),
                custom_weights = get(parameters, "custom_weights", DEFAULT_CUSTOM_WEIGHTS),
                ranker = Symbol(get(parameters, "ranker", DEFAULT_RANKER_NAME)))
end

recommend_req_handler(req::HTTP.Request) = begin
    parameters = httpbody2dict(req)
    _query = parameters["recommend_id"] * " " * join(parameters["filter_fields"], " ")
    return InternalRequest(
                operation = :recommend,
                recommender = Symbol(parameters["recommender"]),  # if missing, throws
                request_id_key = Symbol.(parameters["recommend_id_key"]),  # if missing, throws
                query = _query,
                input_parser = Symbol(parameters["input_parser"]),  # if missing, throws
                return_fields = Symbol.(parameters["return_fields"]),  # if missing, throws
                sort_fields = Symbol.(get(parameters, "sort_fields", DEFAULT_SORT_FIELDS)),
                sort_reverse = get(parameters, "sort_reverse", DEFAULT_SORT_REVERSE),
                search_method = Symbol(get(parameters, "search_method", DEFAULT_SEARCH_METHOD)),
                searchable_filters = Symbol.(get(parameters, "searchable_filters", String[])),
                max_matches = get(parameters, "max_matches", DEFAULT_MAX_MATCHES),
                response_size = get(parameters, "response_size", DEFAULT_RESPONSE_SIZE),
                response_page = get(parameters, "response_page", DEFAULT_RESPONSE_PAGE),
                max_suggestions = get(parameters, "max_suggestions", DEFAULT_MAX_SUGGESTIONS),
                custom_weights = get(parameters, "custom_weights", DEFAULT_CUSTOM_WEIGHTS),
                ranker = Symbol(get(parameters, "ranker", DEFAULT_RANKER_NAME)))
end

rank_req_handler(req::HTTP.Request) = begin
    parameters = httpbody2dict(req)
    _all_ids = strip.(parameters["rank_ids"])  # if missing, throws
    return InternalRequest(
                operation = :rank,
                ranker = Symbol(parameters["ranker"]),  # if missing, throws
                query = join(_all_ids, " "),
                request_id_key = Symbol.(parameters["rank_id_key"]),   # if missing, throws
                return_fields = Symbol.(parameters["return_fields"]),  # if missing, throws
                response_size = get(parameters, "response_size", length(_all_ids)),
                response_page = get(parameters, "response_page", DEFAULT_RESPONSE_PAGE))
end

httpbody2string(req::HTTP.Request) = read(IOBuffer(HTTP.payload(req)), String)
httpbody2dict(req::HTTP.Request) = JSON.parse(httpbody2string(req))


text:
#=
parseBall:
- Julia version:
- Author: Dan
- Date: 2021-11-28
=#

push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL

pdb_filepath = "1crn.pdb"

System(pdb_filepath)

elapsed = @elapsed System(pdb_filepath)

println(elapsed)


text:
if !success(`iconv --version`)
	warn("iconv not found, skipping unicode tests!")
	@windows_only warn("Use WinRPM.install(\"win_iconv\") to run these tests")
else
	# Create unicode test data directory
	unicodedir = mktempdir()

	# Use perl to generate the primary data
	primary_encoding = "UTF-32BE"
	primary_path = replace(joinpath(unicodedir, primary_encoding*".unicode"),"\\","\\\\\\\\")
	run(`perl -e "
	$$fname = \"$primary_path\";
	open(UNICODEF, \">\", \"$$fname\")         or die \"can\'t open $$fname: $$!\";
	binmode(UNICODEF);
	print UNICODEF pack \"N*\", 0xfeff, 0..0xd7ff, 0xe000..0x10ffff;
	close(UNICODEF);"` )

	# Use iconv to generate the other data
	for encoding in ["UTF-32LE", "UTF-16BE", "UTF-16LE", "UTF-8"]
	    output_path = joinpath(unicodedir, encoding*".unicode")
	    f = Base.FS.open(output_path,Base.JL_O_WRONLY|Base.JL_O_CREAT,Base.S_IRUSR | Base.S_IWUSR | Base.S_IRGRP | Base.S_IROTH)
	    run(`iconv -f $primary_encoding -t $encoding $primary_path` |> f)
	    Base.FS.close(f)
	end

	f=open(joinpath(unicodedir,"UTF-32LE.unicode"))
	str1 = UTF32String(reinterpret(Char, read(f, Uint32, 1112065)[2:end]))
	close(f)

	f=open(joinpath(unicodedir,"UTF-8.unicode"))
	str2 = UTF8String(read(f, Uint8, 4382595)[4:end])
	close(f)
	@test str1 == str2

	str1 = "∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"
	str2 = UTF32String(
	    8704, 32, 949, 32, 62, 32, 48, 44, 32, 8707, 32,
	    948, 32, 62, 32, 48, 58, 32, 124, 120, 45, 121, 124,
	    32, 60, 32, 948, 32, 8658, 32, 124, 102, 40, 120,
	    41, 45, 102, 40, 121, 41, 124, 32, 60, 32, 949
	)
	@test str1 == str2

	# Cleanup unicode data
	for encoding in ["UTF-32BE", "UTF-32LE", "UTF-16BE", "UTF-16LE", "UTF-8"]
	    rm(joinpath(unicodedir,encoding*".unicode"))
	end
	rmdir(unicodedir)

	# UTF16
	u8 = "𝕥𝟶f𠂊"
	u16 = utf16(u8)
	@test length(u16.data) == 7
	@test length(u16) == 4
	@test utf8(u16) == u8
	@test collect(u8) == collect(u16)
        @test u16 == utf16(u16.data) == utf16(reinterpret(Uint8, u16.data))
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
function MidPointIntegrator()
    return MidPointIntegrator(())
end

function MidPointIntegrator(arg0::jdouble, arg1::jdouble, arg2::jint, arg3::jint)
    return MidPointIntegrator((jdouble, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function MidPointIntegrator(arg0::jint, arg1::jint)
    return MidPointIntegrator((jint, jint), arg0, arg1)
end


text:
### Plain Monte Carlo jobs do not use tasks

type PlainMCJob{M<:MCModel, S<:MCSampler, R<:MCRunner, T<:MCTuner} <: MCJob
  model::Vector{M}
  sampler::Vector{S}
  runner::Vector{R}
  tuner::Vector{T}
  heap::Vector{MCHeap}
  send::Function
  receive::Function
  reset::Function
  jobtype::Symbol
  dim::Int

  function PlainMCJob(m::Vector{M}, s::Vector{S}, r::Vector{R}, t::Vector{T})
    job::PlainMCJob = new()

    job.dim = length(m)
    @assert job.dim == length(s) == length(r) == length(t) "Number of models, samplers, runners and tuners not equal."

    job.model, job.sampler, job.runner, job.tuner = m, s, r, t
    job.heap = MCHeap[initialize_heap(m[i], s[i], r[i], t[i]) for i = 1:job.dim]
    job.send = identity
    job.receive = (i::Int)->iterate!(job.heap[i], m[i], s[i], r[i], t[i], identity)
    job.reset = (i::Int, x::Vector{Float64})->reset!(job.heap[i], x)
    job.jobtype = :plain

    job
  end
end

PlainMCJob{M<:MCModel, S<:MCSampler, R<:MCRunner, T<:MCTuner}(m::M, s::S, r::R, t::T) =
  PlainMCJob{M, S, R, T}([m], [s], [r], [t])


text:
# # Making an Image of a Black Hole

# `Comrade` has been designed to work with the EHT and ngEHT.
# In this tutorial we will show how to reproduce some of the results
# from [EHTC VI 2019](https://iopscience.iop.org/article/10.3847/2041-8213/ab1141).

# In EHTC VI, they considered fitting simple geometric models to the data
# to estimate the image size, shape, brightness profile etc of the black hole.
# In this page we will construct a similar model and fit it to the data in under
# 50 lines of code (sans comments). To start we load some packages we will need

using Comrade

# The next step is to load the data. For this we will use the publically
# available M 87 data which can be downloaded
# from [cyverse](https://datacommons.cyverse.org/browse/iplant/home/shared/commons_repo/curated/EHTC_FirstM87Results_Apr2019).
# For an introduction to data loading see [Loading Data into Comrade](@ref).

load_ehtim()
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "../assets/SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))
obs.add_scans()
# kill 0-baselines since we don't care about large scale flux and make scan-average data
obs = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true)
# grab data products we want to fit: log closure amplitudes and closure phases
dlcamp = extract_lcamp(obs)
dcphase = extract_cphase(obs)

# For this demo only consider closure products since these are invariant to station specific
# gain systematics. Given these data products we can then form our radio likelihood:

lklhd = RadioLikelihood(dlcamp, dcphase)

# The `lklhd` constructs a measure using `MeasureTheory.jl`. To evaluate the likelihood
# we need to pass it a type that implements the `Comrade` model interface which is described
# in [Model Interface](@ref).

# To finish the construction of our posterior we need to specify an image model and a prior.
# For the image model we will be using a modified `MRing`, which is a
# infinitely thin delta ring with an azimuthal structure given by a Fourier expansion.
# To give the ring some width we will convolve the ring with a gaussian, and add an
# additional gaussian to the image to model any non-ring flux. For the model a user
# must give a function that accepts a named tuple and return the constructed model:

function model(θ)
  (;radius, width, α, β, f, σG, τG, ξG, xG, yG) = θ
  ring = f*smoothed(stretched(MRing((α,), (β,)), radius, radius), width)
  g = (1-f)*shifted(rotated(stretched(Gaussian(), σG, σG*(1+τG)), ξG), xG, yG)
  return ring + g
end

# We now need to specify the priors for our model. The easiest way to do this is to
# specify a NamedTuple of distributions:

using Distributions
prior = (
          radius = Uniform(μas2rad(10.0), μas2rad(30.0)),
          width = Uniform(μas2rad(1.0), μas2rad(10.0)),
          α = Uniform(-0.5, 0.5),
          β = Uniform(-0.5, 0.5),
          f = Uniform(0.0, 1.0),
          σG = Uniform(μas2rad(1.0), μas2rad(40.0)),
          τG = Uniform(0.0, 0.75),
          ξG = Uniform(0.0, 1π),
          xG = Uniform(-μas2rad(80.0), μas2rad(80.0)),
          yG = Uniform(-μas2rad(80.0), μas2rad(80.0))
        )

# To form the posterior we now call

post = Posterior(lklhd, prior, model)

# This constructs a posterior density that can be evaluated by calling `logdensityof`.
# For example

logdensityof(post, (radius = μas2rad(20.0),
                  width = μas2rad(10.0),
                  α = 0.3,
                  β = 0.3,
                  f = 0.6,
                  σG = μas2rad(20.0),
                  τG = 0.1,
                  ξG = 0.5,
                  xG = 0.0,
                  yG = 0.0))

# We can now try to sample from our posterior `post` so that we can make probabilstic
# inferences about our data.

# Now this model is in **parameter** space. Often optimization and sampling algorithms
# want it in some modified space. For example, nested sampling algorithms want the
# parameters in the unit hypercube. To transform the posterior to the unit hypercube we
# can use the `ascube` function

cpost = ascube(post)

# If we want to flatten the parameter space and move to (-∞, ∞) support we can use the
# `asflat` function

fpost = asflat(post)

# These transformed posterior expect a vector of parameters. That is we can evaluate the
# transformed log density by calling

logdensityof(cpost, rand(dimension(cpost)))
logdensityof(fpost, randn(dimension(fpost)))

# Note that this automatically takes care of the jacobian in the parameter transformation.

# ## Sampling the posterior

# Our strategy here will be to use Hamiltonian Monte Carlo. However, to lower burn-in time
# we will first use an optimizer to find a reasonable starting location. Since this is a lower
# dimensional problem we will use BlackboxOptim or the GalacticBBO package

using ComradeGalactic
using GalacticBBO

ndim = dimension(fpost)
f = OptimizationFunction(fpost)
prob = OptimizationProblem(f, randn(ndim), nothing, lb=fill(-5.0, ndim), ub=fill(5.0, ndim))
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(); maxiters=50_000)

# The sol vector is in the transformed space, so first we need to transform back to parameter space

xopt = transform(fpost, sol)

# And we can also plot the map

using Plots
plot(model(xopt), title="MAP image", xlims=(-60.0,50.0), ylims=(-60.0,50.0))

# The main goal of `Comrade` is to explore the posterior of the model parameters.
# Currently the go to tool is [AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl).
# To sample from the posterior you can use the following:

using ComradeAHMC
chain, stats = sample(post, AHMC(metric=DiagEuclideanMetric(ndim)), 2000; nadapts=1000, init_params=xopt)

# That's it! To finish it up we can then plot some simple visual fit diagnostics.

# First to plot the image we call

plot(model(chain[end]), title="Random image", xlims=(-60.0,50.0), ylims=(-60.0,50.0))

# What about the mean image? Well let's grab 100 images from the chain
meanimg = mean(intensitymap.(model.(sample(chain, 100)), μas2rad(120.0), μas2rad(120.0), 128, 128))
plot(sqrt.(meanimg), title="Mean Image") #plot on a sqrt color scale to see the Gaussian

# That looks similar to the EHTC VI, and it took us no time at all!. To see how well the
# model is fitting the data we can plot the model and data products

plot(model(xopt), dlcamp)

# or even better the residuals

residual(model(xopt), dlcamp)

# These residuals are a little ratty which suggests our model is missing some emission.
# In fact, this model is slightly too simple to explain the data.
# Check out [EHTC VI 2019](https://iopscience.iop.org/article/10.3847/2041-8213/ab1141)
# for some ideas what features need to be added to the model to get a better fit!


text:
########################################################################
# Client for synchronous learning
########################################################################
mutable struct Client{T1<:Int64, T2<:Float64, T3<:Vector{T1}, T4<:Matrix{T2}, T5<:SparseMatrixCSC{T2, T1}, T6<:Flux.Chain}
    id::T1                                  # client index
    Xtrain::T5                              # training data
    Xtest::T5                               # test data
    num_classes::T1                         # number of classes
    num_clients::T1                         # number of clients
    num_epoches::T1                         # number of epoches
    batch_size::T1                          # number of batches
    learning_rate::T2                       # learning rate
    W::Union{T4, T6}                        # client model
    batch::T3                               # mini-batch
    grads::T4                               # gradient information
    back                                    # backward operator
    function Client(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Xtest::SparseMatrixCSC{Float64, Int64}, config::Dict{String, Union{Int64, Float64, String}})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        num_epoches = config["num_epoches"]
        batch_size = config["batch_size"]
        learning_rate = config["learning_rate"]
        dm = size(Xtrain, 1)
        if config["local_model"] == "linear"
            W = zeros(Float64, num_classes, dm)
        elseif config["local_model"] == "mlp"
            W = Chain(Dense(dm, 32, σ), Dense(32, num_classes))
        end
        back = undef
        batch = zeros(Int64, batch_size)
        grads = zeros(Float64, num_classes, batch_size)
        new{Int64, Float64, Vector{Int64}, Matrix{Float64}, SparseMatrixCSC{Float64, Int64},Flux.Chain}(id, Xtrain, Xtest, num_classes, num_clients, num_epoches, batch_size, learning_rate, W, batch, grads, back)
    end
end

# update batch
function update_batch(c::Client, batch::Vector{Int64})
    c.batch .= batch
end

# update gradient information
function update_grads(c::Client, grads::Matrix{Float64})
    c.grads .= grads
end

# update client model W
function update_model(c::Client)
    if typeof(c.W) <: Matrix{Float64}
        Xbatch = c.Xtrain[:, c.batch]
        Wgrad = (c.grads * Xbatch') ./ c.batch_size
        c.W .-= c.learning_rate * Wgrad
    else
        gs = c.back(c.grads)
        w = params(c.W)
        for i = 1:length(w)
            w[i] .-= (c.learning_rate / c.batch_size) * gs[w[i]]
        end
    end
end

########################################################################
# Client for asynchronous learning
########################################################################
mutable struct AsynClient{T1<:Int64, T2<:Float64, T3<:Matrix{T2}, T4<:SparseMatrixCSC{T2, T1}, T5<:Flux.Chain}
    id::T1                                  # client index
    Xtrain::T4                              # training data
    Xtest::T4                               # test data
    num_classes::T1                         # number of classes
    num_clients::T1                         # number of clients
    batch_size::T1                          # number of batches
    learning_rate::T2                       # learning rate
    W::Union{T3, T5}                        # client model
    grads::T3                               # gradient information
    ts::T2                                  # time gap between successive communications
    num_commu::T1                           # number of communication rounds
    back                                    # backward operator
    function AsynClient(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Xtest::SparseMatrixCSC{Float64, Int64}, ts::Float64, config::Dict{String, Union{Int64, Float64, String}})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        batch_size = config["batch_size"]
        learning_rate = config["learning_rate"]
        dm = size(Xtrain, 1)
        if config["local_model"] == "linear"
            W = zeros(Float64, num_classes, dm)
        elseif config["local_model"] == "mlp"
            W = Chain(Dense(dm, 32, σ), Dense(32, num_classes))
        end
        back = undef
        grads = zeros(Float64, num_classes, batch_size )
        num_commu = 0
        new{Int64, Float64, Matrix{Float64}, SparseMatrixCSC{Float64, Int64}, Flux.Chain}(id, Xtrain, Xtest, num_classes, num_clients, batch_size, learning_rate, W, grads, ts, num_commu, back)
    end
end

function update_grads(c::AsynClient, grads::Matrix{Float64})
    @printf "Client %i finish downloading gradient \n" c.id
    c.num_commu += 1
    c.grads = grads
end

function update_model(c::AsynClient, batch::Vector{Int64})
    if typeof(c.W) <: Matrix{Float64}
        Xbatch = c.Xtrain[:, batch]
        Wgrad = (c.grads * Xbatch') ./ c.batch_size
        c.W .-= c.learning_rate * Wgrad
    else
        gs = c.back(c.grads)
        w = params(c.W)
        for i = 1:length(w)
            w[i] .-= (c.learning_rate / c.batch_size) * gs[w[i]]
        end
    end
end


text:
"""
    Try{T} = Union{Const{<:Exception}, Identity{T}}
    @Try(error("something happend")) isa Const(<:Thrown{ErrorException})
    @Try(:successfull) == Identity(:successfull)

A specific case of [`Either`](@ref), where the Failure is always an Exception.
This can be used as an alternative to using try-catch syntax and allows for very flexible
error handling, as the error is now captured in a proper defined type.
Often it is really handy to treat errors like other values (without the need of extra try-catch syntax which only
applies to exceptions).

We reuse [`Identity`](@ref) for representing the single-element-container and `Const(<:Exception)` as the Exception
thrown.
"""
const Try{T} = Union{Const{<:Exception}, Identity{T}}
Try(t) = Identity(t)
Try(t::Exception) = Const(t)

"""
    Thrown(exception, stacktrace)

Thrown is like Exception, however can also cary stacktraces
"""
struct Thrown{E, S} <: Exception
  exception::E
  stacktrace::S
end

# == controversy https://github.com/JuliaLang/julia/issues/4648
function Base.:(==)(a::Thrown, b::Thrown)
  a.exception == b.exception && a.stacktrace == b.stacktrace
end

function Base.show(io::IO, x::Thrown)
  print(io, "Thrown($(repr(x.exception)))")
end

# Multiline version, following https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing-1
function Base.show(io::IO, ::MIME"text/plain", exc::Thrown)
  println(io, "Thrown($(repr(exc.exception)))")
  for (exc′, bt′) in exc.stacktrace
    showerror(io, exc′, bt′)
    println(io)
  end
end
function Base.showerror(io::IO, ::MIME"text/plain", exc::Thrown)
  Base.show(io, MIME"text/plain"(), exc)
end

"""
    MultipleExceptions(exception1, exception2, ...)
    MultipleExceptions(tuple_or_vector_of_exceptions)

Little helper type which combines several Exceptions into one new Exception.

In the several arguments version, and only there, if an MultipleExceptions is given, it will be flattened directly for convenience.
"""
struct MultipleExceptions{Es<:Tuple} <: Exception
  exceptions::Es
  function MultipleExceptions(exceptions::Tuple)
    @assert(isempty(exceptions) || eltype(exceptions) <: Exception,
      "expecting tuple of exception types, however got `eltype(exceptions) = $(eltype(exceptions))`")
    new{typeof(exceptions)}(exceptions)
  end
end
MultipleExceptions(exceptions::Vector) = MultipleExceptions(tuple(exceptions...))
function MultipleExceptions(args...)
  exceptions = ()
  for a in args
    exceptions = a isa MultipleExceptions ? tuple(exceptions..., a.exceptions...) : tuple(exceptions..., a)
  end
  MultipleExceptions(exceptions)
end

# == controversy https://github.com/JuliaLang/julia/issues/4648
Base.:(==)(a::MultipleExceptions, b::MultipleExceptions) = a.exceptions == b.exceptions

# Multiline version, following https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing-1
function Base.show(io::IO, ::MIME"text/plain", exc::MultipleExceptions)
  println(io, "MultipleExceptions:")
  for exception in exc.exceptions
    println(repr(exception))
  end
end
function Base.showerror(io::IO, ::MIME"text/plain", exc::MultipleExceptions)
  Base.show(io, MIME"text/plain"(), exc)
end

Base.merge(e1::Exception, e2::Exception) = MultipleExceptions((e1, e2))
Base.merge(es::MultipleExceptions, e::Exception) = MultipleExceptions(tuple(es.exceptions..., e))
Base.merge(e::Exception, es::MultipleExceptions) = MultipleExceptions(tuple(e, es.exceptions...))
Base.merge(es1::MultipleExceptions, es2::MultipleExceptions) = MultipleExceptions(tuple(es1.exceptions..., es2.exceptions...))

# we use a macro instead of dispatching on Try(f::Function) as this interferes e.g. with mapn
# (in mapn anonymous functions are passed through, which should not get executed automatically)
"""
    @Try begin
      your_code
    end

Macro which directly captures an Excpetion into a proper `Try`representation.

It translates to
```julia
try
  r = your_code
  Identity(r)
catch exc
  Const(Thrown(exc, Base.catch_stack()))
end
```
"""
macro Try(expr)
  quote
    try
      r = $(esc(expr))
      Identity(r)
    catch exc
      Const(Thrown(exc, Base.catch_stack()))
    end
  end
end

"""
    @TryCatch YourException begin
      your_code
    end

A version of [`@Try`](@ref) which catches only specific errors.
Every other orrer will be `rethrown`.

It translates to
```julia
try
  r = your_code
  Identity(r)
catch exc
  if exc isa YourException
    Const(Thrown(exc, Base.catch_stack()))
  else
    rethrow()
  end
end
```
"""
macro TryCatch(exception, expr)
  quote
    try
      r = $(esc(expr))
      Identity(r)
    catch exc
      if exc isa $(esc(exception))
        Const(Thrown(exc, Base.catch_stack()))
      else
        rethrow()
      end
    end
  end
end

"""
    isoption(::Const{Nothing}) = true
    isoption(::Identity) = true
    isoption(other) = false

check whether something is a [`Try`](@ref)
"""
istry(::Identity) = true
istry(::Const{<:Exception}) = true
istry(other) = false

"""
    issuccess(::Identity) = true
    issuccess(::Const{<:Exception}) = false

Similar to [`isright`](@ref), but only defined for `Const{<:Exception}`. Will
throw MethodError when applied on other `Const`.
"""
issuccess(::Identity) = true
issuccess(::Const{<:Exception}) = false

"""
    isfailure(::Identity) = false
    isfailure(::Const{<:Exception}) = true

Similar to [`isleft`](@ref), but only defined for `Const{<:Exception}`. Will
throw MethodError when applied on other `Const`.
"""
isfailure(::Identity) = false
isfailure(::Const{<:Exception}) = true

Base.eltype(::Type{Try{T}}) where T = T
Base.eltype(::Type{Try}) = Any


text:
module QSFit

export QSO, parent_recipe, add_spec!, logio, close_logio, PreparedSpectrum

import GFit: Domain, CompEval,
    Parameter, AbstractComponent, prepare!, evaluate!, fit!

using CMPFit, GFit
using Pkg, Pkg.Artifacts
using Statistics, DelimitedFiles, Dierckx, Printf, DataStructures
using Unitful, UnitfulAstro
using Dates
using SpecialFunctions
using Gnuplot
using TextParse, DataFrames

include("cosmology.jl")
include("ccm_unred.jl")
include("components/powerlaw.jl")
include("components/cutoff_powerlaw.jl")
include("components/hostgalaxy.jl")
include("components/ironopt.jl")
include("components/ironuv.jl")
include("components/balmercont.jl")
include("components/SpecLineGauss.jl")
include("components/SpecLineAsymmGauss.jl")
include("components/SpecLineLorentz.jl")
include("components/SpecLineVoigt.jl")
include("utils.jl")
include("convolutions.jl")
include("Spectrum.jl")

abstract type AbstractRecipe end

function default_options(::Type{T}) where T <: AbstractRecipe
    out = OrderedDict{Symbol, Any}()
    out[:wavelength_range] = [1215, 7.3e3]
    out[:min_spectral_coverage] = Dict{Symbol, Float64}(:default => 0.6)
    out[:skip_lines] = Vector{Symbol}()
    return out
end

struct QSO{T <: AbstractRecipe}
    name::String
    z::Float64
    mw_ebv::Float64
    cosmo::Cosmology.AbstractCosmology
    flux2lum::Float64
    logfile::Union{Nothing, String}
    options::OrderedDict{Symbol, Any}
    specs::Vector{Spectrum}
end

function QSO{T}(name, z; ebv=0., logfile=nothing, cosmo=default_cosmology()) where T <: AbstractRecipe
    @assert z > 0
    @assert ebv >= 0
    ld = uconvert(u"cm", luminosity_dist(cosmo, float(z)))
    flux2lum = 4pi * ld^2 * (scale_flux() * unit_flux()) / (scale_lum() * unit_lum())
    return QSO{T}(string(name), float(z), float(ebv), cosmo, flux2lum, logfile,
                  default_options(T), Vector{Spectrum}())
end

parent_recipe(source::QSO{T}) where T <: AbstractRecipe =
    QSO{supertype(T)}(getfield.(Ref(source), fieldnames(typeof(source)))...)

add_spec!(source::QSO, spec::Spectrum) =
    push!(source.specs, spec)

const logio_streams = Dict{String, IOStream}()

function logio(source::QSO)
    if isnothing(source.logfile)
        GFit.showsettings.plain = false
        return stdout
    end
    if !haskey(logio_streams, source.logfile)
        if isfile(source.logfile)
            f = open(source.logfile, "a")
        else
            f = open(source.logfile, "w")
        end
        logio_streams[source.logfile] = f
        println(f, "Timestamp: ", now())
        GFit.showsettings.plain = true
    end
    return logio_streams[source.logfile]
end

function close_logio(source::QSO)
    if haskey(logio_streams, source.logfile)
        close(logio_streams[source.logfile])
        delete!(logio_streams, source.logfile)
        GFit.showsettings.plain = false
    end
end

include("SpectralLines.jl")

struct PreparedSpectrum
    id::Int
    orig::Spectrum
    domain::GFit.Domain{1}
    data::GFit.Measures{1}
    lcs::OrderedDict{Symbol, LineComponent}
end

include("DefaultRecipe.jl")
#include("DefaultRecipe_multi.jl")
include("reduce.jl")
include("viewer.jl")
include("gnuplot.jl")
include("interactive_guess.jl")

end  # module


text:
using Test
using CompScienceMeshes

# G0 is the interior of the metal (no simulation required here)
# G1 is the unbounded exterior
# G2 is the bounded interior

fn = joinpath(@__DIR__, "assets", "thick_cladding.msh")
G01 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma01")
G02 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma02")
G12 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma12")

G = CompScienceMeshes.read_gmsh_mesh(fn)

@test numcells(G01) + numcells(G02) + numcells(G12) == numcells(G)

# Build region boundaries with normals pointing inward:
G1 = weld(-G01, -G12)
@test CompScienceMeshes.isoriented(G1)

G2 = weld(-G02, G12)
@test CompScienceMeshes.isoriented(G2)

isclosed(m) = (numcells(boundary(m)) == 0)

@test isclosed(G1)
@test isclosed(G2)


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct Attribute{S, T} <: AbstractAttribute
    iid::String
    type::AttributeType{S}
    value::T
end

function Attribute(t::Proto.Thing)
    iid = bytes2hex(t.iid)
    isempty(iid) && throw(TypeDBClientException(CONCEPT_MISSING_IID))

    attribute_type = instantiate(t._type)
    value_type = first(typeof(attribute_type).parameters)

    if value_type == Proto.AttributeType_ValueType.BOOLEAN
        return Attribute(iid, attribute_type, t.value.boolean)
    elseif value_type == Proto.AttributeType_ValueType.LONG
        return Attribute(iid, attribute_type, t.value.long)
    elseif value_type == Proto.AttributeType_ValueType.DOUBLE
        return Attribute(iid, attribute_type, t.value.double)
    elseif value_type == Proto.AttributeType_ValueType.STRING
        return Attribute(iid, attribute_type, t.value.string)
    elseif value_type == Proto.AttributeType_ValueType.DATETIME
        return Attribute(iid, attribute_type, t.value.date_time)
    else
        throw(TypeDBClientException(CONCEPT_BAD_VALUE_TYPE))
    end
end

"""
    get_owners(transaction::AbstractCoreTransaction,
        attribute::AbstractAttribute,
        thing_type::Optional{AbstractThingType} = nothing)

gives back all things has this attribute
"""
function get_owners(
    transaction::AbstractCoreTransaction,
    attribute::AbstractAttribute,
    thing_type::Optional{AbstractThingType} = nothing)

    iid = attribute.iid
    owner_req = thing_type === nothing ?
                                ThingRequestBuilder.attribute_get_owners_req(iid) :
                                ThingRequestBuilder.attribute_get_owners_req(iid, proto(thing_type))

    res_owners = stream(transaction, owner_req)
    return instantiate.(collect(Iterators.flatten(
        r.thing_res_part.attribute_get_owners_res_part.things for r in res_owners)))
end


text:
module SwapStreams

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end SwapStreams

using MappedArrays

export BigEndian, LittleEndian, SwapStream

const LittleEndian = 0x01020304
const BigEndian = 0x04030201

"""
    SwapStream([file_endianness],  io)

Given a file endian type and an IO stream returns a type that automatically byte
swaps to and from the appropriate endian type.

## Examples
```jldoctest
julia> using SwapStreams

julia> s = SwapStream(IOBuffer());  # assume byte swapping is necessary

julia> write(s, [1:10...]);         # byte swap each element before writing to buffer

julia> seek(s, 0);

julia> read!(s.io, Vector{Int}(undef, 10))  # raw data from buffer
10-element Array{Int64,1}:
  72057594037927936
 144115188075855872
 216172782113783808
 288230376151711744
 360287970189639680
 432345564227567616
 504403158265495552
 576460752303423488
 648518346341351424
 720575940379279360

julia> seek(s, 0);

julia> read!(s, Vector{Int}(undef, 10))  # byte swapped data from buffer
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

```
"""
struct SwapStream{S,IOType} <: IO
    io::IOType

    function SwapStream{S}(io::IOType) where {S,IOType<:IO}
        if S isa Bool
            return new{S,IOType}(io)
        else
            throw(TypeError(:SwapStream, Bool, S))
        end
    end

    SwapStream(io::IO) = SwapStream{true}(io)

    function SwapStream(file_endianness::UInt32, io::IOType) where {IOType}
        if file_endianness === ENDIAN_BOM
            return new{false,IOType}(io)
        else
            return new{true,IOType}(io)
        end
    end
end

is_swapping(x) = is_swapping(typeof(x))
is_swapping(::Type{<:SwapStream{S}}) where {S} = S

Base.seek(s::SwapStream, n::Integer) = seek(s.io, n)
Base.position(s::SwapStream)  = position(s.io)
Base.skip(s::SwapStream, n::Integer) = skip(s.io, n)
Base.eof(s::SwapStream) = eof(s.io)
Base.isreadonly(s::SwapStream) = isreadonly(s.io)
Base.isreadable(s::SwapStream) = isreadable(s.io)
Base.iswritable(s::SwapStream) = iswritable(s.io)
Base.stat(s::SwapStream) = stat(s.io)
Base.close(s::SwapStream) = close(s.io)
Base.isopen(s::SwapStream) = isopen(s.io)
Base.ismarked(s::SwapStream) = ismarked(s.io)
Base.mark(s::SwapStream) = mark(s.io)
Base.unmark(s::SwapStream) = unmark(s.io)
Base.reset(s::SwapStream) = reset(s.io)
Base.seekend(s::SwapStream) = seekend(s.io)

Base.read(s::SwapStream, n::Int) = read(s.io, n)

Base.read!(s::SwapStream{S}, a::Array{Int8}) where {S} = read!(s.io, a)
Base.read!(s::SwapStream{S}, a::Array{UInt8}) where {S} = read!(s.io, a)

# TODO should this use ArrayInterface.ismutable?
function Base.read!(s::SwapStream{S}, a::AbstractArray{T}) where {S,T}
    read!(s.io, a)
    if S
        @inbounds for i in eachindex(a)
            a[i] = bswap(a[i])
        end
    end
    return a
end

Base.read(s::SwapStream, ::Type{Int8}) = read(s.io, Int8)
Base.read(s::SwapStream, ::Type{UInt8}) = read(s.io, UInt8)

function Base.read(
    s::SwapStream{S},
    T::Union{Type{Float16}, Type{Float32}, Type{Float64}, Type{Int128}, Type{Int16}, Type{Int32}, Type{Int64}, Type{UInt128}, Type{UInt16}, Type{UInt32}, Type{UInt64}}
) where {S}

    if S
        return bswap(Base.read!(s.io, Ref{T}(0))[]::T)
    else
        return Base.read!(s.io, Ref{T}(0))[]::T
    end
end

function Base.read!(s::SwapStream{S}, ref::Base.RefValue{<:NTuple{N,T}}) where {S,N,T}
    if S
        return map(bswap, read!(s.io, ref)[])
    else
        return read!(s.io, ref)[]
    end
end

function Base.write(s::SwapStream{S}, x::Array) where {S}
    if S
        return write(s.io, mappedarray(ntoh, hton, x))
    else
        return write(s.io, x)
    end
end

function Base.write(s::SwapStream{S}, x::AbstractArray) where {S}
    if S
        return write(s.io, mappedarray(ntoh, hton, x))
    else
        return write(s.io, x)
    end
end

function Base.write(
    s::SwapStream{S},
    x::Union{Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128,Float16,Float32,Float64}
) where {S}

    if S
        return write(s.io, bswap(x))
    else
        return write(s.io, x)
    end
end

end # module


text:
using MultivariateSeries

X = @ring x0 x1 x2
n = length(X)
d = 4
r = 4

Xi = rand(n,r)
w = fill(1.0,r)
T = tensor(w,Xi,X, d)

k = 2
H = hankel(T,k)
P = perp(T,k)

@assert length(P)==2


text:
#  Solves a linear elliptic PDE in 1D
  using FEMfunctions

  using LinearAlgebra
  using SparseArrays

  nElements = 40

  include("onedElliptic.jl")

  order = 1
  error = onedElliptic(nElements,order)
  @test error<1e-10

  order = 2
  error = onedElliptic(nElements,order)
  @test error<1e-10

  order = 3
  error = onedElliptic(nElements,order)
  @test error<1e-10


text:
using GridInterpolations
using Statistics
using BenchmarkTools

function benchmark_interpolate(grid, data, rand_points)
    for point in eachcol(rand_points)
        interpolate(grid, data, point)
    end
end

function benchmark(GridType::Type, n_dims::Int; n_points_per_dim::Int=15, n_rand_points::Int=100, verbosity::Int=1)

    # construct grid
    grid = GridType((range(0, 1, length=n_points_per_dim) for _ in 1:n_dims)...)
    rand_points = rand(n_dims, n_rand_points)
    rand_data = rand(length(grid))

    trial = @benchmark benchmark_interpolate($grid, $rand_data, $rand_points) samples=10000

    if verbosity > 0
        println("Benchmarking interpolation on $(label(grid))")
        println("dimensionality:          $(n_dims)")
        println("points per dimension:    $(n_points_per_dim)")
        println("number of random points: $(n_rand_points)")
        if verbosity > 1
            show(stdout, MIME"text/plain"(), trial)  # make sure printing is verbose by using MIME
        else
            println("\n", trial)
        end
        println("\n--------------------------------------------------")
    end

    return trial
end

verbosity = 1
for ndims in 1:6

    println("\n##################################################")
    println("### Benchmark for ndims = $(ndims)")
    println("##################################################")

    trial_rectangle = benchmark(RectangleGrid, ndims; verbosity=verbosity)
    trial_simplex = benchmark(SimplexGrid, ndims; verbosity=verbosity)

    # println(ratio(trial_estimate_rectangle, trial_estimate_simplex))
    println(
        "RectangleGrid vs SimplexGrid: ",
        judge(median(trial_rectangle), median(trial_simplex))
    )
    println("##################################################\n")

end


text:
# Map a rule over the grids it reads from, updating the grids it writes to.
# This is broken into setup methods and application methods,
# for dispatch and to introduce a function barrier for type stability.

maprule!(data::AbstractSimData, rule) = maprule!(data, Val{ruletype(rule)}(), rule)
function maprule!(data::AbstractSimData, ruletype::Val{<:CellRule}, rule)
    rkeys, _ = _getreadgrids(rule, data)
    wkeys, _ = _getwritegrids(rule, data)
    maprule!(RuleData(data), ruletype, rule, rkeys, wkeys)
    return data
end
function maprule!(data::AbstractSimData, ruletype::Val{<:NeighborhoodRule}, rule)
    rkeys, rgrids = _getreadgrids(rule, data)
    wkeys, wgrids = _getwritegrids(rule, data)
    # NeighborhoodRule will read from surrounding cells
    # so may incorporate cells from the masked area
    _maybemask!(rgrids)
    # Copy or zero out boundary where needed
    _updateboundary!(rgrids)
    _cleardest!(data[neighborhoodkey(rule)])
    maprule!(RuleData(data), ruletype, rule, rkeys, wkeys)
    # Swap the dest/source of grids that were written to
    # and combine the written grids with the original simdata
    return _replacegrids(data, wkeys, _swapsource(_to_readonly(wgrids)))
end
function maprule!(data::AbstractSimData, ruletype::Val{<:SetRule}, rule)
    rkeys, _ = _getreadgrids(rule, data)
    wkeys, wgrids = _getwritegrids(rule, data)
    map(_astuple(wkeys, wgrids)) do g
        copyto!(parent(dest(g)), parent(source(g)))
    end
    ruledata = RuleData(_combinegrids(data, wkeys, wgrids))
    maprule!(ruledata, ruletype, rule, rkeys, wkeys)
    return _replacegrids(data, wkeys, _swapsource(_to_readonly(wgrids)))
end
function maprule!(data::AbstractSimData, ruletype::Val{<:SetGridRule}, rule)
    rkeys, rgrids = _getreadgrids(rule, data)
    wkeys, wgrids = _getwritegrids(rule, data)
    _maybemask!(rgrids) # TODO... mask wgrids?
    ruledata = RuleData(_combinegrids(data, wkeys, wgrids))
    # Run the rule
    applyrule!(ruledata, rule)
    return data
end
function maprule!(data::AbstractSimData, ruletype::Val, rule, rkeys, wkeys)
    maprule!(data, proc(data), opt(data), ruletype, rule, rkeys, wkeys)
end

function maprule!(data::AbstractSimData, proc::CPU, opt, ruletype::Val, rule, rkeys, wkeys)
    let data=data, proc=proc, opt=opt, rule=rule,
        rkeys=rkeys, wkeys=wkeys, ruletype=ruletype
        optmap(data, proc, opt, ruletype, rkeys) do i, j
            cell_kernel!(data, ruletype, rule, rkeys, wkeys, i, j)
        end
    end
end
function maprule!(
    data::AbstractSimData, proc::CPU, opt, ruletype::Val{<:NeighborhoodRule},
    rule, args...
)
    hoodgrid = data[neighborhoodkey(rule)]
    _mapneighborhoodrule!(data, hoodgrid, proc, opt, ruletype, rule, args...)
end

function _mapneighborhoodrule!(
    data, hoodgrid::GridData{<:Tuple{Y,X},R}, proc, opt, ruletype, rule, rkeys, wkeys
) where {Y,X,R}
    let data=data, hoodgrid=hoodgrid, proc=proc, opt=opt, ruletyp=ruletype, rule=rule, rkeys=rkeys, wkeys=wkeys
        B = 2R
        # UNSAFE: we must avoid sharing status blocks, it could cause race conditions
        # when setting status from different threads. So we split the grid in 2 interleaved
        # sets of rows, so that we never run adjacent rows simultaneously
        procmap(proc, 1:2:_indtoblock(Y, B)) do bi
            row_kernel!(data, hoodgrid, proc, opt, ruletype, rule, rkeys, wkeys, bi)
        end
        procmap(proc, 2:2:_indtoblock(Y, B)) do bi
            row_kernel!(data, hoodgrid, proc, opt, ruletype, rule, rkeys, wkeys, bi)
        end
    end
    return nothing
end

### Rules that don't need a neighborhood buffer ####################

# optmap
# Map kernel over the grid, specialising on PerformanceOpt.

# Run kernels with SparseOpt, block by block:
function optmap(
    f, simdata::AbstractSimData{S}, proc, ::SparseOpt, ruletype::Val{<:Rule}, rkeys
) where {S<:Tuple{Y,X}} where {Y,X}
    # Only use SparseOpt for single-grid rules with grid radii > 0
    grid = _firstgrid(simdata, rkeys)
    R = radius(grid)
    if R == 0
        optmap(f, simdata, proc, NoOpt(), ruletype, rkeys)
        return nothing
    end
    B = 2R
    status = sourcestatus(grid)
    let f=f, proc=proc, status=status
        procmap(proc, 1:_indtoblock(X+R, B)) do bj
            for  bi in 1:_indtoblock(Y+R, B)
                status[bi, bj] || continue
                # Convert from padded block to init dimensions
                istart, jstart = _blocktoind(bi, B) - R, _blocktoind(bj, B) - R
                # Stop at the init row/column size, not the padding or block multiple
                istop, jstop = min(istart + B - 1, Y), min(jstart + B - 1, X)
                # Skip the padding
                istart, jstart  = max(istart, 1), max(jstart, 1)
                for j in jstart:jstop
                    @simd for i in istart:istop
                        f(i, j)
                    end
                end
            end
        end
    end
    return nothing
end
# Run kernel over the whole grid, cell by cell:
function optmap(
    f, simdata::AbstractSimData{S}, proc, ::NoOpt, ::Val{<:Rule}, rkeys
) where S<:Tuple{Y,X} where {Y,X}
    procmap(proc, 1:X) do j
        for i in 1:Y
            f(i, j) # Run rule for each row in column j
        end
    end
end
function optmap(
    f, simdata::AbstractSimData{S}, proc, ::NoOpt, ::Val{<:CellRule}, rkeys
) where S<:Tuple{Y,X} where {Y,X}
    procmap(proc, 1:X) do j
        @simd for i in 1:Y
            f(i, j) # Run rule for each row in column j
        end
    end
end
# procmap
# Map kernel over the grid, specialising on Processor
# Looping over cells or blocks on CPU
@inline procmap(f, proc::SingleCPU, range) =
    for n in range
        f(n) # Run rule over each column
    end
@inline procmap(f, proc::ThreadedCPU, range) =
    Threads.@threads for n in range
        f(n) # Run rule over each column, threaded
    end

# cell_kernel!
# runs a rule for the current cell
@inline function cell_kernel!(simdata, ruletype::Val{<:Rule}, rule, rkeys, wkeys, i, j)
    readval = _readcell(simdata, rkeys, i, j)
    writeval = applyrule(simdata, rule, readval, (i, j))
    _writecell!(simdata, ruletype, wkeys, writeval, i, j)
    writeval
end
@inline function cell_kernel!(simdata, ::Val{<:SetRule}, rule, rkeys, wkeys, i, j)
    readval = _readcell(simdata, rkeys, i, j)
    applyrule!(simdata, rule, readval, (i, j))
    nothing
end

# row_kernel!
# Run a NeighborhoodRule rule row by row. When we move along a row by one cell, we
# access only a single new column of data with the height of 4R, and move the existing
# data in the neighborhood buffers array across by one column. This saves on reads
# from the main array.
function row_kernel!(
    simdata::AbstractSimData, grid::GridData{<:Tuple{Y,X},R}, proc, opt::NoOpt,
    ruletype::Val, rule::Rule, rkeys, wkeys, bi
) where {Y,X,R}
    B = 2R
    i = _blocktoind(bi, B)
    i > Y && return nothing
    # Loop along the block ROW.
    src = parent(source(grid))
    buffers = _initialise_buffers(src, Val{R}(), i, 1)
    blocklen = min(Y, i + B - 1) - i + 1
    for j = 1:X
        buffers = _update_buffers(buffers, src, Val{R}(), i, j)
        # Loop over the COLUMN of buffers covering the block
        for b in 1:blocklen
            @inbounds bufrule = _setbuffer(rule, buffers[b])
            cell_kernel!(simdata, ruletype, bufrule, rkeys, wkeys, i + b - 1, j)
        end
    end
    return nothing
end
function row_kernel!(
    simdata::AbstractSimData, grid::GridData{<:Tuple{Y,X},R}, proc, opt::SparseOpt,
    ruletype::Val, rule::Rule, rkeys, wkeys, bi
) where {Y,X,R}
    B = 2R
    S = 2R + 1
    nblockcols = _indtoblock(X+R, B)
    src = parent(source(grid))
    srcstatus, dststatus = sourcestatus(grid), deststatus(grid)

    # Blocks ignore padding! the first block contains padding.
    i = _blocktoind(bi, B)
    i > Y && return nothing
    # Get current bloc
    skippedlastblock = true

    # Initialise block status for the start of the row
    # The first column always runs, it's buggy otherwise.
    @inbounds bs11, bs12 = true, true
    @inbounds bs21, bs22 = true, true
    # New block status
    newbs12 = false
    newbs22 = false
    buffers = _initialise_buffers(src, Val{R}(), i, 1)
    for bj = 1:nblockcols
        # Shuffle current buffer status
        bs11, bs21 = bs12, bs22
        @inbounds bs12, bs22 = srcstatus[bi, bj + 1], srcstatus[bi + 1, bj + 1]
        # Skip this block if it and the neighboring blocks are inactive
        if !(bs11 | bs12 | bs21 | bs22)
            skippedlastblock = true
            # Run the rest of the chain if it exists and more than 1 grid is used
            if rule isa Chain && length(rule) > 1 && length(rkeys) > 1
                # Loop over the grid COLUMNS inside the block
                jstart = _blocktoind(bj, B)
                jstop = min(jstart + B - 1, X)
                for j in jstart:jstop
                    # Loop over the grid ROWS inside the block
                    blocklen = min(Y, i + B - 1) - i + 1
                    for b in 1:blocklen
                        cell_kernel!(simdata, ruletype, rule, rkeys, wkeys, i + b - 1, j)
                    end
                end
            end
            continue
        end
        # Define area to loop over with the block.
        # It's variable because the last block may be partial
        jstart = _blocktoind(bj, B)
        jstop = min(jstart + B - 1, X)

        # Reinitialise neighborhood buffers if we have skipped a section of the array
        if skippedlastblock
            buffers = _initialise_buffers(src, Val{R}(), i, jstart)
            skippedlastblock = false
        end
        # Shuffle new buffer status
        newbs11 = newbs12
        newbs21 = newbs22
        newbs12 = newbs22 = false

        # Loop over the grid COLUMNS inside the block
        for j in jstart:jstop
            # Update buffers unless feshly populated
            buffers = _update_buffers(buffers, src, Val{R}(), i, j)
            # Which block column are we in, 1 or 2
            curblockj = (j - jstart) ÷ R + 1
            # Loop over the COLUMN of buffers covering the block
            blocklen = min(Y, i + B - 1) - i + 1
            for b in 1:blocklen
                # Set rule buffer
                bufrule = _setbuffer(rule, buffers[b])
                # Run the rule kernel for the cell
                writeval = cell_kernel!(simdata, ruletype, bufrule, rkeys, wkeys, i + b - 1, j)
                # Update the status for the current block
                cs = _cellstatus(opt, wkeys, writeval)
                curblocki = R == 1 ? b : (b - 1) ÷ R + 1
                if curblocki == 1
                    curblockj == 1 ? (newbs11 |= cs) : (newbs12 |= cs)
                else
                    curblockj == 1 ? (newbs21 |= cs) : (newbs22 |= cs)
                end
            end

            # Combine new block status with deststatus array
            @inbounds dststatus[bi, bj] |= newbs11
            @inbounds dststatus[bi+1, bj] |= newbs21
            @inbounds dststatus[bi, bj+1] |= newbs12
            @inbounds dststatus[bi+1, bj+1] |= newbs22
        end
    end
    return nothing
end

#### Utils

@inline _cellstatus(opt::SparseOpt, wkeys::Tuple, writeval) = _isactive(writeval[1], opt)
@inline _cellstatus(opt::SparseOpt, wkeys, writeval) = _isactive(writeval, opt)

@inline _firstgrid(simdata, ::Val{K}) where K = simdata[K]
@inline _firstgrid(simdata, ::Tuple{Val{K},Vararg}) where K = simdata[K]

_to_readonly(data::Tuple) = map(ReadableGridData, data)
_to_readonly(data::WritableGridData) = ReadableGridData(data)

# _maybemask!
# mask the source grid
_maybemask!(wgrids::Union{Tuple,NamedTuple}) = map(_maybemask!, wgrids)
_maybemask!(wgrid::GridData) = _maybemask!(wgrid, proc(wgrid), mask(wgrid))
_maybemask!(wgrid::GridData, proc, mask::Nothing) = nothing
function _maybemask!(
    wgrid::GridData{<:Tuple{Y,X}}, proc::CPU, mask::AbstractArray
) where {Y,X}
    procmap(proc, 1:X) do j
        @simd for i in 1:Y
            source(wgrid)[i, j] *= mask[i, j]
        end
    end
end
function _maybemask!(
    wgrid::GridData{<:Tuple{Y,X}}, proc, mask::AbstractArray
) where {Y,X}
    sourceview(wgrid) .*= mask
end

# _cleardest!
# Clear the destination grid and its status for SparseOpt.
# NoOpt writes to every cell, so this is not required
_cleardest!(grid) = _cleardest!(grid, opt(grid))
function _cleardest!(grid, opt::SparseOpt)
    dest(grid) .= source(grid)
    deststatus(grid) .= false
end
_cleardest!(grid, opt) = nothing


text:
@testset "map" begin
    ds = KeyedDataset(
        flatten([
            :g1 => [
                :a => KeyedArray(zeros(3); time=1:3),
                :b => KeyedArray(ones(3, 2); time=1:3, loc=[:x, :y]),
            ],
            :g2 => [
                :a => KeyedArray(ones(3); time=1:3),
                :b => KeyedArray(zeros(3, 2); time=1:3, loc=[:x, :y]),
            ]
        ])...
    );

    # Since all of the tests are reusing `ds` we can assum it isn't being mutated.
    # Maybe there's a better way to test the combinations, but this seemed pretty terse?

    # Test mapping over every component
    r = map(a -> a .+ 1.0, ds)
    @test mean(r[(:g1, :a)]) == mean(r[(:g2, :b)]) == 1.0
    @test mean(r[(:g1, :b)]) == mean(r[(:g2, :a)]) == 2.0

    # Test mapping over components
    r = map(a -> a .+ 1.0, ds, (:g1, :__))
    @test mean(r[(:g2, :b)]) == 0.0
    @test mean(r[(:g1, :a)]) == mean(r[(:g2, :a)]) == 1.0
    @test mean(r[(:g1, :b)]) == 2.0

    # Test mapping over more complex patterns
    r = map(a -> a .+ 1.0, ds, (:__, :a, :_))
    @test mean(r[(:g2, :b)]) == 0.0
    @test mean(r[(:g1, :a)]) == mean(r[(:g1, :b)]) == 1.0
    @test mean(r[(:g2, :a)]) == 2.0

    # Test mapping just dims
    r = map(a -> a .+ 1.0, ds, (:__, :loc))
    @test mean(r[(:g1, :a)]) == 0.0
    @test mean(r[(:g2, :b)]) == mean(r[(:g2, :a)]) == 1.0
    @test mean(r[(:g1, :b)]) == 2.0

    # Test mapping over component and dim pattern
    r = map(a -> a .+ 1.0, ds, (:g1, :__, :loc))
    @test mean(r[(:g1, :a)]) == mean(r[(:g2, :b)]) == 0.0
    @test mean(r[(:g2, :a)]) == 1.0
    @test mean(r[(:g1, :b)]) == 2.0
end

@testset "mapslices" begin
    ds = KeyedDataset(
        flatten([
            :g1 => [
                :a => KeyedArray(zeros(3); time=1:3),
                :b => KeyedArray(ones(3, 2); time=1:3, loc=[:x, :y]),
            ],
            :g2 => [
                :a => KeyedArray(ones(3); time=1:3),
                :b => KeyedArray(zeros(3, 2); time=1:3, loc=[:x, :y]),
            ]
        ])...
    )

    r = mapslices(sum, ds; dims=:time)
    expected = [
        (:g1, :a) => [0.0],
        (:g1, :b) => [3.0 3.0],
        (:g2, :a) => [3.0],
        (:g2, :b) => [0.0 0.0],
    ]
    for (k, v) in expected
        @test v == r[k]
    end

    r = mapslices(sum, ds, (:__, :b, :_); dims=:loc)
    @test r == mapslices(sum, ds; dims=:loc)
    expected = [
        (:g1, :a) => zeros(3),
        (:g1, :b) => fill(2.0, (3, 1)),
        (:g2, :a) => ones(3),
        (:g2, :b) => zeros((3, 1)),
    ]
    for (k, v) in expected
        @test v == r[k]
    end

    # Reducing over time will violate our :time dimension constraint
    @test_throws KeyAlignmentError mapslices(sum, ds, (:__, :b, :_); dims=:time)
    # Dimension doesn't exist in the selection
    @test_throws ArgumentError mapslices(sum, ds, (:__, :b, :_); dims=:foo)
end

@testset "merge" begin
    # Test destructive merging of valid and invalid datasets
    @testset "additive" begin
        ds1 = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            )
        )

        ds2 = KeyedDataset(
            :val2 => KeyedArray(
                rand(4, 3, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            )
        )
        ds = merge(ds1, ds2)
        @test issetequal(keys(ds.data), [(:val1,), (:val2,)])

        # Make sure merge works on >2 datasets
        ds3 = KeyedDataset(
            :val3 => KeyedArray(
                rand(4, 3, 2) .+ 2.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            )
        )

        ds = merge(ds1, ds2, ds3)
        @test issetequal(keys(ds.data), [(:val1,), (:val2,), (:val3,)])
    end
    @testset "replace" begin
        ds1 = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            )
        )

        ds2 = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            )
        )
        ds = merge(ds1, ds2)
        @test issetequal(keys(ds.data), [(:val1,)])
        @test all(x -> x >= 1.0, ds.val1)
    end
    @testset "key mismatch" begin
        ds1 = KeyedDataset(
            :val1 => KeyedArray(
                rand(4, 3, 2);
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:3,
                obj=[:a, :b],
            )
        )

        ds2 = KeyedDataset(
            :val2 => KeyedArray(
                rand(4, 2, 2) .+ 1.0;
                time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
                loc=1:2,
                obj=[:a, :b],
            )
        )
        @test_throws KeyAlignmentError merge(ds1, ds2)
    end
end

@testset "rekey" begin
    ds = KeyedDataset(
        :val1 => KeyedArray(
            rand(4, 3, 2);
            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
            loc=1:3,
            obj=[:a, :b],
        ),
        :val2 => KeyedArray(
            rand(4, 3, 2) .+ 1.0;
            time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
            loc=1:3,
            obj=[:a, :b],
        );
    )

    r = AxisSets.rekey(k -> ZonedDateTime.(k, tz"UTC"), ds, :time)
    @test eltype(r.time) <: ZonedDateTime
    @test eltype(r.val1.time) <: ZonedDateTime
    @test eltype(r.val2.time) <: ZonedDateTime

    r = AxisSets.rekey(k -> k .+ 1, ds, (:__, :loc))
    @test r.loc == 2:4
    @test r.loc == 2:4
    @test r.loc == 2:4
end


text:
"""
Functionality for handling the defautl conent of the various Genie files (migrations, models, controllers, etc).
"""
module FileTemplates

using Genie.Inflector

"""
    new_task(module_name::String) :: String

Default content for a new Genie task.
"""
function newtask(module_name::String) :: String
  """
  module $module_name

  using Genie, Genie.Toolbox

  \"\"\"
  Description of the task here
  \"\"\"
  function run_task()
    # Build something great
  end

  end
  """
end
const new_task = newtask

"""
    newcontroller(controller_name::String) :: String

Default content for a new Genie controller.
"""
function newcontroller(controller_name::String) :: String
  """
  module $(controller_name)Controller
  # Build something great
  end
  """
end
const new_controller = newcontroller

"""
    newchannel(channel_name::String) :: String

Default content for a new Genie channel.
"""
function newchannel(channel_name::String) :: String
  """
  module $(channel_name)Channel

  using Genie.WebChannels, Genie.Router

  function subscribe()
    WebChannels.subscribe(wsclient(@params), :$(lowercase(channel_name)))
    "OK"
  end

  end
  """
end
const new_channel = newchannel

"""
    newtest(plural_name::String, singular_name::String) :: String

Default content for a new test file.
"""
function newtest(plural_name::String, singular_name::String) :: String
  """
  using Genie, App.$(plural_name)

  ### Your tests here
  @test 1 == 1
  """
end
const new_test = newtest

"""
"""
function appmodule(path::String)
  path = replace(path, "-"=>"_") |> strip
  appname = split(path, "/", keepempty = false)[end] |> String |> Inflector.from_underscores

  content = """
  module $appname

  using Genie, Genie.Router, Genie.Renderer, Genie.AppServer

  function main()
    Base.eval(Main, :(const UserApp = $appname))

    include("genie.jl")

    Base.eval(Main, :(const Genie = $appname.Genie))
    Base.eval(Main, :(using Genie))
  end

  main()

  end
  """

  (appname, content)
end

end


text:
using FANN

# took some inspiration from Julia NLopt package
typealias _DataSet Ptr{FANN.fann_train_data}

# Create DataSet from input data.
#
# Observations are organised in columns in the matrices X and y, to
# reflect the data format required by FANN, which uses an array of
# pointers to arrays to store the observations. Hence, we make use
# of Julia column-major order to save a copy of the data.
#
# Parameters
# ----------
# X : n_inp x n_obs matrix
# Y : n_out x n_obs vector
type DataSet
    data::_DataSet
    function DataSet(d::_DataSet)
        d = new(d)
        finalizer(d, destroy)
        d
    end
    function DataSet(X::Matrix, Y::Matrix)
   		if size(X, 2) != size(Y, 2)
   			error("sizes of X and Y do not match")
   		end
		num_input, num_data = size(X)
		num_output = size(Y, 1)
		d = ccall((:fann_create_train_array, libfann), _DataSet,
				  (Uint32, Uint32, Ptr{fann_type}, Uint32, Ptr{fann_type}),
				  num_data,
				  num_input,
				  convert(Ptr{fann_type}, pointer(X)),
				  num_output,
				  convert(Ptr{fann_type}, pointer(Y)))
        if d == C_NULL
            error("Error in fann_create_train")
        end
        DataSet(d)
    end
end

Base.convert(::Type{_DataSet}, d::DataSet) = d.data
Base.show(io::IO, d::DataSet) = print(io, "DataSet()")
destroy(d::DataSet) = ccall((:fann_destroy_train, libfann), Void, (_DataSet,), d)

# Save a dataset to file. Used for debugging.
function savedset(dset::DataSet, filename::ASCIIString)
	retval = ccall((:fann_save_train, libfann),
			       Cint,
			       (Ptr{fann_train_data}, Ptr{Uint8}),
			       dset, filename)
	retval != 0 && error("Error saving data")
    return retval
end

function loaddset(filename::ASCIIString)
    dset = ccall((:fann_read_train_from_file, libfann), Ptr{fann_train_data}, (Ptr{Uint8},), filename)
    dset == C_NULL && error("Error reading dataset")
    return DataSet(dset)
end