text:
using DataFrames
using LinearAlgebra

"""
    UniqueTypes

A subtype of AbstractTypes where all individuals are completely
distinct. This type is the simplest AbstractTypes subtype, which
identifies all individuals as unique and completely distinct from each
other.

"""
struct UniqueTypes <: Diversity.API.AbstractTypes
    num::Int64
    names::Vector{String}

    function UniqueTypes(num::Integer)
        num > 0 || error("Too few species")
        new(num, map(x -> "$x", 1:num))
    end

    function UniqueTypes(names::Vector{String})
        num = length(names)
        num > 0 || error("Too few species")
        new(num, names)
    end
end

import Diversity.API._hassimilarity
_hassimilarity(::UniqueTypes) = false

import Diversity.API._counttypes
function _counttypes(ut::UniqueTypes, ::Bool)
    return ut.num
end

import Diversity.API._gettypenames
function _gettypenames(ut::UniqueTypes, ::Bool)
    return ut.names
end

import Diversity.API._calcsimilarity
function _calcsimilarity(ut::UniqueTypes, ::Real)
    return I
end

import Diversity.API._calcordinariness
function _calcordinariness(::UniqueTypes, abundances::AbstractArray, ::Real)
    return abundances
end

import Diversity.API._getdiversityname
_getdiversityname(::UniqueTypes) = "Unique"

"""
    Species

A subtype of AbstractTypes where all species are completely distinct.
This type is the simplest AbstractTypes subtype, which identifies all
species as unique and completely distinct from each other.

"""
const Species = UniqueTypes

"""
    Taxonomy

A subtype of AbstractTypes with similarity between related taxa,
creating taxonomic similarity matrices.

"""
struct Taxonomy{FP <: AbstractFloat} <: Diversity.API.AbstractTypes
    speciesinfo::DataFrame
    taxa::Dict{Symbol, FP}
    typelabel::Symbol

    function Taxonomy{FP}(speciesinfo::DataFrame,
                          taxa::Dict{Symbol, FP},
                          typelabel::Symbol) where FP <: AbstractFloat
        sort(describe(speciesinfo)[!,:variable]) == sort([keys(taxa)...]) ||
            error("Taxon labels do not match similarity values")
        typelabel ∈ describe(speciesinfo)[!,:variable] ||
            error("$typelabel not found in DataFrame column names")
        new{FP}(speciesinfo, taxa, typelabel)
    end
end

function Taxonomy(speciesinfo::DataFrame, taxa::Dict,
                  typelabel::Symbol = :Species)
    Taxonomy{valtype(taxa)}(speciesinfo, taxa, typelabel)
end

import Diversity.API.floattypes
function floattypes(::Taxonomy{FP}) where FP <: AbstractFloat
    return Set([FP])
end

function _counttypes(tax::Taxonomy, ::Bool)
    return nrow(tax.speciesinfo)
end

function _gettypenames(tax::Taxonomy, ::Bool)
    return tax.speciesinfo[!,tax.typelabel]
end

function _calcsimilarity(::Taxonomy, ::Real)
    error("Can't generate a taxonomic similarity matrix yet")
end

_getdiversityname(::Taxonomy) = "Taxonomy"

"""
    GeneralTypes{FP, M}

An AbstractTypes subtype with a general similarity matrix. This
subtype simply holds a matrix with similarities between individuals.

# Members:

- `z` A two-dimensional matrix representing similarity between
individuals.
"""
struct GeneralTypes{FP <: AbstractFloat,
                    M <: AbstractMatrix{FP}} <: Diversity.API.AbstractTypes
    """
        z

    A two-dimensional matrix representing similarity between
    individuals.
    """
    z::M

    """
        names

    Optional vector of type names.
    """
    names::Vector{String}

    """
    # Constructor for GeneralTypes

    Creates an instance of the GeneralTypes class, with an arbitrary
    similarity matrix.
    """
    function GeneralTypes(zmatrix::M) where {FP <: AbstractFloat,
                                             M <: AbstractMatrix{FP}}
        size(zmatrix, 1) == size(zmatrix, 2) ||
            throw(DimensionMismatch("Similarity matrix is not square"))

        minimum(zmatrix) ≥ 0 || throw(DomainError(minimum(zmatrix),
                                      "Similarities must be ≥ 0"))
        maximum(zmatrix) ≤ 1 || @warn "Similarity matrix has values above 1"

        new{FP, M}(zmatrix, map(x -> "$x", 1:size(zmatrix, 1)))
    end

    function GeneralTypes(zmatrix::M, names::Vector{String}) where
        {FP <: AbstractFloat, M <: AbstractMatrix{FP}}
        size(zmatrix, 1) == size(zmatrix, 2) ||
            throw(DimensionMismatch("Similarity matrix is not square"))

        minimum(zmatrix) ≥ 0 || throw(DomainError(minimum(zmatrix),
                                      "Similarities must be ≥ 0"))
        maximum(zmatrix) ≤ 1 || @warn "Similarity matrix has values above 1"

        length(names) == size(zmatrix, 1) ||
            error("Species name vector does not match similarity matrix")

        new{FP, M}(zmatrix, names)
    end
end

function floattypes(::GeneralTypes{FP, M}) where {FP <: AbstractFloat,
                                                  M <: AbstractMatrix{FP}}
    return Set([FP])
end

function _counttypes(gt::GeneralTypes, ::Bool)
    return size(gt.z, 1)
end

function _gettypenames(gt::GeneralTypes, ::Bool)
    return gt.names
end

function _calcsimilarity(gt::GeneralTypes, ::Real)
    return gt.z
end

_getdiversityname(::GeneralTypes) = "Arbitrary Z"


text:
__precompile__(true)
"""
Case folding tables for Unicode characters

Copyright 2018-2020 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
"""
module CaseTables

struct CaseTable
    l_tab::NTuple{128,UInt8}
    u_tab::NTuple{128,UInt8}
    t_tab::NTuple{64,UInt8}

    is_l_tab::NTuple{256, UInt8}
    is_u_tab::NTuple{256, UInt8}
    can_l_tab::NTuple{256, UInt8}
    can_u_tab::NTuple{256, UInt8}

    can_l_flg::UInt128
    can_u_flg::UInt128
    can_t_flg::UInt128
    can_sl_flg::UInt128
    can_su_flg::UInt128

    is_l_flg::UInt128
    is_u_flg::UInt128
    is_sl_flg::UInt128
    is_su_flg::UInt128

    siz_l_flg::UInt128
    siz_u_flg::UInt128
    max_siz_l::UInt32
    max_siz_u::UInt32
end

# Calculate tables (later move these to library built based on Unicode tables with BinaryBuilder,
# loaded by BinaryProvider)

# For speed, we want to know which characters can be lowercased, uppercased, or titlecased
# We want to also know which ones increase or decrease in size in UTF-8 encoding

@inline utf8_len(cp::Unsigned) = ifelse(cp <= 0x7f, 1, ifelse(cp <= 0x7ff, 2, 3 + (cp > 0xffff)))
@inline utf8_len(ch::Char) = utf8_len(ch%UInt32)

const _z4 = 0%UInt32

const zerotup = (_z4, _z4, _z4, _z4, _z4, _z4, _z4, _z4,
                 _z4, _z4, _z4, _z4, _z4, _z4, _z4, _z4)

function _add_to_bv!(vec, tmp)
    tt = tuple(tmp...)
    tt == zerotup && return 0x0
    fill!(tmp, 0%UInt32)
    for (i, t) in enumerate(vec)
        (t == tt) && return UInt8(i)
    end
    push!(vec, tt)
    UInt8(length(vec))
end

# case tables should be:
# 128-bit bitmap to say whether a particular range of 512 code points can be folded
# 64-byte ntuple{32,UInt16}
function case_tables()
    bitvec = NTuple{16,UInt32}[]  # each element represents 512 characters
    tupvec = NTuple{32,UInt16}[]  # each element represents 32 characters
    offvec = NTuple{32,UInt8}[]   # each element maps 1024 character range to offsets in tupvec

    sizvecl = Pair{UInt16,UInt64}[]
    sizvecu = Pair{UInt16,UInt64}[]

    # Top level tables
    l_tab  = fill(0x0, 128)
    u_tab  = fill(0x0, 128)
    t_tab  = fill(0x0, 64)

    # The elements of these are all offsets into bitvec
    is_l_tab  = fill(0x0, 256)
    is_u_tab  = fill(0x0, 256)
    can_l_tab = fill(0x0, 256)
    can_u_tab = fill(0x0, 256)
    is_l_flg = is_u_flg = is_sl_flg = is_su_flg = 0%UInt128

    # These tables get reused, stored in offvec
    l_off = fill(0x0, 32)
    u_off = fill(0x0, 32)
    t_off = fill(0x0, 32)

    # These tables get reused, stored in tupvec
    tmp_l = fill(0x0000, 32)
    tmp_u = fill(0x0000, 32)
    tmp_t = fill(0x0000, 32)

    # These tables get reused, stored in bitvec
    bit_is_l  = fill(0%UInt32, 16)
    bit_is_u  = fill(0%UInt32, 16)
    bit_can_l = fill(0%UInt32, 16)
    bit_can_u = fill(0%UInt32, 16)

    # Handle BMP
    l_mid = u_mid = t_mid = false
    can_l_flg = can_u_flg = can_t_flg = 0%UInt128
    siz_l_flg = siz_u_flg = 0%UInt128
    max_siz_l = max_siz_u = 0%UInt32
    tmp_is_l = tmp_is_u = tmp_can_l = tmp_can_u = 0%UInt32
    for rng in (0x0080:0x20:0xd7e0, 0xe000:0x20:0xffe0), base in rng
        hipos  = (base >>> 10) + 1
        hibit  = UInt128(1) << (base >>> 9)
        midbits = ((base >>> 5) & 0x1f) + 1

        # Handle block of 32 characters
        l_flg = u_flg = t_flg = false
        diff_l = diff_u = 0%UInt64
        for off = 0x00:0x1f
            cp = UInt16(base + off)
            lowbit = UInt32(1) << off
            ch = Char(cp)
            cl = UInt16(lowercase(ch))
            cu = UInt16(uppercase(ch))
            ct = UInt16(titlecase(ch))
            tmp_l[off+1] = cl
            tmp_u[off+1] = cu
            tmp_t[off+1] = ct
            islowercase(ch) && (tmp_is_l |= lowbit)
            isuppercase(ch) && (tmp_is_u |= lowbit)
            cp === cl === cu === ct && continue
            sizc = utf8_len(ch)
            sizu = utf8_len(cu)
            sizt = utf8_len(ct)
            if cl !== cp
                l_flg = true
                tmp_can_l |= lowbit
                diff = utf8_len(cl) - sizc
                -2 <= diff <= 1 || error("Size difference for $cp -> $cl is $diff")
                diff == 0 || (diff_l |= (UInt64(diff & 3)<<(off<<1)); max_siz_l = cp)
            end
            if cu !== cp
                u_flg = true
                tmp_can_u |= lowbit
                diff = sizu - sizc
                -2 <= diff <= 1 || error("Size difference for $cp -> $cu is $diff")
                diff == 0 || (diff_u |= (UInt64(diff & 3)<<(off<<1)); max_siz_u = cp)
            end
            if ct !== cu
                t_flg = true
                sizu == sizt ||
                    error("Titlecase and Uppercase are not same size in UTF-8: " *
                          "$cp $cu:$sizu $ct:$sizt")
            end
        end

        bitoff = ((base >>> 5) & 0xf) + 1
        bit_is_l[bitoff]  = tmp_is_l
        bit_is_u[bitoff]  = tmp_is_u
        bit_can_l[bitoff] = tmp_can_l
        bit_can_u[bitoff] = tmp_can_u

        tmp_is_l = tmp_is_u = tmp_can_l = tmp_can_u = 0%UInt32

        if l_flg
            can_l_flg |= hibit
            l_mid = true # Have at least one in this set of blocks
            push!(tupvec, tuple(tmp_l...))
            l_off[midbits] = UInt8(length(tupvec))
        end
        if u_flg
            can_u_flg |= hibit
            u_mid = true
            push!(tupvec, tuple(tmp_u...))
            u_off[midbits] = UInt8(length(tupvec))
        end
        if t_flg
            can_t_flg |= hibit
            t_mid = true
            push!(tupvec, tuple(tmp_t...))
            t_off[midbits] = UInt8(length(tupvec))
        else
            t_off[midbits] = u_off[midbits]
        end

        diff_l == 0 || (siz_l_flg |= hibit; push!(sizvecl, base => diff_l))
        diff_u == 0 || (siz_u_flg |= hibit; push!(sizvecu, base => diff_u))

        if bitoff == 16
            # Reset bits
            pos  = (base >>> 9) + 1
            (is_l_tab[pos]  = _add_to_bv!(bitvec, bit_is_l)) == 0  || (is_l_flg |= hibit)
            (is_u_tab[pos]  = _add_to_bv!(bitvec, bit_is_u)) == 0  || (is_u_flg |= hibit)
            (can_l_tab[pos] = _add_to_bv!(bitvec, bit_can_l)) == 0 || (can_l_flg |= hibit)
            (can_u_tab[pos] = _add_to_bv!(bitvec, bit_can_u)) == 0 || (can_u_flg |= hibit)
        end

        # Check for end of chunk
        midbits == 32 || continue

        if l_mid
            push!(offvec, tuple(l_off...))
            l_tab[hipos] = UInt8(length(offvec))
            fill!(l_off, 0x0)
            l_mid = false
        end
        if u_mid
            push!(offvec, tuple(u_off...))
            u_tab[hipos] = UInt8(length(offvec))
            fill!(u_off, 0x0)
            u_mid = false
        end
        if t_mid
            push!(offvec, tuple(t_off...))
            t_tab[hipos] = UInt8(length(offvec))
            t_mid = false
        else
            t_tab[hipos] = u_tab[hipos]
        end
    end

    # Handle SLP
    can_sl_flg = can_su_flg = 0%UInt128

    for base in 0x0000:0x20:0xffe0
        hipos  = (base >>> 10) + 65
        hibit  = (1%UInt128) << ((base >>> 9) & 0x7f)
        midbits = ((base >>> 5) & 0x1f) + 1

        # Handle block of 32 characters
        l_flg = u_flg = false
        for off = 0x00:0x1f
            cp = UInt32(0x10000 + base + off)
            ch = Char(cp)
            cl = UInt32(lowercase(ch))
            cu = UInt32(uppercase(ch))
            ct = UInt32(titlecase(ch))
            lowbit = UInt32(1) << off
            islowercase(ch) && (tmp_is_l |= lowbit)
            isuppercase(ch) && (tmp_is_u |= lowbit)
            tmp_l[off+1] = cl%UInt16
            tmp_u[off+1] = cu%UInt16
            cp === cl === cu === ct && continue
            ct == cu || error("titlecase: $cp -> $ct not same as uppercase $cu")
            sizc = 4
            sizl = utf8_len(cl)
            sizu = utf8_len(cu)
            sizl == sizu == 4 || error("UTF-8 sizes not 4: $cp: $cl=$sizl, $cu=$sizu")
            if cl !== cp
                tmp_can_l |= lowbit
                l_flg = true
            end
            if cu !== cp
                tmp_can_u |= lowbit
                u_flg = true
            end
        end

        if l_flg
            l_mid = true
            push!(tupvec, tuple(tmp_l...))
            l_off[midbits] = UInt8(length(tupvec))
        end
        if u_flg
            u_mid = true
            push!(tupvec, tuple(tmp_u...))
            u_off[midbits] = UInt8(length(tupvec))
        end

        bitoff = ((base >>> 5) & 0xf) + 1
        bit_is_l[bitoff] = tmp_is_l
        bit_is_u[bitoff] = tmp_is_u
        bit_can_l[bitoff] = tmp_can_l
        bit_can_u[bitoff] = tmp_can_u

        tmp_is_l = tmp_is_u = tmp_can_l = tmp_can_u = 0%UInt32

        if bitoff == 16
            # Reset bits
            pos  = (base >>> 9) + 129
            (is_l_tab[pos]  = _add_to_bv!(bitvec, bit_is_l)) == 0  || (is_sl_flg |= hibit)
            (is_u_tab[pos]  = _add_to_bv!(bitvec, bit_is_u)) == 0  || (is_su_flg |= hibit)
            (can_l_tab[pos] = _add_to_bv!(bitvec, bit_can_l)) == 0 || (can_sl_flg |= hibit)
            (can_u_tab[pos] = _add_to_bv!(bitvec, bit_can_u)) == 0 || (can_su_flg |= hibit)
        end

        # Check for end of chunk
        midbits == 32 || continue

        if l_mid
            push!(offvec, tuple(l_off...))
            l_tab[hipos] = UInt8(length(offvec))
            fill!(l_off, 0x0)
            l_mid = false
        end
        if u_mid
            push!(offvec, tuple(u_off...))
            u_tab[hipos] = UInt8(length(offvec))
            fill!(u_off, 0x0)
            u_mid = false
        end
    end

    # Check that there are no upper / lower / title case characters above SLP
    for cp in 0x20000:0x10ffff
        ch = Char(cp)
        ch == lowercase(ch) == uppercase(ch) == titlecase(ch) ||
            error("$cp has lower/upper/titlecase")
    end

    (CaseTable(tuple(l_tab...), tuple(u_tab...), tuple(t_tab...),
               tuple(is_l_tab...), tuple(is_u_tab...),
               tuple(can_l_tab...), tuple(can_u_tab...),
               can_l_flg, can_u_flg, can_t_flg, can_sl_flg, can_su_flg,
               is_l_flg, is_u_flg, is_sl_flg, is_su_flg,
               siz_l_flg, siz_u_flg, max_siz_l, max_siz_u),
     tuple(tupvec...), tuple(offvec...), tuple(bitvec...),
     tuple(sizvecl...), tuple(sizvecu...))
end

const ct, tupvec, offvec, bitvec, sizvecl, sizvecu = case_tables()

end # module CaseTables

ext: jl


text:
# From @jbaron https://github.com/denizyuret/Knet.jl/issues/587

using Serialization

"""
Enable saving and loading of models by specialized KnetArray methods for Julia serialization
This will effectively move a GPU weight to the CPU before serializing it and move it back to
the GPU when deserializing.
"""
function Serialization.serialize(s::Serialization.AbstractSerializer, p::KnetArray)
    Serialization.serialize_type(s, typeof(p))
    Serialization.serialize(s, Array(p))
end

function Serialization.deserialize(s::Serialization.AbstractSerializer, t::Type{<:KnetArray})
    arr = Serialization.deserialize(s)
    return KnetArray(arr)
end


text:
###test local default functions
@test GeneralizedMetropolisHastings._defaultkey() == :param
@test GeneralizedMetropolisHastings._defaultvalue(Float64) == 0.0
@test GeneralizedMetropolisHastings._defaultvalue(Float32) == 0.0f0
@test GeneralizedMetropolisHastings._defaultvalue(Distributions.Uniform(0.0,11.0)) == 5.5
@test GeneralizedMetropolisHastings._defaultvalue(Distributions.DiscreteUniform(0,11)) == 5
@test GeneralizedMetropolisHastings._defaultunivariate(0.0,1.0) == Distributions.Uniform(0.0,1.0)
@test GeneralizedMetropolisHastings._defaultunivariate(0,1) == Distributions.DiscreteUniform(0,1)
@test_throws MethodError GeneralizedMetropolisHastings._defaultunivariate(0.0,1)

###test the parameter factory functions for ParameterDefault
d1 = parameter(:d,3.0)
d2 = parameter(:d,Float64)
d3 = parameter(:d,Float32)
d4 = parameter(3.0f0)
d5 = parameter(Float64)
d6 = parameter(Float32)
@test typeof(d1) <: ParameterDefault{Float64} && d1.key == :d && d1.default == 3.0
@test typeof(d2) <: ParameterDefault{Float64} && d2.key == :d && d2.default == 0.0
@test typeof(d3) <: ParameterDefault{Float32} && d3.key == :d && d3.default == 0.0f0
@test typeof(d4) <: ParameterDefault{Float32} && d4.default == 3.0f0
@test typeof(d5) <: ParameterDefault{Float64} && d5.default == 0.0
@test typeof(d6) <: ParameterDefault{Float32} && d6.default == 0.0
@test d2 == parameter(:d,Float64) #equality operator
@test d6 == parameter(Float32)
@test d1 != d2
@test d1 != d5

###test parameter factory functions for ParameterUnivariate
p1 = parameter(:p,Distributions.Normal(),3.0)
p2 = parameter(:p,Distributions.Binomial(),0)
p3 = parameter(:p,Distributions.Uniform(0.0,11.0))
p4 = parameter(:p,Distributions.DiscreteUniform(0,9))
p5 = parameter(:p,0.0,11.0,3.0)
p6 = parameter(:p,0.0,11.0)
p7 = parameter(Distributions.Uniform(0.0,11.0))
p8 = parameter(Distributions.DiscreteUniform(0,9))
p9 = parameter(0.0,11.0,3.0)
p10 = parameter(0.0,11.0)
@test typeof(p1) <: ParameterUnivariate && p1.key == :p && p1.default == 3.0 && p1.prior == Distributions.Normal()
@test typeof(p2) <: ParameterUnivariate && p2.key == :p && p2.default == 0 && p2.prior == Distributions.Binomial()
@test typeof(p3) <: ParameterUnivariate && p3.key == :p && p3.default == 5.5 && p3.prior == Distributions.Uniform(0.0,11.0)
@test typeof(p4) <: ParameterUnivariate && p4.key == :p && p4.default == 4 && p4.prior == Distributions.DiscreteUniform(0,9)
@test typeof(p5) <: ParameterUnivariate && p5.key == :p && p5.default == 3.0 && p5.prior == Distributions.Uniform(0.0,11.0)
@test typeof(p6) <: ParameterUnivariate && p6.key == :p && p6.default == 5.5 && p6.prior == Distributions.Uniform(0.0,11.0)
@test typeof(p7) <: ParameterUnivariate && p7.default == 5.5 && p7.prior == Distributions.Uniform(0.0,11.0)
@test typeof(p8) <: ParameterUnivariate && p8.default == 4 && p8.prior == Distributions.DiscreteUniform(0,9)
@test typeof(p9) <: ParameterUnivariate && p9.default == 3.0 && p9.prior == Distributions.Uniform(0.0,11.0)
@test typeof(p10) <: ParameterUnivariate && p10.default == 5.5 && p10.prior == Distributions.Uniform(0.0,11.0)
@test p3 == p6
@test p7 == p10
@test p1 != p2
@test p1 != p7
@test_throws MethodError parameter(:u,Distributions.Normal(),3)
@test_throws MethodError parameter(:u,Distributions.Binomial(),1.0)
@test_throws MethodError parameter(:u,1.0,2.0,3)
@test_throws AssertionError parameter(:u,Distributions.Uniform(1.0,2.0),0.0)
@test_throws MethodError parameter(Distributions.Normal(),3)
@test_throws MethodError parameter(Distributions.Binomial(),1.0)
@test_throws MethodError parameter(1.0,3)
@test_throws AssertionError parameter(Distributions.Uniform(1.0,2.0),0.0)

@test length(parameters([:d1,:d2],[3.0,1.0])) == 2
@test length(parameters([:d1,:d2],Float32)) == 2
@test length(parameters([1.0,2.0,3.0])) == 3
@test length(parameters(2,Float32)) == 2
vec1 = parameters([:d1,:d2],[Distributions.Normal(),Distributions.Uniform(0.0,11.0)],[0.0,2.0])
vec2 = parameters([:d1,:d2],[0.0,5.0],[6.0,10.0],[3.0,7.5])
vec3 = parameters([:d1,:d2],[Distributions.Normal(),Distributions.Uniform(0.0,11.0)])
vec4 = parameters([:d1,:d2],[0.0,5.0],[6.0,10.0])
vec5 = parameters([Distributions.Normal(),Distributions.Uniform(0.0,11.0)],[0.0,2.0])
vec6 = parameters([0.0,5.0],[6.0,10.0],[3.0,7.5])
vec7 = parameters([Distributions.Normal(),Distributions.Uniform(0.0,11.0)])
vec8 = parameters([0.0,5.0],[6.0,10.0])
vec9 = parameters([3.0,7.5])
@test vec2 == vec4
@test vec6 == vec8

###test the initvalues functions
v1 = GeneralizedMetropolisHastings._initvalue(Val{:default},d1,Float64)
v3 = GeneralizedMetropolisHastings._initvalue(Val{:default},d3,Float32)
v4 = GeneralizedMetropolisHastings._initvalue(Val{:default},d4,Float64)
v6 = GeneralizedMetropolisHastings._initvalue(Val{:default},d6,Float32)
@test typeof(v1) == Float64 && v1 == 3.0
@test typeof(v3) == Float32 && v3 == 0.0f0
@test typeof(v4) == Float64 && v4 == 3.0
@test typeof(v6) == Float32 && v6 == 0.0f0
@test GeneralizedMetropolisHastings._initvalue(Val{:default},p1,Float64) == 3.0
@test GeneralizedMetropolisHastings._initvalue(Val{:default},p3,Float64) == 5.5
@test GeneralizedMetropolisHastings._initvalue(Val{:default},p5,Float64) == 3.0
@test GeneralizedMetropolisHastings._initvalue(Val{:default},p8,Float64) == 4.0
@test (srand(0) ; GeneralizedMetropolisHastings._initvalue(Val{:prior},p3,Float64)) == (srand(0) ; rand(Distributions.Uniform(0.0,11.0)))
@test (srand(0) ; GeneralizedMetropolisHastings._initvalue(Val{:prior},p4,Float32)) == (srand(0) ; Float32(rand(Distributions.DiscreteUniform(0,9))))
@test (srand(0) ; GeneralizedMetropolisHastings._initvalue(Val{:prior},p4,Int)) == (srand(0) ; Int(rand(Distributions.DiscreteUniform(0,9))))
@test GeneralizedMetropolisHastings._initvalues!(Val{:default},AbstractParameter[d1,d3,p1,p4],Array{Float32}(4)) == [3.0f0,0.0f0,3.0f0,4.0f0]
@test initvalues!(trait(:initialize,:default),AbstractParameter[d1,d3,p1,p4],Array{Float64}(4)) == [3.0,0.0,3.0,4.0]
@test initvalues(trait(:initialize,:default),AbstractParameter[d1,d3,p1,p4],Float64) == [3.0,0.0,3.0,4.0]
@test (srand(0) ; initvalues(trait(:initialize,:prior),AbstractParameter[d1,d3,p1,p4],Float64)) == (srand(0) ; [3.0,0.0,rand(Distributions.Normal()),Float64(rand(Distributions.DiscreteUniform(0,9)))])
@test_throws AssertionError initvalues!(trait(:initialize,:default),AbstractParameter[p1],Array{Float64}(2))

###test the logprior functions
l1 = parameter(Distributions.Uniform(1.0,3.0))
l2 = parameter(2.0)
l3 = parameter(Distributions.DiscreteUniform(1.0,3.0))
@test GeneralizedMetropolisHastings._logprior(l1,2.0,Float64) == log(1.0/2.0)
@test GeneralizedMetropolisHastings._logprior(l1,4.0,Float64) == -Inf
@test GeneralizedMetropolisHastings._logprior(l2,2.0,Float64) == 0.0
@test GeneralizedMetropolisHastings._logprior(l1,2.0,Float32) == -log(2.0f0)
@test GeneralizedMetropolisHastings._logprior(l1,4.0,Float32) == -Inf32
@test GeneralizedMetropolisHastings._logprior(l2,2.0,Float32) == 0.0f0
@test GeneralizedMetropolisHastings._logprior(l3,2.0,Float64) == log(1.0/3.0)
@test GeneralizedMetropolisHastings._logprior([l1,l2,l3],[2.0,2.0,2.0],Float64) == log(1.0/2.0) + log(1.0/3.0)
@test logprior(l1,2.0,Float64) == log(1.0/2.0)
@test logprior(AbstractParameter[l1,l2,l3],[2.0,2.0,2.0],Float64) == log(1.0/2.0) + log(1.0/3.0)
@test logprior(AbstractParameter[l1,l1],[2.0,4.0],Float64) == -Inf
@test logprior(AbstractParameter[],Float64[],Float64) == 0.0
@test logprior(ParameterDefault[l2,l2],[2.0,4.0],Float64) == 0.0
@test logprior(ParameterUnivariate[l1,l3],[2.0,2.0],Float64) == log(1.0/2.0) + log(1.0/3.0)
@test logprior!(zeros(2),[l1,l2,l3],[2.0 2.0;2.0 2.0;2.0 4.0]) == [log(1.0/2.0)+log(1.0/3.0),-Inf]
@test logprior([l1,l2,l3],[2.0 2.0;2.0 2.0;2.0 4.0],Float64) == [log(1.0/2.0)+log(1.0/3.0),-Inf]
@test logprior(ParameterDefault[l2,l2],[2.0 2.0;2.0 2.0],Float64) == zeros(2)

println("====================")
println("Test show() function")
println()
show(d1)
show(d4)
show(p1)
show(p3)
show(p4)
show(AbstractParameter[l1,l2,l3])
show(ParameterUnivariate[l1,l3])
show(vec1)
show(vec4)
show([vec1;vec4])
show([vec1;vec9])
println()
println("End  show() function")
println("====================")
println()

nothing


text:
@with_kw mutable struct LShapedData{T <: AbstractFloat}
    Q::T = 1e10
    θ::T = -1e10
    ncuts::Int = 0
    iterations::Int = 1
    consolidations::Int = 0
end

@with_kw mutable struct LShapedParameters{T <: AbstractFloat}
    τ::T = 1e-6
    cut_scaling::T = 1.0
    debug::Bool = false
    log::Bool = true
    keep::Bool = true
    offset::Int = 0
    indent::Int = 0
end

"""
    LShaped

Functor object for the L-shaped algorithm. Create using the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Algorithm parameters
- `τ::AbstractFloat = 1e-6`: Relative tolerance for convergence checks.
- `debug::Bool = false`: Specifies if extra information should be saved for debugging purposes. Defaults to false for memory efficiency.
- `cut_scaling::AbstractFloat = 1.0`: Rescaling factor for cutting planes to improve numerical stability.
- `log::Bool = true`: Specifices if L-shaped procedure should be logged on standard output or not.
...
"""
struct LShaped{T <: AbstractFloat,
               A <: AbstractVector,
               SP <: StochasticProgram,
               M <: LQSolver,
               S <: LQSolver,
               E <: AbstractExecution,
               F <: AbstractFeasibility,
               R <: AbstractRegularization,
               Agg <: AbstractAggregation,
               C <: AbstractConsolidation} <: AbstractLShapedSolver
    stochasticprogram::SP
    data::LShapedData{T}
    parameters::LShapedParameters{T}

    # Master
    mastersolver::M
    mastervector::A
    x::A
    Q_history::A

    # Subproblems
    nscenarios::Int

    # Execution
    execution::E

    # Feasibility
    feasibility::F

    # Regularization
    regularization::R

    # Cuts
    θs::A
    cuts::Vector{AnySparseOptimalityCut{T}}
    aggregation::Agg
    consolidation::C
    θ_history::A

    progress::ProgressThresh{T}

    function LShaped(stochasticprogram::StochasticProgram,
                     x₀::AbstractVector,
                     mastersolver::MPB.AbstractMathProgSolver,
                     subsolver::MPB.AbstractMathProgSolver,
                     feasibility_cuts::Bool,
                     executer::Execution,
                     regularizer::AbstractRegularizer,
                     aggregator::AbstractAggregator,
                     consolidator::AbstractConsolidator; kw...)
        if nworkers() > 1 && executer isa Serial
            @warn "There are worker processes, consider using distributed version of algorithm"
        end
        executer = if nworkers() == 1 && !(executer isa Serial)
            @warn "There are no worker processes, defaulting to serial version of algorithm"
            Serial()
        else
            executer
        end
        first_stage = StochasticPrograms.get_stage_one(stochasticprogram)
        length(x₀) != first_stage.numCols && error("Incorrect length of starting guess, has ", length(x₀), " should be ", first_stage.numCols)

        T = promote_type(eltype(x₀), Float32)
        x₀_ = convert(AbstractVector{T}, copy(x₀))
        mastervector = convert(AbstractVector{T}, copy(x₀))
        A = typeof(x₀_)
        SP = typeof(stochasticprogram)
        msolver = LQSolver(first_stage, mastersolver)
        M = typeof(msolver)
        S = LQSolver{typeof(MPB.LinearQuadraticModel(subsolver)),typeof(subsolver)}
        n = StochasticPrograms.nscenarios(stochasticprogram)
        # Feasibility
        feasibility = feasibility_cuts ? HandleFeasibility(T) : IgnoreFeasibility()
        F = typeof(feasibility)
        # Execution policy
        execution = executer(n,F,T,A,S)
        E = typeof(execution)
        # Regularization policy
        regularization = regularizer(x₀_)
        R = typeof(regularization)
        # Aggregation policy
        aggregation = aggregator(n, T)
        Agg = typeof(aggregation)
        # Consolidation policy
        consolidation = consolidator(T)
        C = typeof(consolidation)
        # Algorithm parameters
        params = LShapedParameters{T}(; kw...)

        lshaped = new{T,A,SP,M,S,E,F,R,Agg,C}(stochasticprogram,
                                              LShapedData{T}(),
                                              params,
                                              msolver,
                                              mastervector,
                                              x₀_,
                                              A(),
                                              n,
                                              execution,
                                              feasibility,
                                              regularization,
                                              A(),
                                              Vector{SparseOptimalityCut{T}}(),
                                              aggregation,
                                              consolidation,
                                              A(),
                                              ProgressThresh(T(1.0), 0.0, "$(indentstr(params.indent))L-Shaped Gap "))
        # Initialize solver
        initialize!(lshaped, subsolver)
        return lshaped
    end
end
LShaped(stochasticprogram::StochasticProgram,
        mastersolver::MPB.AbstractMathProgSolver,
        subsolver::MPB.AbstractMathProgSolver,
        feasibility_cuts::Bool,
        executer::Execution,
        regularizer::AbstractRegularizer,
        aggregator::AbstractAggregator,
        consolidator::AbstractConsolidator; kw...) = LShaped(stochasticprogram,
                                                             rand(decision_length(stochasticprogram)),
                                                             mastersolver,
                                                             subsolver,
                                                             feasibility_cuts,
                                                             executer,
                                                             regularizer,
                                                             aggregator,
                                                             consolidator; kw...)

function show(io::IO, lshaped::LShaped)
    println(io, typeof(lshaped).name.name)
    println(io, "State:")
    show(io, lshaped.data)
    println(io, "Parameters:")
    show(io, lshaped.parameters)
end

function show(io::IO, ::MIME"text/plain", lshaped::LShaped)
    show(io, lshaped)
end

function (lshaped::LShaped)()
    # Reset timer
    lshaped.progress.tfirst = lshaped.progress.tlast = time()
    # Start workers (if any)
    start_workers!(lshaped)
    # Start procedure
    while true
        status = iterate!(lshaped)
        if status != :Valid
            close_workers!(lshaped)
            return status
        end
    end
end


text:
type NFib{T<:Integer}
    n::T
    klim::T
    seeder::Function
end

type FState
    a::Array{BigInt,1}
    adex::Integer
    k::Integer
end

function Base.start{T<:Integer}(nf::NFib{T})
    a = nf.seeder(nf.n)
    adex = 1
    k = 1
    return FState(a, adex, k)
end

function Base.done{T<:Integer}(nf::NFib{T}, fs::FState)
    fs.k > nf.klim
end

function Base.next{T<:Integer}(nf::NFib{T}, fs::FState)
    f = sum(fs.a)
    fs.a[fs.adex] = f
    fs.adex = rem1(fs.adex+1, nf.n)
    fs.k += 1
    return (f, fs)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Compat
using Test
using UUIDs: UUID, uuid1, uuid_version

@test isempty(detect_ambiguities(Base, Core, Compat))

@testset "CartesianIndex" begin
    # https://github.com/JuliaLang/julia/pull/29440
    ci = CartesianIndex(1, 1)
    @test length(-ci:ci) == 9
    # https://github.com/JuliaLang/julia/pull/29442
    @test oneunit(ci) === ci
    # https://github.com/JuliaLang/julia/pull/30268
    A = randn(1,2,3)
    @test get(A, CartesianIndex(1,2,3), :some_default) === A[1,2,3]
    @test get(A, CartesianIndex(2,2,3), :some_default) === :some_default
    @test get(11:15, CartesianIndex(6), nothing) === nothing
    @test get(11:15, CartesianIndex(5), nothing) === 15
end

# julia#29679
@test !isnothing(1)
@test isnothing(nothing)

# https://github.com/JuliaLang/julia/pull/29749
@testset "row/column/slice iterators" begin
    # Simple ones
    M = [1 2 3; 4 5 6; 7 8 9]
    @test collect(eachrow(M)) == collect(eachslice(M, dims = 1)) == [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
    @test collect(eachcol(M)) == collect(eachslice(M, dims = 2)) == [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
    @test_throws DimensionMismatch eachslice(M, dims = 4)

    # Higher-dimensional case
    M = reshape([(1:16)...], 2, 2, 2, 2)
    @test_throws MethodError collect(eachrow(M))
    @test_throws MethodError collect(eachcol(M))
    @test collect(eachslice(M, dims = 1))[1][:, :, 1] == [1 5; 3 7]
end

# Support for positional `stop`
@test range(0, 5, length = 6) == 0.0:1.0:5.0
@test range(0, 10, step = 2) == 0:2:10
@test_throws ArgumentError range(0, 10)

mutable struct TLayout
    x::Int8
    y::Int16
    z::Int32
end
tlayout = TLayout(5,7,11)
@test hasfield(TLayout, :y)
@test !hasfield(TLayout, :a)
@test hasproperty(tlayout, :x)
@test !hasproperty(tlayout, :p)

@test merge((a=1,b=1)) == (a=1,b=1)
@test merge((a=1,), (b=2,), (c=3,)) == (a=1,b=2,c=3)

@testset "only" begin
    @test only([3]) === 3
    @test_throws ArgumentError only([])
    @test_throws ArgumentError only([3, 2])

    @test @inferred(only((3,))) === 3
    @test_throws ArgumentError only(())
    @test_throws ArgumentError only((3, 2))

    @test only(Dict(1=>3)) === (1=>3)
    @test_throws ArgumentError only(Dict{Int,Int}())
    @test_throws ArgumentError only(Dict(1=>3, 2=>2))

    @test only(Set([3])) === 3
    @test_throws ArgumentError only(Set(Int[]))
    @test_throws ArgumentError only(Set([3,2]))

    @test @inferred(only((;a=1))) === 1
    @test_throws ArgumentError only(NamedTuple())
    @test_throws ArgumentError only((a=3, b=2.0))

    @test @inferred(only(1)) === 1
    @test @inferred(only('a')) === 'a'
    @test @inferred(only(Ref([1, 2]))) == [1, 2]
    @test_throws ArgumentError only(Pair(10, 20))

    @test only(1 for ii in 1:1) === 1
    @test only(1 for ii in 1:10 if ii < 2) === 1
    @test_throws ArgumentError only(1 for ii in 1:10)
    @test_throws ArgumentError only(1 for ii in 1:10 if ii > 2)
    @test_throws ArgumentError only(1 for ii in 1:10 if ii > 200)
end

# https://github.com/JuliaLang/julia/pull/32628
@testset "mod with ranges" begin
    for n in -10:10
        @test mod(n, 0:4) == mod(n, 5)
        @test mod(n, 1:5) == mod1(n, 5)
        @test mod(n, 2:6) == 2 + mod(n-2, 5)
        @test mod(n, Base.OneTo(5)) == mod1(n, 5)
    end
    @test mod(Int32(3), 1:5) == 3
    @test mod(big(typemax(Int))+99, 0:4) == mod(big(typemax(Int))+99, 5)
    @test_throws MethodError mod(3.141, 1:5)
    @test_throws MethodError mod(3, UnitRange(1.0,5.0))
    @test_throws MethodError mod(3, 1:2:7)
    @test_throws DivideError mod(3, 1:0)
end

using LinearAlgebra

@testset "generalized dot #32739" begin
    # stdlib/LinearAlgebra/test/generic.jl
    for elty in (Int, Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat})
        n = 10
        if elty <: Int
            A = rand(-n:n, n, n)
            x = rand(-n:n, n)
            y = rand(-n:n, n)
        elseif elty <: Real
            A = convert(Matrix{elty}, randn(n,n))
            x = rand(elty, n)
            y = rand(elty, n)
        else
            A = convert(Matrix{elty}, complex.(randn(n,n), randn(n,n)))
            x = rand(elty, n)
            y = rand(elty, n)
        end
        @test dot(x, A, y) ≈ dot(A'x, y) ≈ *(x', A, y) ≈ (x'A)*y
        @test dot(x, A', y) ≈ dot(A*x, y) ≈ *(x', A', y) ≈ (x'A')*y
        elty <: Real && @test dot(x, transpose(A), y) ≈ dot(x, transpose(A)*y) ≈ *(x', transpose(A), y) ≈ (x'*transpose(A))*y
        B = reshape([A], 1, 1)
        x = [x]
        y = [y]
        @test dot(x, B, y) ≈ dot(B'x, y)
        @test dot(x, B', y) ≈ dot(B*x, y)
        elty <: Real && @test dot(x, transpose(B), y) ≈ dot(x, transpose(B)*y)
    end

    # stdlib/LinearAlgebra/test/symmetric.jl
    n = 10
    areal = randn(n,n)/2
    aimg  = randn(n,n)/2
    @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        asym = transpose(a) + a                 # symmetric indefinite
        aherm = a' + a                 # Hermitian indefinite
        apos  = a' * a                 # Hermitian positive definite
        aposs = apos + transpose(apos)        # Symmetric positive definite
        ε = εa = eps(abs(float(one(eltya))))
        x = randn(n)
        y = randn(n)
        b = randn(n,n)/2
        x = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(x, zeros(n)) : x)
        y = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(y, zeros(n)) : y)
        b = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(b, zeros(n,n)) : b)

        @testset "generalized dot product" begin
            for uplo in (:U, :L)
                @test dot(x, Hermitian(aherm, uplo), y) ≈ dot(x, Hermitian(aherm, uplo)*y) ≈ dot(x, Matrix(Hermitian(aherm, uplo)), y)
                @test dot(x, Hermitian(aherm, uplo), x) ≈ dot(x, Hermitian(aherm, uplo)*x) ≈ dot(x, Matrix(Hermitian(aherm, uplo)), x)
            end
            if eltya <: Real
                for uplo in (:U, :L)
                    @test dot(x, Symmetric(aherm, uplo), y) ≈ dot(x, Symmetric(aherm, uplo)*y) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), y)
                    @test dot(x, Symmetric(aherm, uplo), x) ≈ dot(x, Symmetric(aherm, uplo)*x) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), x)
                end
            end
        end
    end

    # stdlib/LinearAlgebra/test/uniformscaling.jl
    @testset "generalized dot" begin
        x = rand(-10:10, 3)
        y = rand(-10:10, 3)
        λ = rand(-10:10)
        J = UniformScaling(λ)
        @test dot(x, J, y) == λ*dot(x, y)
    end

    # stdlib/LinearAlgebra/test/bidiag.jl
    # The special method for this is not in Compat #683, so this tests the generic fallback
    @testset "generalized dot" begin
        for elty in (Float64, ComplexF64)
            dv = randn(elty, 5)
            ev = randn(elty, 4)
            x = randn(elty, 5)
            y = randn(elty, 5)
            for uplo in (:U, :L)
                B = Bidiagonal(dv, ev, uplo)
                @test dot(x, B, y) ≈ dot(B'x, y) ≈ dot(x, Matrix(B), y)
            end
        end
    end

    # Diagonal -- no such test in Base.
    @testset "diagonal" begin
        x = rand(-10:10, 3) .+ im
        y = rand(-10:10, 3) .+ im
        d = Diagonal(rand(-10:10, 3) .+ im)
        @test dot(x,d,y) == dot(x,collect(d),y) == dot(x, d*y)
    end
end

# https://github.com/JuliaLang/julia/pull/33568
@testset "function composition" begin
    @test ∘(x -> x-2, x -> x-3, x -> x+5)(7) == 7
    fs = [x -> x[1:2], uppercase, lowercase]
    @test ∘(fs...)("ABC") == "AB"

    # https://github.com/JuliaLang/julia/pull/34251
    @testset "unary" begin
        @test ∘(identity) === identity
        @test ∘(inv) === inv
    end
end

# https://github.com/JuliaLang/julia/pull/33128
@testset "pkgdir" begin
    @test pkgdir(Main) === nothing
    @test joinpath(pkgdir(Compat), "") == abspath(joinpath(@__DIR__, ".."))
end

# https://github.com/JuliaLang/julia/pull/33736/
@testset "ReverseOrdering constructor" begin
    @test Base.Order.ReverseOrdering() == Base.Order.Reverse
end

# https://github.com/JuliaLang/julia/pull/32968
@testset "filter on Tuples" begin
    @test filter(isodd, (1,2,3)) == (1, 3)
    @test filter(isequal(2), (true, 2.0, 3)) === (2.0,)
    @test filter(i -> true, ()) == ()
    @test filter(identity, (true,)) === (true,)
    longtuple = ntuple(identity, 20)
    @test filter(iseven, longtuple) == ntuple(i->2i, 10)
    @test filter(x -> x<2, (longtuple..., 1.5)) === (1, 1.5)
end

# https://github.com/JuliaLang/julia/pull/34652
@testset "ismutable" begin
    @test ismutable(1) == false
    @test ismutable([]) == true
end

# https://github.com/JuliaLang/julia/pull/28761
@testset "uuid5" begin
    u1 = uuid1()
    u5 = uuid5(u1, "julia")
    @test uuid_version(u5) == 5
    @test u5 == UUID(string(u5)) == UUID(GenericString(string(u5)))
    @test u5 == UUID(UInt128(u5))

    following_uuids = [
        UUID("22b4a8a1-e548-4eeb-9270-60426d66a48e"),
        UUID("30ea6cfd-c270-569f-b4cb-795dead63686"),
        UUID("31099374-e3a0-5fde-9482-791c639bf29b"),
        UUID("6b34b357-a348-53aa-8c71-fb9b06c3a51e"),
        UUID("fdbd7d4d-c462-59cc-ae6a-0c3b010240e2"),
        UUID("d8cc6298-75d5-57e0-996c-279259ab365c"),
    ]

    for (idx, init_uuid) in enumerate(following_uuids[1:end-1])
        next_id = uuid5(init_uuid, "julia")
        @test next_id == following_uuids[idx+1]
    end

    # Some UUID namespaces provided in the appendix of RFC 4122
    # https://tools.ietf.org/html/rfc4122.html#appendix-C
    namespace_dns  = UUID(0x6ba7b8109dad11d180b400c04fd430c8) # 6ba7b810-9dad-11d1-80b4-00c04fd430c8
    namespace_url  = UUID(0x6ba7b8119dad11d180b400c04fd430c8) # 6ba7b811-9dad-11d1-80b4-00c04fd430c8
    namespace_oid  = UUID(0x6ba7b8129dad11d180b400c04fd430c8) # 6ba7b812-9dad-11d1-80b4-00c04fd430c8
    namespace_x500 = UUID(0x6ba7b8149dad11d180b400c04fd430c8) # 6ba7b814-9dad-11d1-80b4-00c04fd430c8

    # Python-generated UUID following each of the standard namespaces
    standard_namespace_uuids = [
        (namespace_dns,  UUID("00ca23ad-40ef-500c-a910-157de3950d07")),
        (namespace_oid,  UUID("b7bf72b0-fb4e-538b-952a-3be296f07f6d")),
        (namespace_url,  UUID("997cd5be-4705-5439-9fe6-d77b18d612e5")),
        (namespace_x500, UUID("993c6684-82e7-5cdb-bd46-9bff0362e6a9")),
    ]

    for (init_uuid, next_uuid) in standard_namespace_uuids
        result = uuid5(init_uuid, "julia")
        @test next_uuid == result
    end
end

@testset "Irrational zero and one" begin
    @test one(pi) === true
    @test zero(pi) === false
    @test one(typeof(pi)) === true
    @test zero(typeof(pi)) === false
end

# https://github.com/JuliaLang/julia/pull/32753
@testset "evalpoly real" begin
    for x in -1.0:2.0, p1 in -3.0:3.0, p2 in -3.0:3.0, p3 in -3.0:3.0
        evpm = @evalpoly(x, p1, p2, p3)
        @test evalpoly(x, (p1, p2, p3)) == evpm
        @test evalpoly(x, [p1, p2, p3]) == evpm
    end
end
@testset "evalpoly complex" begin
    for x in -1.0:2.0, y in -1.0:2.0, p1 in -3.0:3.0, p2 in -3.0:3.0, p3 in -3.0:3.0
        z = x + im * y
        evpm = @evalpoly(z, p1, p2, p3)
        @test evalpoly(z, (p1, p2, p3)) == evpm
        @test evalpoly(z, [p1, p2, p3]) == evpm
    end
    @test evalpoly(1+im, (2,)) == 2
    @test evalpoly(1+im, [2,]) == 2
end

# https://github.com/JuliaLang/julia/pull/35298
begin
    # A custom linear slow sparse-like array that relies upon Dict for its storage
    struct TSlow{T,N} <: AbstractArray{T,N}
        data::Dict{NTuple{N,Int}, T}
        dims::NTuple{N,Int}
    end
    TSlow(::Type{T}, dims::Int...) where {T} = TSlow(T, dims)
    TSlow(::Type{T}, dims::NTuple{N,Int}) where {T,N} = TSlow{T,N}(Dict{NTuple{N,Int}, T}(), dims)

    TSlow{T,N}(X::TSlow{T,N})         where {T,N  } = X
    TSlow(     X::AbstractArray{T,N}) where {T,N  } = TSlow{T,N}(X)
    TSlow{T  }(X::AbstractArray{_,N}) where {T,N,_} = TSlow{T,N}(X)
    TSlow{T,N}(X::AbstractArray     ) where {T,N  } = begin
        A = TSlow(T, size(X))
        for I in CartesianIndices(X)
            A[Tuple(I)...] = X[Tuple(I)...]
        end
        A
    end
    Base.size(A::TSlow) = A.dims
    Base.similar(A::TSlow, ::Type{T}, dims::Dims) where {T} = TSlow(T, dims)
    Base.IndexStyle(::Type{A}) where {A<:TSlow} = IndexCartesian()
    Base.getindex(A::TSlow{T,N}, i::Vararg{Int,N}) where {T,N} = get(A.data, i, zero(T))
    Base.setindex!(A::TSlow{T,N}, v, i::Vararg{Int,N}) where {T,N} = (A.data[i] = v)
end

# https://github.com/JuliaLang/julia/pull/35304
@testset "similar(PermutedDimsArray)" begin
    x = PermutedDimsArray([1 2; 3 4], (2, 1))
    @test similar(x, 3,3) isa Array
    z = TSlow([1 2; 3 4])
    x_slow = PermutedDimsArray(z, (2, 1))
    @test similar(x_slow, 3,3) isa TSlow
end

# https://github.com/JuliaLang/julia/pull/34548
@testset "@NamedTuple" begin
    @test (@NamedTuple {a::Int, b::String}) === NamedTuple{(:a, :b),Tuple{Int,String}} ===
        @NamedTuple begin
            a::Int
            b::String
        end
    @test (@NamedTuple {a::Int, b}) === NamedTuple{(:a, :b),Tuple{Int,Any}}
end

struct NonFunctionCallable end
(::NonFunctionCallable)(args...) = +(args...)

@testset "mergewith" begin
    d1 = Dict("A" => 1, "B" => 2)
    d2 = Dict("B" => 3.0, "C" => 4.0)
    @test mergewith(+, d1, d2) == Dict("A" => 1, "B" => 5, "C" => 4)
    @test mergewith(*, d1, d2) == Dict("A" => 1, "B" => 6, "C" => 4)
    @test mergewith(-, d1, d2) == Dict("A" => 1, "B" => -1, "C" => 4)
    @test mergewith(NonFunctionCallable(), d1, d2) == Dict("A" => 1, "B" => 5, "C" => 4)
    @test foldl(mergewith(+), [d1, d2]; init=Dict{Union{},Union{}}()) ==
        Dict("A" => 1, "B" => 5, "C" => 4)
end

@testset "mergewith!" begin
    d1 = Dict("A" => 1, "B" => 3, "C" => 4)
    d2 = Dict("B" => 3, "C" => 4)
    mergewith!(+, d1, d2)
    @test d1 == Dict("A" => 1, "B" => 6, "C" => 8)
    mergewith!(*, d1, d2)
    @test d1 == Dict("A" => 1, "B" => 18, "C" => 32)
    mergewith!(-, d1, d2)
    @test d1 == Dict("A" => 1, "B" => 15, "C" => 28)
    mergewith!(NonFunctionCallable(), d1, d2)
    @test d1 == Dict("A" => 1, "B" => 18, "C" => 32)
    @test foldl(mergewith!(+), [d1, d2]; init=empty(d1)) ==
        Dict("A" => 1, "B" => 21, "C" => 36)
end

# https://github.com/JuliaLang/julia/pull/34427
@testset "isdisjoint" begin
    for S in (Set, BitSet, Vector)
        for (l,r) in ((S([1,2]),     S([3,4])),
                      (S([5,6,7,8]), S([7,8,9])),
                      (S([1,2]),     S([3,4])),
                      (S([5,6,7,8]), S([7,8,9])),
                      (S([1,2,3]),   S()),
                      (S(),          S()),
                      (S(),          S([1,2,3])),
                      (S([1,2,3]),   S([1])),
                      (S([1,2,3]),   S([1,2])),
                      (S([1,2,3]),   S([1,2,3])),
                      (S([1,2,3]),   S([4])),
                      (S([1,2,3]),   S([4,1])))
            @test isdisjoint(l,l) == isempty(l)
            @test isdisjoint(l,r) == isempty(intersect(l,r))
        end
    end
end

# https://github.com/JuliaLang/julia/pull/35577
@testset "union on OneTo" begin
    @test union(Base.OneTo(3), Base.OneTo(4)) === Base.OneTo(4)
end

nothing


text:
# Predefined quantities

const qₑ = -1.60217662e-19  # electron charge, [C]
const mₑ = 9.10938356e-31   # electron mass, [kg]
const qᵢ = 1.60217662e-19   # proton mass, [C]
const mᵢ = 1.673557546e-27  # proton mass, [kg]
const c  = 299792458.       # speed of light, [m/s]
const μ₀ = 4π*1e-7          # Vacuum permeability, [H/m]
const ϵ₀ = 1/(c^2*μ₀)       # Vacuum permittivity, [F/m]
const kB = 1.38064852e-23   # Boltzmann constant, [m²kg/(s²K)]
const Re = 6.371e6          # Earth radius, [m]

# Define units, LaTeX markup names, and LaTeX markup units for intrinsic values
const units_predefined = Dict(
   :rhom => ("kg/m3", L"$\rho_m$",  L"$\mathrm{kg}\,\mathrm{m}^{-3}$"),
   :rhoq => ("C/m3", L"$\rho_q$", L"$\mathrm{C}\,\mathrm{m}^{-3}$"),
   :rho => ("1/m3", L"$n_\mathrm{p}$", L"$\mathrm{m}^{-3}$"),
   :rhobackstream => ("1/m3",  L"$n_\mathrm{p,st}$", L"$\mathrm{m}^{-3}$"),
   :rhononbackstream => ("1/m3", L"$n_\mathrm{p,th}$", L"$\mathrm{m}^{-3}$"),
   :rho_v => ("1/m2s", L"$\Gamma_\mathrm{p}$", L"$\mathrm{m}^{-2}$s"),
   :rhovbackstream => ("1/m2s", L"$\Gamma_\mathrm{p,st}$", L"$\mathrm{m}^{-2}$s"),
   :rhovnonbackstream => ("1/m2s", L"$\Gamma_\mathrm{p,th}$", L"$\mathrm{m}^{-2}$s"),
   :v => ("m/s", L"$V$", L"$\mathrm{m}\,\mathrm{s}^{-1}$"),
   :vbackstream => ("m/s", L"$V_\mathrm{p,st}$", L"$\mathrm{m}\,\mathrm{s}^{-1}$"),
   :vNonbackstream => ("m/s", L"$V_\mathrm{p,th}$", L"$\mathrm{m}\,\mathrm{s}^{-1}$"),
   :b => ("T", L"$B$", "T"),
   :b_vol => ("T",  L"$B_\mathrm{vol}$", "T"),
   :background_b => ("T",  L"$B_\mathrm{bg}$", "T"),
   :perturbed_b => ("T", L"$B_\mathrm{pert}$", "T"),
   :bgb => ("T", L"$B_\mathrm{bg}$", "T"),
   :perb => ("T", L"B_\mathrm{pert}$", "T"),
   :perb_vol => ("T", L"B_\mathrm{vol,pert}$", "T"),
   :e => ("V/m", L"$E$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :e_vol => ("V/m",  L"$E_\mathrm{vol}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :exhall_000_100 => ("V/m", L"$E_\mathrm{Hall,000,100}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :exhall_001_101 => ("V/m", L"$E_\mathrm{Hall,001,101}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :exhall_010_110 => ("V/m", L"$E_\mathrm{Hall,010,110}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :exhall_011_111 => ("V/m", L"$E_\mathrm{Hall,011,111}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :eyhall_000_010 => ("V/m", L"$E_\mathrm{Hall,000,010}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :eyhall_001_011 => ("V/m", L"$E_\mathrm{Hall,001,011}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :eyhall_100_110 => ("V/m", L"$E_\mathrm{Hall,100,110}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :eyhall_101_111 => ("V/m", L"$E_\mathrm{Hall,101,111}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :ezhall_000_001 => ("V/m", L"$E_\mathrm{Hall,000,001}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :ezhall_010_011 => ("V/m", L"$E_\mathrm{Hall,010,011}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :ezhall_100_101 => ("V/m", L"$E_\mathrm{Hall,100,101}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :ezhall_110_111 => ("V/m", L"$E_\mathrm{Hall,110,111}$", L"$\mathrm{V}\,\mathrm{m}^{-1}$"),
   :pressure => ("Pa",  L"$P$", "Pa"),
   :pressure_dt2 => ("Pa", L"$P_{\mathrm{d}t/2}}$", "Pa"),
   :pressure_r => ("Pa",  L"$P_r$", "Pa"),
   :pressure_v => ("Pa", L"$P_v$", "Pa"),
   :ptensordiagonal => ("Pa", L"$\mathcal{P}_\mathrm{diag}$", "Pa"),
   :ptensoroffdiagonal => ("Pa", L"$\mathcal{P}_\mathrm{off-diag}$", "Pa"),
   :ptensorbackstreamdiagonal => ("Pa", L"$\mathcal{P}_\mathrm{st,diag}$", "Pa"),
   :ptensorbackstreamoffdiagonal => ("Pa",  L"$\mathcal{P}_\mathrm{st,off-diag}$", "Pa"),
   :ptensornonbackstreamdiagonal => ("Pa", L"$\mathcal{P}_\mathrm{th,diag}$", "Pa"),
   :ptensornonbackstreamoffdiagonal => ("Pa", L"$\mathcal{P}_\mathrm{th,off-diag}$", "Pa"),
   :t => ("K", L"$T$", "K"),
   :tpar => ("K", L"$T$", "K"),
   :tperp => ("K", L"$T$", "K"),
   :max_v_dt => ("s", L"$\Delta t_{\mathrm{max},v}$", "s"),
   :max_r_dt => ("s", L"$\Delta t_{\mathrm{max},r}$", "s"),
   :max_fields_dt => ("s", L"$\Delta t_\mathrm{max,FS}$", "s"),
   :minvalue => ("s3/m6", L"$f_\mathrm{Min}$", L"$\mathrm{m}^{-6}\,\mathrm{s}^{3}$"),
   :effectivesparsitythreshold => ("s3/m6", L"$f_\mathrm{Min}$", L"$\mathrm{m}^{-6}\,\mathrm{s}^{3}$"),
   :rho_loss_adjust => ("1/m3", L"$\Delta_\mathrm{loss} n_\mathrm{p}$", L"$\mathrm{m}^{-3}$"),
   :energydensity => ("eV/cm3", L"$\rho_{\mathrm{energy}}$", L"$\mathrm{eV}\,\mathrm{cm}^{-3}$"),
   :precipitationdiffflux => ("1/(cm2 sr s eV)", L"$'Delta F_\mathrm{precipitation}$", L"$\mathrm{cm}^{-2} \,\mathrm{sr}^{-1}\,\mathrm{s}^{-1}\,\mathrm{eV}^{-1}$"),
)

# Define derived parameters
const variables_predefined = Dict(
   :Bmag => function (meta, ids=UInt64[])
      rho_ = findfirst(endswith("rho"), meta.variable)
      ρ = isempty(ids) ?
         readvariable(meta, meta.variable[rho_]) :
         readvariable(meta, meta.variable[rho_], ids)
      if hasvariable(meta, "vg_b_vol")
         Bmag = isempty(ids) ?
            sqrt.(sum(readvariable(meta, "vg_b_vol").^2, dims=1)) :
            sqrt.(sum(readvariable(meta, "vg_b_vol", ids).^2, dims=1))
      else
         @assert isempty(ids) "Do not support reading selected cells from FSGrid!"
         Bmag = sqrt.(sum(readvariable(meta, "fg_b").^2, dims=1))
      end
      @inbounds for i = eachindex(ρ) # sparsity/inner boundary
         Bmag[i] == 0.0 && (Bmag[i] = NaN)
      end
      Bmag
   end,
   :Emag => function (meta, ids=UInt64[])
      rho_ = findfirst(endswith("rho"), meta.variable)
      ρ = isempty(ids) ?
         readvariable(meta, meta.variable[rho_]) :
         readvariable(meta, meta.variable[rho_], ids)
      if hasvariable(meta, "vg_e_vol")
         Emag = isempty(ids) ?
            sqrt.(sum(readvariable(meta, "vg_e_vol").^2, dims=1)) :
            sqrt.(sum(readvariable(meta, "vg_e_vol", ids).^2, dims=1))
      else
         @assert isempty(ids) "Do not support reading selected cells from FSGrid!"
         Emag = sqrt.(sum(readvariable(meta, "fg_e").^2, dims=1))
      end
      @inbounds for i = eachindex(ρ) # sparsity/inner boundary
         Emag[i] == 0.0 && (Emag[i] = NaN)
      end
      Emag
   end,
   :Vmag => function (meta, ids=UInt64[])
      rho_ = findfirst(endswith("rho"), meta.variable)
      ρ = isempty(ids) ?
         readvariable(meta, meta.variable[rho_]) :
         readvariable(meta, meta.variable[rho_], ids)
      Vmag = isempty(ids) ?
         vec(sqrt.(sum(readvariable(meta, "proton/vg_v").^2, dims=1))) :
         vec(sqrt.(sum(readvariable(meta, "proton/vg_v", ids).^2, dims=1)))
      @inbounds for i = eachindex(ρ) # sparsity/inner boundary
         Vmag[i] == 0.0 && (Vmag[i] = NaN)
      end
      Vmag
   end,
   :Rhom => function (meta, ids=UInt64[])
      if hasvariable(meta, "vg_rhom")
         ρm = isempty(ids) ?
            readvariable(meta, "vg_rhom") :
            readvariable(meta, "vg_rhom", ids)
      elseif hasvariable(meta, "proton/vg_rho")
         ρm = isempty(ids) ?
            readvariable(meta, "proton/vg_rho") .* mᵢ :
            readvariable(meta, "proton/vg_rho", ids) .* mᵢ
      end
      ρm
   end,
   :P => function (meta, ids=UInt64[]) # scalar pressure
      if hasvariable(meta, "vg_pressure")
         P = isempty(ids) ?
            readvariable(meta, "vg_pressure") :
            readvariable(meta, "vg_pressure", ids)
      else
         Pdiag = isempty(ids) ?
            readvariable(meta, "proton/vg_ptensor_diagonal") :
            readvariable(meta, "proton/vg_ptensor_diagonal", ids)
         P = vec(mean(Pdiag, dims=1))
      end
      P
   end,
   :VS => function (meta, ids=UInt64[]) # sound speed
      P = readvariable(meta, "P", ids)
      ρm = readvariable(meta, "Rhom", ids)
      @inbounds for i = eachindex(ρm) # sparsity/inner boundary
         ρm[i] == 0.0 && (ρm[i] = NaN)
      end
      vs = @. √( (P*5.0f0/3.0f0) / ρm )
   end,
   :VA => function (meta, ids=UInt64[]) # Alfvén speed
      ρm = readvariable(meta, "Rhom", ids)
      @inbounds for i = eachindex(ρm) # sparsity/inner boundary
         ρm[i] == 0.0 && (ρm[i] = NaN)
      end
      Bmag = readvariable(meta, "Bmag", ids)
      VA = @. $vec(Bmag) / √(ρm*μ₀)
   end,
   :MA => function (meta, ids=UInt64[]) # Alfvén Mach number
      V = readvariable(meta, "Vmag", ids)
      @inbounds for i = eachindex(V) # sparsity/inner boundary
         V[i] == 0.0 && (V[i] = NaN)
      end
      VA = readvariable(meta, "VA", ids)
      V ./ VA
   end,
   :MS => function (meta, ids=UInt64[]) # Sonic Mach number
      V = readvariable(meta, "Vmag", ids)
      @inbounds for i = eachindex(V) # sparsity/inner boundary
         V[i] == 0.0 && (V[i] = NaN)
      end
      VS = readvariable(meta, "VS", ids)
      V ./ VS
   end,
   :Vpar => function (meta, ids=UInt64[]) # velocity ∥ B
      V = isempty(ids) ?
         readvariable(meta, "proton/vg_v") :
         readvariable(meta, "proton/vg_v", ids)
      if isempty(ids)
         b = readvariable(meta, "vg_b_vol") ./ readvariable(meta, "Bmag")
      else
         b = readvariable(meta, "vg_b_vol", ids) ./ readvariable(meta, "Bmag", ids)
      end
      [V[:,i] ⋅ b[:,i] for i in 1:size(V,2)]
   end,
   :Vperp => function (meta, ids=UInt64[]) # velocity ⟂ B
      V = isempty(ids) ?
         readvariable(meta, "proton/vg_v") :
         readvariable(meta, "proton/vg_v", ids)
      if isempty(ids)
         b = readvariable(meta, "vg_b_vol") ./ readvariable(meta, "Bmag")
      else
         b = readvariable(meta, "vg_b_vol", ids) ./ readvariable(meta, "Bmag", ids)
      end
      Vperp = zeros(eltype(V), size(V, 2))
      # Avoid sqrt of negative values, but does not guarantee orthogonality.
      @inbounds for i in eachindex(Vperp)
         Vpar = V[:,i] ⋅ b[:,i] .* b[:,i]
         Vperp[i] = norm(V[:,i] - Vpar)
      end
      Vperp
   end,
   :T => function (meta, ids=UInt64[]) # scalar temperature
      P = readvariable(meta, "P", ids)
      n = isempty(ids) ?
         readvariable(meta, "proton/vg_rho") :
         readvariable(meta, "proton/vg_rho", ids)
      @inbounds for i = eachindex(n) # sparsity/inner boundary
         n[i] == 0.0 && (n[i] = NaN)
      end
      T = @. P / (n*kB)
   end,
   :Vth => function (meta, ids=UInt64[]) # thermal velocity
      T = readvariable(meta, "T", ids)
      @. √(3 * kB * T / mᵢ) # assume proton
   end,
   :Ppar => function (meta, ids=UInt64[]) # P component ∥ B
      P = readvariable(meta, "Protated", ids)
      @views P[3,3,:]
   end,
   :Pperp => function (meta, ids=UInt64[]) # P component ⟂ B
      P = readvariable(meta, "Protated", ids)
      Pperp = [0.5f0(P[1,1,i] + P[2,2,i]) for i in 1:size(P,3)]
   end,
   :Tpar => function (meta, ids=UInt64[]) # T component ∥ B
      P = readvariable(meta, "Protated", ids)
      n = isempty(ids) ?
         readvariable(meta, "proton/vg_rho") :
         readvariable(meta, "proton/vg_rho", ids)
      @inbounds for i = eachindex(n) # sparsity/inner boundary
         n[i] == 0.0 && (n[i] = NaN)
      end
      @. P[3,3,:] / (n*kB)
   end,
   :Tperp => function (meta, ids=UInt64[]) # scalar T component ⟂ B
      P = readvariable(meta, "Protated", ids)
      n = isempty(ids) ?
         readvariable(meta, "proton/vg_rho") :
         readvariable(meta, "proton/vg_rho", ids)
      @inbounds for i = eachindex(n) # sparsity/inner boundary
         n[i] == 0.0 && (n[i] = NaN)
      end
      Pperp = [0.5(P[1,1,i] + P[2,2,i]) for i in 1:size(P,3)]
      @. Pperp / (n*kB)
   end,
   :Tanisotropy => function (meta, ids=UInt64[]) # T⟂ / T∥
      Tperp = readvariable(meta, "Tperp", ids)
      Tpar = readvariable(meta, "Tpar", ids)
      @. Tperp / Tpar
   end,
   :J => function (meta, ids=UInt64[])
      @assert isempty(ids) "Do not support current calculation for selected cells!"
      B = readvariable(meta, "vg_b_vol")
      B = reshape(B, 3, meta.ncells...)
      J = curl(meta.dcoord, B) ./ μ₀
      J = reshape(J, 3, :) # To be consistent with shape assumptions
   end,
   :Protated => function (meta, ids=UInt64[])
      # Rotate the pressure tensor to align the 3rd direction with B
      B = isempty(ids) ?
         readvariable(meta, "vg_b_vol") :
         readvariable(meta, "vg_b_vol", ids)
      Pdiag = isempty(ids) ?
         readvariable(meta, "proton/vg_ptensor_diagonal") :
         readvariable(meta, "proton/vg_ptensor_diagonal", ids)
      Podiag = isempty(ids) ?
         readvariable(meta, "proton/vg_ptensor_offdiagonal") :
         readvariable(meta, "proton/vg_ptensor_offdiagonal", ids)
      P = zeros(Float32, 3, 3, size(Pdiag, 2))
      @inbounds for i = 1:size(P, 3)
         P[1,1,i] = Pdiag[1,i]
         P[2,2,i] = Pdiag[2,i]
         P[3,3,i] = Pdiag[3,i]
         P[1,2,i] = P[2,1,i] = Podiag[1,i]
         P[2,3,i] = P[3,2,i] = Podiag[2,i]
         P[3,1,i] = P[1,3,i] = Podiag[3,i]
         @views rotateWithB!(P[:,:,i], B[:,i])
      end
      P
   end,
   :Panisotropy => function (meta, ids=UInt64[]) # P⟂ / P∥
      PR = readvariable(meta, "Protated", ids)
      @. 0.5f0*(PR[1,1,:] + PR[2,2,:]) / PR[3,3,:]
   end,
   :Pdynamic => function (meta, ids=UInt64[])
      V = readvariable(meta, "Vmag", ids)
      ρm = readvariable(meta, "Rhom", ids)
      @. ρm * V * V
   end,
   :Pb => function (meta, ids=UInt64[])
      B² = isempty(ids) ?
         vec(sum(readvariable(meta, "vg_b_vol").^2, dims=1)) :
         vec(sum(readvariable(meta, "vg_b_vol", ids).^2, dims=1))
      @inbounds for i = eachindex(B²) # sparsity/inner boundary
         B²[i] == 0.0 && (B²[i] = NaN)
      end
      mu2inv = 0.5/μ₀
      @. B² * mu2inv
   end,
   :Poynting => function (meta, ids=UInt64[])
      if hasvariable(meta, "vg_b_vol") && hasvariable(meta, "vg_e_vol")
         E = isempty(ids) ?
            readvariable(meta, "vg_e_vol") :
            readvariable(meta, "vg_e_vol", ids)
         B = isempty(ids) ?
            readvariable(meta, "vg_b_vol") :
            readvariable(meta, "vg_b_vol", ids)
      elseif hasvariable(meta, "fg_e") && hasvariable(meta, "fg_b")
         @assert isempty(ids) "Do not support selecting cells in FSGrid!"
         E = readvariable(meta, "fg_e")
         B = readvariable(meta, "fg_b")
      else
         E = readvariable(meta, "E")
         B = readvariable(meta, "B")
      end
      F = similar(E)
      Rpost = CartesianIndices(size(E)[2:end])
      @inbounds for i in Rpost
         F[:,i] = E[:,i] × B[:,i] ./ μ₀
      end
      F
   end,
   :Agyrotropy => function (meta, ids=UInt64[])
      # non-gyrotropy measure Q [Swisdak 2016]
      # Original derivation for electrons. Here we do protons first.
      if isempty(ids)
         Pdiag = readvariable(meta, "proton/vg_ptensor_diagonal")
         Podiag = readvariable(meta, "proton/vg_ptensor_offdiagonal")
         B = readvariable(meta, "vg_b_vol")
      else
         Pdiag = readvariable(meta, "proton/vg_ptensor_diagonal", ids)
         Podiag = readvariable(meta, "proton/vg_ptensor_offdiagonal", ids)
         B = readvariable(meta, "vg_b_vol", ids)
      end
      Pxx = selectdim(Pdiag, 1, 1)
      Pyy = selectdim(Pdiag, 1, 2)
      Pzz = selectdim(Pdiag, 1, 3)
      Pxy = selectdim(Podiag, 1, 1)
      Pyz = selectdim(Podiag, 1, 2) # Warning: the order may be wrong!
      Pxz = selectdim(Podiag, 1 ,3) # Warning: the order may be wrong!

      Bx = selectdim(B, 1, 1)
      By = selectdim(B, 1, 2)
      Bz = selectdim(B, 1, 3)

      I₁ = @. Pxx + Pyy + Pzz
      I₂ = @. Pxx*Pyy + Pxx*Pzz + Pyy*Pzz - Pxy*Pxy - Pyz*Pyz - Pxz*Pxz

      Ppar = @. (Bx*Bx*Pxx + By*By*Pyy + Bz*Bz*Pzz +
	      2*(Bx*By*Pxy + Bx*Bz*Pxz + By*Bz*Pyz))/B²
      Qsqr = @. √(1 - 4I₂/((I₁ - Ppar)*(I₁ + 3Ppar)))
   end,
   :Beta => function (meta, ids=UInt64[])
      P = readvariable(meta, "P", ids)
      B² = isempty(ids) ?
         vec(sum(readvariable(meta, "vg_b_vol").^2, dims=1)) :
         vec(sum(readvariable(meta, "vg_b_vol", ids).^2, dims=1))
      @inbounds for i = eachindex(B²) # sparsity/inner boundary
         B²[i] == 0.0 && (B²[i] = NaN)
      end
      @. 2 * μ₀ * P / B²
   end,
   :IonInertial => function (meta, ids=UInt64[])
      n = isempty(ids) ?
         readvariable(meta, "proton/vg_rho") :
         readvariable(meta, "proton/vg_rho", ids)
      Z = 1
      ωi = @. √(n/(mᵢ*μ₀)) * Z * qᵢ
      di = @. c / ωi
   end,
   :Larmor => function (meta, ids=UInt64[])
      Vth = readvariable(meta, "Vth", ids)
      B = readvariable(meta, "Bmag", ids)
      rg = @. mᵢ * Vth / (qᵢ * B)
   end,
   :Gyroperiod => function (meta, ids=UInt64[])
      B = readvariable(meta, "Bmag", ids)
      T = @. 2π * mᵢ / (qᵢ * B)
   end,
   :Gyrofrequency => function (meta, ids=UInt64[]) # [1/s]
      B = readvariable(meta, "Bmag", ids)
      f = @. qᵢ * B / (mᵢ * 2π)
   end,
   :Plasmaperiod => function (meta, ids=UInt64[])
      n = isempty(ids) ?
         readvariable(meta, "proton/vg_rho") :
         readvariable(meta, "proton/vg_rho", ids)
      T = @. 2π / (qᵢ * √(n  / (mᵢ * ϵ₀)))
   end,
   :Omegap => function (meta, ids=UInt64[]) # plasma frequency, [1/s]
      n = isempty(ids) ?
         readvariable(meta, "proton/vg_rho") :
         readvariable(meta, "proton/vg_rho", ids)
      ωₚ = @. qᵢ * √(n  / (mᵢ * ϵ₀)) / 2π
   end,
)


text:
"""
These are the simplified version of communication methods, only using MPI.COMM_WORLD.
"""

function Barrier()
    MPI.Barrier(COMM_WORLD)
end

function Bcast!(arr::AbstractArray; root::Integer=0)
    MPI.Bcast!(arr, root, COMM_WORLD)
end

function Send(arr::AbstractArray, dest::Integer; tag::Integer=0)
    MPI.Send(arr, dest, tag, COMM_WORLD)
end

function Recv!(arr::AbstractArray, src::Integer; tag::Integer=0)
    MPI.Recv!(arr, src, tag, COMM_WORLD)
end

function Isend(arr::AbstractArray, dest::Integer; tag::Integer=0)
    MPI.Isend(arr, dest, tag, COMM_WORLD)
end

function Irecv!(arr::AbstractArray, src::Integer; tag::Integer=0)
    MPI.Irecv!(arr, src, tag, COMM_WORLD)
end

function Reduce!(sendarr::AbstractArray, recvarr::AbstractArray; op=MPI.SUM, root::Integer=0)
    MPI.Reduce!(sendarr, recvarr, op, root, COMM_WORLD)
end

function Allreduce!(arr::AbstractArray; op=MPI.SUM)
    MPI.Allreduce!(arr, op, COMM_WORLD)
end

function Allgather!(sendarr::AbstractArray, recvarr::AbstractArray)
    @assert length(recvarr) >= Size() * length(sendarr)
    MPI.Allgather!(sendarr, recvarr, length(sendarr), COMM_WORLD)
end

function Allgather!(sendarr::AbstractArray, recvarr::AbstractArray, count::Integer)
    @assert length(recvarr) >= Size() * count
    MPI.Allgather!(sendarr, recvarr, count, COMM_WORLD)
end

function Allgatherv!(sendarr::AbstractArray, recvarr::AbstractArray, counts::Vector{<:Integer})
#### TODO: auto-detect counts (is it doable?)
    MPI.Allgatherv!(sendarr, recvarr, convert(Vector{Cint}, counts), COMM_WORLD)
end

function Scatter!(sendarr::AbstractArray, recvarr::AbstractArray; root::Integer=0)
    MPI.Scatter!(sendarr, recvarr, length(recvarr), root, COMM_WORLD)
end

function Scatterv!(sendarr::AbstractArray, recvarr::AbstractArray, counts::Vector{<:Integer}; root::Integer=0)
#### TODO: auto-detect counts (is it doable?)
    MPI.Scatterv!(sendarr, recvarr, convert(Vector{Cint}, counts), root, COMM_WORLD)
end

function Gather!(sendarr::AbstractArray, recvarr::AbstractArray; root::Integer=0)
    MPI.Gather!(sendarr, recvarr, length(sendbuf), root, COMM_WORLD)
end

function Gatherv!(sendarr::AbstractArray, recvarr::AbstractArray, counts::Vector{<:Integer}; root::Integer=0)
#### TODO: auto-detect counts (is it doable?)
    MPI.Gatherv!(sendarr, recvarr, convert(Vector{Cint}, counts), root, COMM_WORLD)
end


text:
@reexport module CliffWalking

export CliffWalkingEnv

using ReinforcementLearningBase

const NX = 4
const NY = 12
const Start = CartesianIndex(4, 1)
const Goal = CartesianIndex(4, 12)
const Actions = [
    CartesianIndex(0, -1),  # left
    CartesianIndex(0, 1),   # right
    CartesianIndex(-1, 0),  # up
    CartesianIndex(1, 0),   # down
]

const LinearInds = LinearIndices((NX, NY))

function iscliff(p::CartesianIndex{2})
    x, y = Tuple(p)
    x == 4 && y > 1 && y < NY
end

Base.@kwdef mutable struct CliffWalkingEnv <: AbstractEnv
    position::CartesianIndex{2} = Start
end

RLBase.get_observation_space(env::CliffWalkingEnv) = DiscreteSpace(length(LinearInds))
RLBase.get_action_space(env::CliffWalkingEnv) = DiscreteSpace(length(Actions))

function (env::CliffWalkingEnv)(a::Int)
    x, y = Tuple(env.position + Actions[a])
    env.position = CartesianIndex(min(max(x, 1), NX), min(max(y, 1), NY))
    nothing
end

RLBase.observe(env::CliffWalkingEnv) =
    (
        reward = env.position == Goal ? 0.0 : (iscliff(env.position) ? -100.0 : -1.0),
        terminal = env.position == Goal || iscliff(env.position),
        state = LinearInds[env.position],
    )

function RLBase.reset!(env::CliffWalkingEnv)
    env.position = Start
    nothing
end

end


text:
#!/usr/bin/env julia

"""
Utilities for producing and operating on numerical sequences.
"""
module Sequences

"""
    Fibonacci([T::Type{<:Integer}])

Iterable object that returns subsequent terms of the Fibonacci sequence.

## Fields
- `T`: The integer type of each term in the sequence. Defaults to `Int`

## Examples
```julia-repl
julia> first(Fibonacci())
1

julia> first(Fibonacci(UInt8))
0x01

julia> for n in Fibonacci()
           println(n)
           if n > 10
               break
           end
       end
1
2
3
5
8
13

julia> for (i, n) in Iterators.enumerate(Fibonacci())
           println(n)
           if i == 5
               break
           end
       end
1
2
3
5
8
```

See also: [`fibonacci`](@ref)
"""
struct Fibonacci
    T::Type{<:Integer}
end
Fibonacci() = Fibonacci(Int)

function Base.iterate(
    iter::Fibonacci,
    state::Tuple{Integer, Integer} = (zero(iter.T), oneunit(iter.T))
)
    fib_old, fib_new = state
    fib_old, fib_new = fib_new, (fib_old + fib_new)
    return fib_new, (fib_old, fib_new)
end

"""
    arithmetic_series(a::Integer, n::Integer, d::Integer = 1)

Returns the sum of the arithmetic sequence with parameters `a`, `n`, `d`.

# Arguments
- `a`: The first term in the sequence
- `n`: The total number of terms in the sequence
- `d`: The difference between two consecutive terms in the sequence

# Preconditions
- `n ≥ 0`

# Examples
```jldoctest
julia> arithmetic_series(1, 4)
10

julia> arithmetic_series(2, 6, -1)
-3

julia> arithmetic_series(3, 5, 4)
55
```
"""
function arithmetic_series(a::Integer, n::Integer, d::Integer = 1)
    return n * (2a + (n - 1)d) ÷ 2
end

"""
    fibonacci([T::Type{<:Integer}], n::Integer)

Computes the specified term `F(n)` of the Fibonacci sequence, defined by the
relation `F(m) = F(m - 1) + F(m - 2)` with `F(0) = F(1) = 1`.

## Arguments
- `T`: The integer type of each term in the sequence. Defaults to `typeof(n)`
- `n`: Specifies the term `F(n)` in the Fibonacci sequence to return

## Preconditions
- `n ≥ 0`

## Examples
```jldoctest
julia> fibonacci(0)
1

julia> fibonacci(4)
5

julia> fibonacci(UInt8, 5)
0x08
```

See also: [`Fibonacci`](@ref)
"""
function fibonacci(n::Integer)
    return fibonacci(typeof(n), n)
end
function fibonacci(T::Type{<:Integer}, n::Integer)
    # Fibonacci iterates from F(1), so handle F(0) ourselves
    if n == 0
        return oneunit(T)
    end

    for (i, fib_num) in enumerate(Fibonacci(T))
        if i == n
            return fib_num
        end
    end

    # Fibonacci iterator should never terminate
    throw(ErrorException("iteration terminated unexpectedly"))
end

"""
    next_multiple(n::T, min_val::T) where T <: Integer

TODO: Function description goes here

## Arguments
- TODO: Arguments go here

## Preconditions
- TODO: Preconditions go here

## Examples
```jldoctest
julia> next_multiple()

```
"""
function next_multiple(n::T, min_val::T) where T <: Integer
    return min_val + (n - min_val % n) % n
end

"""
    sum_of_squares(n::Integer)

Computes the sum of the squares of the first `n` natural numbers.

## Preconditions
- `n ≥ 0`

## Examples
```jldoctest
julia> sum_of_squares(1)
1

julia> sum_of_squares(0x02)
0x05

julia> sum_of_squares(4)
30
```
"""
function sum_of_squares(n::Integer)
    T = typeof(n)
    cube_term = convert(T, 2) * n^3
    square_term = convert(T, 3) * n^2
    return (cube_term + square_term + n) ÷ convert(T, 6)
end

"""
    triangular(n::Integer)

Computes the `n`th triangular number `T(n)`, with `T(1) = 1`. The resulting
value is also the sum of the first `n` natural numbers.

## Preconditions
- `n ≥ 0`

## Examples
```jldoctest
julia> triangular(1)
1

julia> triangular(0x03)
0x06

julia> triangular(4)
10
```
"""
triangular(n::Integer) = n * (n + oneunit(n)) ÷ convert(typeof(n), 2)

end

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
module GMeans

using Query
using ..Imports

export call

function call(values)
  gmeans_instance = gmeans.gmeans(values).process()
  clusters = gmeans_instance.get_clusters()
  centers = gmeans_instance.get_centers()

  cluster_results = size(clusters, 1) <= 1 ? mapslices(x -> [x], clusters, dims=2)[:] : clusters

  Dict(:clusters => cluster_results |> @map(x -> x |> @map(values[_+1]) |> collect) |> collect,
       :centers => mapslices(x -> [x], centers, dims=2)[:])
end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Materials, FEMMaterials, Test

m1 = 1.0e-3*[-0.3, -0.3, 1.0, 0.0, 0.0, 0.0]
m2 = 1.0e-3*[-0.5, -0.5, 1.0, 0.0, 0.0, 0.0]
# mat = Material(IdealPlastic, ())
# mat.properties.youngs_modulus = 200.0e3
# mat.properties.poissons_ratio = 0.3
# mat.properties.yield_stress = 100.0
# mat.stress = zeros(6)

parameters = IdealPlasticParameterState(youngs_modulus = 200.0e3,
                                        poissons_ratio = 0.3,
                                        yield_stress = 100.0)

mat = IdealPlastic(parameters=parameters)
times = [0.0]
strains = [zeros(6)]
dt = 0.5
# Go to elastic border
push!(times, times[end]+dt)
push!(strains, strains[end] + m1*dt)

# Proceed to plastic flow
push!(times, times[end]+dt)
push!(strains, strains[end] + m2*dt)

# Reverse direction
push!(times, times[end]+dt)
push!(strains, strains[end] - m1*dt)

# Continue and pass yield criterion
push!(times, times[end]+dt)
push!(strains, strains[end] - (m1 + m2)*dt)

sim = Simulator(mat)
FEMMaterials.initialize!(sim, strains, times)
FEMMaterials.run!(sim)

s33s = [s[3] for s in sim.stresses]
@test isapprox(s33s, [0.0, 100.0, 100.0, 0.0, -100.0])


text:
module CoordinateSystems

export CoordinateSystem, from_global, to_global, rotate_to_global

using ImmutableArrays

immutable CoordinateSystem
    center::Vector3{Float64}
    rotation_matrix::Matrix3x3
end

normalize(a::Vector3{Float64}) = a/norm(a)

function CoordinateSystem(pc::Vector3{Float64},pz::Vector3{Float64},pxz::Vector3{Float64})
    vecz = normalize(pz - pc)
    vecy = normalize(cross(vecz,pxz-pc))
    vecx = normalize(cross(vecy,vecz))
    CoordinateSystem(pc,Matrix3x3(vecx,vecy,vecz))
end

from_global(csys::CoordinateSystem,p::Vector3{Float64}) = inv(csys.rotation_matrix) * (p - csys.center)
rotate_to_global(csys::CoordinateSystem,p::Vector3{Float64}) = csys.rotation_matrix * p
to_global(csys::CoordinateSystem,p::Vector3{Float64}) = rotate_to_global(csys,p) + csys.center

end # module


text:
using FractionalCalculus
using Test

@testset "Test Caputo Direct Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 0, 1, 1e-8, CaputoDirect())[1], 2/sqrt(pi); atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 0, 3.2, 1e-8, CaputoDirect())[1], 4.300306216488329e2; atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.25, 0, 3.2, 1e-8, CaputoDirect())[1], 382.1228113951680; atol = 1e-5)
    @test isapprox(fracdiff(x->(x-3)^3+(x+2)^2+5, 0.5, 0, 0.83, 1e-8, CaputoDirect())[1], 23.899930581118699; atol = 1e-5)
    @test isapprox(fracdiff(x->(x-3)^3+(x+2)^2+5, 0.25, 0, 0.83, 1e-8, CaputoDirect())[1], 22.96354626572943; atol = 1e-5)
    @test isapprox(fracdiff(x->cos(x), 0.5, 0, 0.34, 1e-8, CaputoDirect())[1], -0.147174774607683; atol = 1e-5)
    @test isapprox(fracdiff(x->cos(x), 0.25, 0, 0.34, 1e-8, CaputoDirect())[1], -0.093074344902431; atol = 1e-5)
    @test isapprox(fracdiff(x->exp(-x^2), 0.5, 0, 2.3, 1e-8, CaputoDirect())[1], -0.505891017154289; atol = 1e-5)
    @test isapprox(fracdiff(x->exp(-x^2), 0.25, 0, 2.3, 1e-8, CaputoDirect())[1], -0.762927553656252; atol = 1e-5)

    @test isapprox(fracdiff(x->x, 0.5, 0, [1, 2, 3], 0.0001, CaputoDirect()), [2/sqrt(pi), 1.5957691216057408, 1.9544100476116828]; atol=1e-3)
end

@testset "Test Caputo Piecewise Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 1e-4, CaputoPiecewise()), 2/sqrt(pi); atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 3.2, 1e-4, CaputoPiecewise()), 4.300306216488329e2; atol = 1e-5)
    @test isapprox(fracdiff(cos, 0.5, 0.34, 1e-4, CaputoPiecewise()), -0.147174774607683; atol = 1e-5)
    @test isapprox(fracdiff(x->exp(-x^2), 0.25, 2.3, 1e-3, CaputoPiecewise()), -0.762927553656252; atol = 1e-5)
end

@testset "Test Caputo Piecewise Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, [1, 2, 3], 1e-4, CaputoPiecewise()), [1.1283791670557781, 1.5957691213491374, 1.9544100474864603])
end

@testset "Test Caputo Diethelm Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 1e-4, CaputoDiethelm()), 2/sqrt(pi); atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 3.2, 1e-5, CaputoDiethelm()), 4.300306216488329e2; atol = 1e-5)
end

@testset "Test Caputo Diethelm Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, [1, 2, 3], 0.0001, CaputoDiethelm()), [1.1283791670955023, 1.5957691216057386, 1.954410047611745])
end

@testset "Test Caputo High Order Fractional Derivative" begin
    @test isapprox(fracdiff(x->x^8, 0.5, 1, 0.05, 3, CaputoHighOrder()), 2.8542562386883445; atol=1e-4)
end

@testset "Test Grunwald-Letnikov Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 0, 0.5, GLDirect())[1], 0.7978845583186518; atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 0, 3.2, GLDirect())[1], 4.300306216488329e2; atol = 1e-5)
    @test isapprox(fracdiff(x->x, 0.5, collect(0:0.01:1), 2, GLHighPrecision())[end], 2/sqrt(pi); atol = 1e-4)
end

@testset "Test GL_Multiplicative_Additive()" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.009, GLMultiplicativeAdditive()), 2/sqrt(pi); atol=1e-2)
end

@testset "Test GL_Lagrange_Three_Point_Interp" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.009, GLLagrangeThreePointInterp()), 2/sqrt(pi); atol=1e-2)
end

@testset "Test GL_Finite_Difference" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.009, GLFiniteDifference()), 2/sqrt(pi); atol=1e-2)
end

@testset "Test Grunwald Letnikov Lagrange_Three_Point_Interp" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.01, GLLagrangeThreePointInterp()), 2/sqrt(pi); atol=1e-4)
end

@testset "Test RLDiff Matrix Fractional Deriavtive" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.001, RLDiffMatrix())[end], 2/sqrt(pi); atol = 1e-3)
end

@testset "Test RLDiff Approx" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.001, RLDiffApprox()), 2/sqrt(pi); atol = 1e-4)

    # Test Array input end point
    result = fracdiff(x->x, 0.5, [1, 2, 3], 0.001, RLDiffApprox())
    @test isapprox(result[1], 2/sqrt(pi); atol=1e-4)
    @test isapprox(result[2], 1.5957691216057306; atol=1e-4)
    @test isapprox(result[3], 1.954410047611678; atol=1e-4)
end

@testset "Test RL Linear Spline Interpolation" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.0001, RLLinearSplineInterp()), 2/sqrt(pi); atol = 1e-4)

    # Test Array input end point
    result = fracdiff(x->x, 0.5, [1, 2, 3], 0.001, RLLinearSplineInterp())
    @test isapprox(result[1], 2/sqrt(pi); atol=1e-4)
    @test isapprox(result[2], 1.5957691216057306; atol=1e-4)
    @test isapprox(result[3], 1.954410047611678; atol=1e-4)
end

@testset "Test RL_G1 method" begin
    @test isapprox(fracdiff(x->x, 0.5, 0, 1, 0.006, RLG1()), 2/sqrt(pi); atol=1e-2)
end

@testset "Test RL_D method" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.001, RLD()), 2/sqrt(pi); atol=1e-4)
end

@testset "Test Hadamard Fractional Derivative" begin
    @test isapprox(fracdiff(log, 0.5, 1, 2, 1/4000000, HadamardLRect()), 0.9391460; atol=1e-4)
    @test isapprox(fracdiff(log, 0.3, 1, 2, 1/4000000, HadamardLRect()), 0.8514935; atol=1e-4)
    @test isapprox(fracdiff(log, 0.5, 1, 2, 1/4000000, HadamardRRect()), 0.9391460; atol=1e-4)
    @test isapprox(fracdiff(log, 0.3, 1, 2, 1/4000000, HadamardRRect()), 0.8514936; atol=1e-4)
    @test isapprox(fracdiff(log, 0.5, 1, 2, 1/4000000, HadamardTrap()), 0.9391460; atol=1e-4)
    @test isapprox(fracdiff(log, 0.3, 1, 2, 1/4000000, HadamardTrap()), 0.8514935; atol=1e-4)
end

@testset "Test Riesz Symmetric fractional derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.5, RieszSymmetric()), [0.2651650429449553; 0.2651650429449553; -0.39774756441743303]; atol=1e-2)
end

@testset "Test macros" begin
    @test isapprox(@fracdiff(x->x, 0.5, 1), 1.1283791670955126; atol=1e-4)
    @test isapprox(@fracdiff(x->x^5, 0.5, 3.2), 4.300306216488329e2; atol=1e-2)
end

@testset "Test auxiliary functions" begin
    @test isapprox(RieszMatrix(0.5, 2, 0.1), [-1.58114 1.3835; 1.3835 -1.58114]; atol=1e-3)
    @test isapprox(omega(3, 2), [1, -2, 1, 0]; atol=1e-3)
    @test isapprox(B(2, 2), [1.0 0.0; -2.0 1.0]; atol=1e-3)
    @test isapprox(genfun(2), [1.5 -2 0.5]; atol=1e-3)
    @test isapprox(first_order(exp, 2, 1e-6), exp(2); atol=1e-5)
    @test isapprox(first_order(sin, 2, 1e-6), cos(2); atol=1e-5)
end


text:
using SparseArrays
using LinearAlgebraicRepresentation
LAR = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

V1, (VV,EV,FV,CV) = LAR.simplex(3, true)
tetra1 = V1, EV, FV, CV
V2 = [
  0.0 0.5 0.1 0.1
  0.0 0.1 0.5 0.1
  0.0 0.1 0.1 0.5
]
tetra2 = V2, EV, FV, CV
twotetra = LAR.Struct([ tetra, tetra2 ])
V,EV,FV,CV = LAR.struct2lar(twotetra)
GL.VIEW([ GL.GLGrid(V,FV, GL.Point4d(1,1,1,0.2)) ]);
m = Lar.Model(V);
Lar.addModelCells!(m, 1, convert(LAR.ChainOp, LAR.coboundary_0(EV::LAR.Cells)))
Lar.addModelCells!(m, 2, LAR.coboundary_1(V, FV::LAR.Cells, EV::LAR.Cells))
Lar.viewModel(m,1)

sp_idx = LarA.spaceIndex(m, 2)

de_models = [
	LarA.face_decomposition(m, face_idx, sp_idx[face_idx])
	for face_idx = 1 : size(m, 2, 1)
]

mfrag = Lar.mergeMultipleModels(de_models)

triangulated_faces = LAR.triangulate(convert(Lar.Points, mfrag.G'), mfrag.T[1:2])
FVs = convert(Array{Lar.Cells}, triangulated_faces);
GL.VIEW( GL.GLExplode(mfrag.G,FVs,expl,expl,expl,99,1) );

mArranged = Lar.deepcopy(mfrag)
Lar.addModelCells!(mArranged, 3, LAR.Arrangement.minimal_3cycles(
		convert(Lar.Points, mfrag.G'), mfrag.T[1], mfrag.T[2])
)

V,CVs,FVs,EVs = LAR.pols2tria(mArranged.G, mArranged.T[1], mArranged.T[2], mArranged.T[3])

GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
meshes = GL.GLExplode(V,CVs[1:end],8,4,6,99,0.5);
push!( meshes, GL.GLFrame)
GL.VIEW( meshes );

#=
cop_EV = LAR.coboundary_0(EV::LAR.Cells);
cop_EW = convert(LAR.ChainOp, cop_EV);
cop_FE = LAR.coboundary_1(V, FV::LAR.Cells, EV::LAR.Cells);
W = convert(LAR.Points, V');

V, copEV, copFE, copCF = LAR.Arrangement.spatial_arrangement( W::LAR.Points, cop_EW::LAR.ChainOp, cop_FE::LAR.ChainOp)

EV = LAR.cop2lar(copEV)
FE = [findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
FV = [collect(Set(cat(EV[e] for e in FE[f]))) for f=1:length(FE)]
FV = convert(LAR.Cells, FV)
W = convert(LAR.Points, V')
WW = [[k] for k=1:size(W,2)]

GL.VIEW(GL.numbering(.5)((W,[WW,EV]) ));

triangulated_faces = LAR.triangulate(V, [copEV, copFE])
FVs = convert(Array{LAR.Cells}, triangulated_faces)
V = convert(LAR.Points, V')
GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99));

EVs = LAR.FV2EVs(copEV, copFE) # polygonal face fragments
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));


text:
module ReadJson

using JSON

function get_data_from_json(filename)
    return JSON.parsefile(filename)
end

function set_pallet_length(json)
    return get(json, "pallet_width", 0)
end

function set_pallet_width(json)
    return get(json, "pallet_depth", 0)
end

function set_boxes(json)
    dict_boxes = get(json, "items_to_stack", 0)
    list = [[get(dict_boxes[i], "width", 0), get(dict_boxes[i], "depth", 0), get(dict_boxes[i], "height", 0)] for i in 1:length(dict_boxes)]
    push!(list, [set_pallet_length(json), set_pallet_width(json), 0])
    return list
end

function set_masses(json)
    dict_boxes = get(json, "items_to_stack", 0)
    list = [get(dict_boxes[i], "weight", 0) for i in 1:length(dict_boxes)]
    push!(list, 0)
    return list
end

end


text:
# NOTE: These tests are not part of the active test suite, because they require Distributions.jl.
# DO NOT INCORPORATE INTO runtests.jl.

using Distributions
using LightGraphs
using StatsBase
using Test
import Random

import Base: -
import LightGraphs: randbn
import StatsBase: SummaryStats

function -(s::SummaryStats, t::SummaryStats)
    return SummaryStats(s.mean - t.mean,
        s.min - t.min,
        s.q25 - t.q25,
        s.median - t.median,
        s.q75 - t.q75,
        s.max - t.max)
end
function binomial_test(n, p, s)
    drand = rand(Binomial(n, p), s)
    lrand = Int64[randbn(n, p) for i in 1:s]

    ds = @show summarystats(drand)
    ls = @show summarystats(lrand)
    dσ = @show std(drand)
    lσ = @show std(lrand)

    summarydiff = @show ds - ls
    @test abs(summarydiff.mean) / ds.mean < .10
    @test abs(summarydiff.median) / ds.median < .10
    @test abs(summarydiff.q25) / ds.q25 < .10
    @test abs(summarydiff.q75) / ds.q75 < .10

    @show dσ - lσ
    @test abs(dσ - lσ) / dσ < .10
end
seed!(1234)
n = 10000
p = 0.3
s = 100000

for (n, p, s) in [(100, 0.3, 1000), (1000, 0.8, 1000), (10000, 0.25, 1000)]
    binomial_test(n, p, s)
end


text:
# This file is part of Kpax3. License is MIT.

# A1 = [1 0 0 0 0 0;
#       1 0 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 0 0 1]
#
# A2 = [1 0 0 0 0 0;
#       0 1 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 1 0 0;
#       0 0 0 0 1 0;
#       0 0 0 0 1 0]
#
# Each cutpoint greater than 1 has probability 1 / 5 of being selected
#
# COLUMN 2
# If cutpoint == 7, Kpax3.crossover breaks the second child:
#
# B1 = [1 0 0 0 0 0;
#       0 1 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 1 0 0;
#       0 0 0 0 1 0;
#       0 0 0 0 1 0]
#
# B2 = [1 0 0 0 0 0;
#       0 0 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 0 0 1]
#
# Pr(B1 = [1; 2; 3; 4; 5; 5] | cutpoint = 7) = 1
#
# Pr(B2 = [1; 1; 3; 3; 3; 6] | cutpoint = 7) = 1 / 2
# Pr(B2 = [1; 2; 3; 3; 3; 6] | cutpoint = 7) = 1 / 2
#
# COLUMN 3
# If cutpoint == 14, Kpax3.crossover does not break anything
#
# B1 = [1 0 0 0 0 0;
#       1 0 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 1 0 0;
#       0 0 0 0 1 0;
#       0 0 0 0 1 0]
#
# B2 = [1 0 0 0 0 0;
#       0 1 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 0 0 1]
#
# Pr(B1 = [1; 1; 3; 4; 5; 5] | cutpoint = 14) = 1
#
# Pr(B2 = [1; 2; 3; 3; 3; 6] | cutpoint = 14) = 1
#
# COLUMN 4
# If cutpoint == 21, Kpax3.crossover breaks the second child:
#
# B1 = [1 0 0 0 0 0;
#       1 0 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 1 0 0;
#       0 0 0 0 1 0;
#       0 0 0 0 1 0]
#
# B2 = [1 0 0 0 0 0;
#       0 1 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 0 0 0;
#       0 0 0 0 0 0;
#       0 0 0 0 0 1]
#
# Pr(B1 = [1; 1; 3; 4; 5; 5] | cutpoint = 21) = 1
#
# Pr(B2 = [1; 2; 3; 1; 1; 6] | cutpoint = 21) = (1 / 4) * (2 / 5) = 1 / 10
# Pr(B2 = [1; 2; 3; 2; 1; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 3; 1; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 4; 1; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
#
# Pr(B2 = [1; 2; 3; 1; 2; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 2; 2; 6] | cutpoint = 21) = (1 / 4) * (2 / 5) = 1 / 10
# Pr(B2 = [1; 2; 3; 3; 2; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 4; 2; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
#
# Pr(B2 = [1; 2; 3; 1; 3; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 2; 3; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 3; 3; 6] | cutpoint = 21) = (1 / 4) * (2 / 5) = 1 / 10
# Pr(B2 = [1; 2; 3; 4; 3; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
#
# Pr(B2 = [1; 2; 3; 4; 4; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
#
# Pr(B2 = [1; 2; 3; 1; 5; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 2; 5; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 3; 5; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 4; 5; 6] | cutpoint = 21) = (1 / 4) * (1 / 5) = 1 / 20
#
# COLUMN 5
# If cutpoint == 28, Kpax3.crossover breaks the second child:
#
# B1 = [1 0 0 0 0 0;
#       1 0 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 0 1 0;
#       0 0 0 0 1 0]
#
# B2 = [1 0 0 0 0 0;
#       0 1 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 1 0 0;
#       0 0 0 0 0 0;
#       0 0 0 0 0 1]
#
# Pr(B1 = [1; 1; 3; 3; 5; 5] | cutpoint = 28) = 1
#
# Pr(B2 = [1; 2; 3; 4; 1; 6] | cutpoint = 28) = 1 / 5
# Pr(B2 = [1; 2; 3; 4; 2; 6] | cutpoint = 28) = 1 / 5
# Pr(B2 = [1; 2; 3; 4; 3; 6] | cutpoint = 28) = 1 / 5
# Pr(B2 = [1; 2; 3; 4; 4; 6] | cutpoint = 28) = 1 / 5
# Pr(B2 = [1; 2; 3; 4; 5; 6] | cutpoint = 28) = 1 / 5
#
# COLUMN 6
# If cutpoint == 35, Kpax3.crossover breaks the first child:
#
# A1 = [1 0 0 0 0 0;
#       1 0 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 0 0 0]
#
# A2 = [1 0 0 0 0 0;
#       0 1 0 0 0 0;
#       0 0 1 0 0 0;
#       0 0 0 1 0 0;
#       0 0 0 0 1 0;
#       0 0 0 0 0 1]
#
# Pr(B1 = [1; 1; 3; 3; 3; 1] | cutpoint = 35) = 2 / 6
# Pr(B1 = [1; 1; 3; 3; 3; 3] | cutpoint = 35) = 3 / 6
# Pr(B1 = [1; 1; 3; 3; 3; 6] | cutpoint = 35) = 1 / 6
#
# Pr(B2 = [1; 2; 3; 4; 5; 6] | cutpoint = 35) = 1
#
# SUMMARY
# Unique partitions that are possible to be obtained from this Kpax3.crossover:
#
# Pr(B1 = [1; 2; 3; 4; 5; 5]) = Pr(B1 = R2) = 1 / 5
# Pr(B1 = [1; 1; 3; 4; 5; 5]) = 2 / 5
# Pr(B1 = [1; 1; 3; 3; 5; 5]) = 1 / 5
# Pr(B1 = [1; 1; 3; 3; 3; 1]) = (1 / 5) * (2 / 6) = 1 / 15
# Pr(B1 = [1; 1; 3; 3; 3; 3]) = (1 / 5) * (3 / 6) = 1 / 10
# Pr(B1 = [1; 1; 3; 3; 3; 6]) = (1 / 5) * (1 / 6) = 1 / 30
#
# Pr(B2 = [1; 1; 3; 3; 3; 6]) = (1 / 5) * (1 / 2) = 1 / 10
# Pr(B2 = [1; 2; 3; 3; 3; 6]) = (1 / 5) * (1 / 2) + (1 / 5) * 1 + (1 / 5) * (1 / 10) = 8 / 25
# Pr(B2 = [1; 2; 3; 1; 1; 6]) = (1 / 5) * (1 / 10) = 1 / 50
# Pr(B2 = [1; 2; 3; 2; 1; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 3; 1; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 4; 1; 6]) = (1 / 5) * (1 / 20) + (1 / 5) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 1; 2; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 2; 2; 6]) = (1 / 5) * (1 / 10) = 1 / 50
# Pr(B2 = [1; 2; 3; 3; 2; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 4; 2; 6]) = (1 / 5) * (1 / 20) + (1 / 5) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 1; 3; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 2; 3; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 4; 3; 6]) = (1 / 5) * (1 / 20) + (1 / 5) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 4; 4; 6]) = (1 / 5) * (1 / 20) + (1 / 5) * (1 / 5) = 1 / 20
# Pr(B2 = [1; 2; 3; 1; 5; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 2; 5; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 3; 5; 6]) = (1 / 5) * (1 / 20) = 1 / 100
# Pr(B2 = [1; 2; 3; 4; 5; 6]) = (1 / 5) * (1 / 20) + (1 / 5) * (1 / 5) + (1 / 5) * 1 = 1 / 4

function test_crossover()
  support = Kpax3.GASupport(10, 6)

  R1 = [1; 1; 2; 2; 2; 3]
  R2 = [1; 2; 3; 4; 5; 5]

  N = 1000000

  pr = [1 / 5; 2 / 5; 1 / 5; 1 / 15; 1 / 10; 1 / 30; 1 / 10; 8 / 25; 1 / 50; 1 / 100; 1 / 100; 1 / 20; 1 / 100; 1 / 50; 1 / 100; 1 / 20; 1 / 100; 1 / 100; 1 / 20; 1 / 20; 1 / 100; 1 / 100; 1 / 100; 1 / 4]

  tmp = zeros(Float64, length(pr))

  for i in 1:N
    Kpax3.crossover!(R1, R2, support)

    if support.oi.R == [1; 2; 3; 4; 5; 5]
      tmp[ 1] += 1
    elseif support.oi.R == [1; 1; 3; 4; 5; 5]
      tmp[ 2] += 1
    elseif support.oi.R == [1; 1; 3; 3; 5; 5]
      tmp[ 3] += 1
    elseif support.oi.R == [1; 1; 3; 3; 3; 1]
      tmp[ 4] += 1
    elseif support.oi.R == [1; 1; 3; 3; 3; 3]
      tmp[ 5] += 1
    elseif support.oi.R == [1; 1; 3; 3; 3; 6]
      tmp[ 6] += 1
    end

    if support.oj.R == [1; 1; 3; 3; 3; 6]
      tmp[ 7] += 1
    elseif support.oj.R == [1; 2; 3; 3; 3; 6]
      tmp[ 8] += 1
    elseif support.oj.R == [1; 2; 3; 1; 1; 6]
      tmp[ 9] += 1
    elseif support.oj.R == [1; 2; 3; 2; 1; 6]
      tmp[10] += 1
    elseif support.oj.R == [1; 2; 3; 3; 1; 6]
      tmp[11] += 1
    elseif support.oj.R == [1; 2; 3; 4; 1; 6]
      tmp[12] += 1
    elseif support.oj.R == [1; 2; 3; 1; 2; 6]
      tmp[13] += 1
    elseif support.oj.R == [1; 2; 3; 2; 2; 6]
      tmp[14] += 1
    elseif support.oj.R == [1; 2; 3; 3; 2; 6]
      tmp[15] += 1
    elseif support.oj.R == [1; 2; 3; 4; 2; 6]
      tmp[16] += 1
    elseif support.oj.R == [1; 2; 3; 1; 3; 6]
      tmp[17] += 1
    elseif support.oj.R == [1; 2; 3; 2; 3; 6]
      tmp[18] += 1
    elseif support.oj.R == [1; 2; 3; 4; 3; 6]
      tmp[19] += 1
    elseif support.oj.R == [1; 2; 3; 4; 4; 6]
      tmp[20] += 1
    elseif support.oj.R == [1; 2; 3; 1; 5; 6]
      tmp[21] += 1
    elseif support.oj.R == [1; 2; 3; 2; 5; 6]
      tmp[22] += 1
    elseif support.oj.R == [1; 2; 3; 3; 5; 6]
      tmp[23] += 1
    elseif support.oj.R == [1; 2; 3; 4; 5; 6]
      tmp[24] += 1
    end
  end

  tmp /= N

  @test isapprox(tmp, pr, rtol=0.01)

  nothing
end

test_crossover()


text:
module TestUtils
using GaussianProcesses
using Test, Random

Random.seed!(1)

@testset "Utils" begin
    d, n = 5, 10
    logℓ, logσ = rand(d), rand()
    X = rand(d,n)
    kern = SEArd(logℓ, logσ)
    data = GaussianProcesses.KernelData(kern, X)

    @test GaussianProcesses.distance(kern, X, data) ≈ GaussianProcesses.distance(kern, X)
end
end


text:
#wave = cwts[4]; bc = bcs[1]; β=βs[2]; ave = averagingLengths[end]; n = ns[1]; testF = typesOfTestFunctions[1]; inverseType = inversionMethods[3]
#wave = morl; bc = ZPBoundary(); ave = -1; n = 1382; testF = "just Core"; inverseType = DualFrames(); β = 1.5; eOct = 0
@testset "Inversion" begin
    bcs =(ContinuousWavelets.PerBoundary(), ContinuousWavelets.ZPBoundary(),)
    cwts = (ContinuousWavelets.dog2,ContinuousWavelets.morl)
    βs =(1,)
    averagingLengths=(0,)
    extraOctaves=(0,)
    typesOfTestFunctions = ["Doppler"]
    inversionMethods = [NaiveDelta(), PenroseDelta(), DualFrames()]
    ns = (128, 2039)
    @testset "length $n, with type $wave, bc $bc, β=$β, ave=$(ave) ex=$(testF), inv=$(inverseType)" for n in ns, wave in cwts, bc in bcs, β in βs, ave in averagingLengths, testF in typesOfTestFunctions, inverseType in inversionMethods, eOct in extraOctaves
        if !(wave isa Morlet && eOct < 0)
            if testF =="just Core"
                x̂ = zeros(ComplexF64, n>>1+1); x̂[n>>2 .+ (-10:10)] = 5*(11.5 .- abs.(range(-10/2,10/2,length=21))) .* randn(ComplexF64, 21)
                x = irfft(x̂, n) # test function that is definitely solidly in the "safe" range
            else
                x = testfunction(n,testF)
            end

            wav = ContinuousWavelets.wavelet(wave, β=β,boundary=bc,averagingLength = ave, extraOctaves = eOct)
            with_logger(ConsoleLogger(stderr, Logging.Error)) do
                res = ContinuousWavelets.cwt(x, wav)
                xRecon = real.(ContinuousWavelets.icwt(res, wav, inverseType))
                err = norm(xRecon-x)/norm(x)
                @test err < 3 # none of them are wildly off
            end
        end
    end
end


text:
#Test for Gadfly with different chunk options and figure formats
using Weave
using Base.Test

weave("documents/gadfly_formats_test.txt", doctype="tex", plotlib="gadfly")
result = readall(open("documents/gadfly_formats_test.tex"))
ref = readall(open("documents/gadfly_formats_test_ref.tex"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="pandoc", plotlib="gadfly")
result = readall(open("documents/gadfly_formats_test.md"))
ref = readall(open("documents/gadfly_formats_test_pandoc_ref.md"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="pandoc", plotlib="gadfly", fig_ext=".svg")
result = readall(open("documents/gadfly_formats_test.md"))
ref = readall(open("documents/gadfly_formats_test_svg_ref.md"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="github", plotlib="gadfly", fig_ext=".js.svg")
result = readall(open("documents/gadfly_formats_test.md"))
ref = readall(open("documents/gadfly_formats_test_jssvg_ref.md"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="rst", plotlib="gadfly")
result = readall(open("documents/gadfly_formats_test.rst"))
ref = readall(open("documents/gadfly_formats_test_ref.rst"))
@test result == ref


text:
# currently mysortperm is faster than sortperm in julia base
# check also test_sort.jl in the testing folder

export mysortperm,mysortperm!

function my_sortperm_you_should_regularily_check_if_sortperm_in_base_has_become_more_efficient!(ii, A)
    n = length(A)
   # ii = Array(Int,n)

    for i = 1:n
        ii[i] = i
    end

    B = copy(A)
    quicksort!(B, ii)

   # return ii, B
    return nothing
end # function mysortperm!

function mysortperm(A)
    n = length(A)
    ii = Array(Int, n)

    for i = 1:n
        ii[i] = i
    end

    B = copy(A)
    quicksort!(B, ii)

    return ii, B
end # function mysortperm

function quicksort!(A, order, i=1, j=length(A))
# modified from:
# http://rosettacode.org/wiki/Sorting_algorithms/Quicksort#Julia
    if j > i
        if j - i <= 50
              # Insertion sort for small groups is faster than Quicksort
            InsertionSort!(A, order, i, j)
            return A
        end

        # pivot = A[rand(i:j)] # random element of A
        pivot = A[ div(i + j, 2) ]
        left, right = i, j
        while left <= right
            while A[left] < pivot
                left += 1
            end
            while A[right] > pivot
                right -= 1
            end
            if left <= right
                A[left], A[right] = A[right], A[left]
                order[left], order[right] = order[right], order[left]

                left += 1
                right -= 1
            end
        end  # left <= right

        quicksort!(A, order, i,   right)
        quicksort!(A, order, left, j)
    end  # j > i

    return A
end # function quicksort!

function InsertionSort!(A, order, ii=1, jj=length(A))

    for i = ii + 1:jj
        j = i - 1
        temp  = A[i]
        itemp = order[i]

        while true
            if j == ii - 1
                break
            end
            if A[j] <= temp
                break
            end
            A[j + 1] = A[j]
            order[j + 1] = order[j]
            j -= 1
        end

        A[j + 1] = temp
        order[j + 1] = itemp
    end  # i

    return
end # function InsertionSort!