text:
using LinearAlgebra, GeometricMultigrid, Test, Statistics

function setup_2D(n)
    L = zeros(n+2,n+2,2); L[3:n+1,2:n+1,1] .= 1; L[2:n+1,3:n+1,2] .= 1;
    A = Poisson(L)
    x = FieldVector(Float64[i-1 for i ∈ 1:n+2, j ∈ 1:n+2])
    A,x
end
error(x) = norm(x.-Statistics.mean(x))

@testset "FieldVector.jl" begin
    A,x = setup_2D(3)
    @test x[1] == 1
    @test x[CartesianIndex(1,1)] == 0
    @test norm(x) == √42
    b = zero(x)
    @test typeof(b) <: FieldVector
    @test b == zeros(9)
    @loop b[I] = 1/x[I]
    @test !any(isnan.(b.data))
    @test dot(b,x) == 9
end

@testset "FieldMatrix.jl" begin
    A,x = setup_2D(3)
    @test diag(A) == -Float64[2,3,2,3,4,3,2,3,2]
    @test det(A) < 2eps(eltype(A))
    b = A*x
    @test b == repeat(Float64[1,0,-1],3)
    @test dot(b,A,x) == dot(b,b)
    @test pinv(Matrix(A))*(A*x) ≈ -b

    using IterativeSolvers
    A,x = setup_2D(37)
    @test error(x.-cg(A,A*x))<1e-8
end

@testset "SolveState.jl" begin
    A,x = setup_2D(32)
    @test GeometricMultigrid.residual(A,x,A*x) == zero(x)
    y,_ = GeometricMultigrid.gs(A,A*x,reltol=1e-12)
    @test error(x.-y)<1e-8
end

@testset "MultiGrid.jl" begin
    I = CartesianIndex(4,3,2)
    @test all(GeometricMultigrid.down(J)==I for J ∈ GeometricMultigrid.up(I))
    A,x = setup_2D(4)
    st = mg_state(A,x,x)
    @test diag(st.child.A) == -2ones(4)
    @test isnothing(st.child.child)
    A,x = setup_2D(32)
    y,_ = mg(A,A*x,reltol=1e-12)
    @test error(x.-y)<1e-8
end


text:
var a = "global a";
var b = "global b";
var c = "global c";
{
    var a = "outer a";
    var b = "outer b";
    {
        var a = "inner a";
        print a;
        print b;
        print c;
    }
    print a;
    print b;
    print c;
}
print a;
print b;
print c;
if(true) {
    print "hello";
    print "This is inside a block!";
}
if(false)
    print "Branch 1 executed";
else
    print "Branch 2 executed";

print "hi" or 2; // "hi".
print nil or "yes"; // "yes".

var i = 1;
while (i < 5) {
    print i;
    i = i + 1;
}

print "Beginning for loop...";
for (var i = 0; i < 10; i = i + 1) print i;

fun sayHi(first, last) {
    print "Hi, " + first + " " + last + "!";
}
var zz = "Dear";
var yy = "Reader";
sayHi(zz, yy);
sayHi("Dear", "Reader");

//
// Fibonacci Numbers!!
//
fun fibonacci(n) {
    if (n <= 1) return n;
    return fibonacci(n - 2) + fibonacci(n - 1);
}

for (var i = 0; i < 3; i = i + 1) {
    print fibonacci(i);
}

ext: jl


text:
export cgs, cgs!, mgs, mgs!

"""
KrylovMethods.cgs = KrylovMethods.cgs!(copy(V))

	Classical Gram Schmidt orthogonalization.

"""
cgs(V) = cgs!(copy(V))

"""
function KrylovMethods.cgs!

Inplace Classical Gram Schmidt orthogonalization.

Reference: page 254 in Golub and Van Loan, Matrix Computation, 4th edition.

Input:

	V::Array  - m by n matrix of full rank m<=n

Output:

	V::Array  - m-by-m unitary matrix
	R::Array  - m-by-n upper triangular matrix

"""
function cgs!{T}(V::Array{T})

	m,n    = size(V)
	R      = zeros(T,n,n)
	R[1,1] = norm(V[:,1])
	V[:,1] ./= R[1,1]
	for k=2:n
	    R[1:k-1,k] = V[:,1:k-1]'*V[:,k]
	    V[:,k]     = V[:,k]-V[:,1:k-1]*R[1:k-1,k]
	    R[k,k]     = norm(V[:,k])
	    V[:,k]   ./= R[k,k]
	end
	return V,R
end

"""
KrylovMethods.mgs = KrylovMethods.mgs!(copy(V))

	Modified Gram Schmidt orthogonalization.

"""
mgs(V) = mgs!(copy(V))

"""
function KrylovMethods.mgs!

Inplace Modified Gram Schmidt orthogonalization.

Reference: page 255 in Golub and Van Loan, Matrix Computation, 4th edition.

Input:

	V::Array  - m by n matrix of full rank m<=n

Output:

	V::Array  - m-by-m unitary matrix
	R::Array  - m-by-n upper triangular matrix

"""
function mgs!{T}(V::Array{T})
	m,n    = size(V)
	R      = zeros(T,n,n)
	for k=1:n
		R[k,k]    = norm(V[:,k])
		V[:,k]    ./= R[k,k]
		for j=k+1:n
			R[k,j] = dot(V[:,k],V[:,j])
			V[:,j] -= R[k,j]*V[:,k]
		end
	end
	return V,R
end


text:
# This file was generated, do not modify it. # hide
X = rand(Float32, dim, 10^5)
Q = rand(Float32, dim, 10^3)
k = 10
nothing # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
module InventoryManagementWebApp

using Dash
using CSV, DataFrames, Base64, JSON
using InventoryManagement, Distributions

export build_app

include("utils.jl")
include("callbacks.jl")
include("app.jl")

end # module


text:
# truncation.jl
#
# Implements truncation schemes for truncating a tensor with svd, leftorth or rightorth
abstract type TruncationScheme end

struct NoTruncation <: TruncationScheme
end
notrunc() = NoTruncation()

struct TruncationError{T<:Real} <: TruncationScheme
    ϵ::T
end
truncerr(epsilon::Real) = TruncationError(epsilon)

struct TruncationDimension <: TruncationScheme
    dim::Int
end
truncdim(d::Int) = TruncationDimension(d)

struct TruncationSpace{S<:ElementarySpace} <: TruncationScheme
    space::S
end
truncspace(space::ElementarySpace) = TruncationSpace(space)

struct TruncationCutoff{T<:Real} <: TruncationScheme
    ϵ::T
    add_back::Int
end
truncbelow(epsilon::Real, add_back::Int=0) = TruncationCutoff(epsilon, add_back)

# For a single vector
function _truncate!(v::AbstractVector, ::NoTruncation, p::Real = 2)
    return v, zero(real(eltype(v)))
end

function _truncate!(v::AbstractVector, trunc::TruncationError, p::Real = 2)
    S = real(eltype(v))
    truncerr = zero(S)
    dmax = length(v)
    dtrunc = dmax
    while dtrunc > 0
        dtrunc -= 1
        prevtruncerr = truncerr
        truncerr = norm(view(v, dtrunc+1:dmax), p)
        if truncerr > trunc.ϵ
            dtrunc += 1
            truncerr = prevtruncerr
            break
        end
    end
    resize!(v, dtrunc)
    return v, truncerr
end

function _truncate!(v::AbstractVector, trunc::TruncationDimension, p::Real = 2)
    S = real(eltype(v))
    dtrunc = min(length(v), trunc.dim)
    truncerr = norm(view(v, dtrunc+1:length(v)), p)
    resize!(v, dtrunc)
    return v, truncerr
end

_truncate!(v::AbstractVector, trunc::TruncationSpace, p::Real = 2) =
    _truncate!(v, truncdim(dim(trunc.space)), p)

function _truncate!(v::AbstractVector, trunc::TruncationCutoff, p::Real = 2)
    S = real(eltype(v))
    dtrunc = findlast(Base.Fix2(>, trunc.ϵ), v)
    if dtrunc === nothing
        dtrunc = 0
    end
    dtrunc = min(dtrunc + trunc.add_back, length(v))
    truncerr = norm(view(v, dtrunc+1:length(v)), p)
    resize!(v, dtrunc)
    return v, truncerr
end

# For SectorDict
const SectorVectorDict{I<:Sector} = SectorDict{I, <:AbstractVector}

function _findnexttruncvalue(V::SectorVectorDict{I},
                                truncdim::SectorDict{I, Int}, p::Real) where {I<:Sector}
    S = real(eltype(valtype(V)))
    q = convert(S, p)
    it = keys(V)
    next = iterate(it)
    next === nothing && nothing
    c, s = next
    while truncdim[c] == 0
        next = iterate(it, s)
        next === nothing && return nothing
        c, s = next
    end
    cmin = c
    vmin::S = convert(S, dim(c))^inv(q)*V[c][truncdim[c]]
    next = iterate(it, s)
    while next !== nothing
        c, s = next
        if truncdim[c] > 0
            v = dim(c)^inv(q)*V[c][truncdim[c]]
            if v < vmin
                cmin, vmin = c, v
            end
        end
        next = iterate(it, s)
    end
    return cmin
end

function _truncate!(V::SectorVectorDict, ::NoTruncation, p = 2)
    S = real(eltype(valtype(V)))
    return V, zero(S)
end
function _truncate!(V::SectorVectorDict, trunc::TruncationError, p = 2)
    I = keytype(V)
    S = real(eltype(valtype(V)))
    truncdim = SectorDict{I, Int}(c=>length(v) for (c, v) in V)
    truncerr = zero(S)
    while true
        cmin = _findnexttruncvalue(V, truncdim, p)
        cmin === nothing && break
        truncdim[cmin] -= 1
        prevtruncerr = truncerr
        truncerr = _norm((c=>view(v, truncdim[c]+1:length(v)) for (c, v) in V), p, zero(S))
        if truncerr > trunc.ϵ
            truncdim[cmin] += 1
            truncerr = prevtruncerr
            break
        end
    end
    for (c, v) in V
        resize!(v, truncdim[c])
    end
    return V, truncerr
end
function _truncate!(V::SectorVectorDict, trunc::TruncationDimension, p = 2)
    I = keytype(V)
    S = real(eltype(valtype(V)))
    truncdim = SectorDict{I, Int}(c=>length(v) for (c, v) in V)
    while sum(dim(c)*d for (c, d) in truncdim) > trunc.dim
        cmin = _findnexttruncvalue(V, truncdim, p)
        cmin === nothing && break
        truncdim[cmin] -= 1
    end
    truncerr = _norm((c=>view(v, truncdim[c]+1:length(v)) for (c, v) in V), p, zero(S))
    for (c, v) in V
        resize!(v, truncdim[c])
    end
    return V, truncerr
end

function _truncate!(V::SectorVectorDict, trunc::TruncationSpace, p = 2)
    I = keytype(V)
    S = real(eltype(valtype(V)))
    truncdim = SectorDict{I, Int}(c=>min(length(v), dim(trunc.space, c)) for (c, v) in V)
    truncerr = _norm((c=>view(v, truncdim[c]+1:length(v)) for (c, v) in V), p, zero(S))
    for c in keys(V)
        resize!(V[c], truncdim[c])
    end
    return V, truncerr
end

function _truncate!(V::SectorVectorDict, trunc::TruncationCutoff, p = 2)
    I = keytype(V)
    S = real(eltype(valtype(V)))
    truncdim = SectorDict{I, Int      }(c=>length(v) for (c, v) in V)
    next_val = SectorDict{I, Array{S} }(c=>[zero(S)] for (c, v) in V)
    for (c, v) in V
        newdim = findlast(Base.Fix2(>, trunc.ϵ), v)
        if newdim === nothing
            truncdim[c] = 0
	    next_val[c] = v
	else
	    truncdim[c] = newdim
            next_val[c] = v[newdim+1:end]
        end
    end
    for i in 1:trunc.add_back
        key_max = argmax(next_val)
	length(next_val[key_max]) == 0 && break
        truncdim[key_max] += 1
        next_val[key_max] = next_val[key_max][2:end]
    end
    truncerr = _norm((c=>view(v, truncdim[c]+1:length(v)) for (c, v) in V), p, zero(S))
    for (c, v) in V
        resize!(v, truncdim[c])
    end
    return V, truncerr
end

# Combine truncations
struct MultipleTruncation{T<:Tuple{Vararg{<:TruncationScheme}}} <: TruncationScheme
    truncations::T
end
Base.:&(a::MultipleTruncation, b::MultipleTruncation) =
    MultipleTruncation((a.truncations..., b.truncations...))
Base.:&(a::MultipleTruncation, b::TruncationScheme) =
    MultipleTruncation((a.truncations..., b))
Base.:&(a::TruncationScheme, b::MultipleTruncation) =
    MultipleTruncation((a, b.truncations...))
Base.:&(a::TruncationScheme, b::TruncationScheme) = MultipleTruncation((a, b))

function _truncate!(v, trunc::MultipleTruncation, p::Real = 2)
    v, truncerrs = __truncate!(v, trunc.truncations, p)
    return v, norm(truncerrs, p)
end
function __truncate!(v, trunc::Tuple{Vararg{<:TruncationScheme}}, p::Real = 2)
    v, truncerr1 = _truncate!(v, first(trunc), p)
    v, truncerrtail = __truncate!(v, tail(trunc), p)
    return v, (truncerr1, truncerrtail...)
end
function __truncate!(v, trunc::Tuple{<:TruncationScheme}, p::Real = 2)
    v, truncerr1 = _truncate!(v, first(trunc), p)
    return v, (truncerr1,)
end


text:
#=
Tabulated Gaussian quadrature rules from GMSH.
Obtained from `https://gitlab.onelab.info/gmsh/gmsh/-/blob/gmsh_4_7_0/Numeric/GaussQuadratureTri.cpp`
and `https://gitlab.onelab.info/gmsh/gmsh/-/blob/gmsh_4_7_0/Numeric/GaussQuadratureTet.cpp`.
=#

include("quadrulestables_tetrahedron.jl")
include("quadrulestables_triangle.jl")

##
# Dictionaries that contains the quadrature rules
# for various `AbstractReferenceShape`. The dictionary
# key represents the number of quadrature nodes.

const TRIANGLE_GAUSS_QRULES = Dict(
    1 => TRIANGLE_G1N1,
    3 => TRIANGLE_G2N3,
    4 => TRIANGLE_G3N4,
    6 => TRIANGLE_G4N6,
    7 => TRIANGLE_G5N7,
    12 => TRIANGLE_G6N12,
    13 => TRIANGLE_G7N13,
    16 => TRIANGLE_G8N16,
    19 => TRIANGLE_G9N19,
    25 => TRIANGLE_G10N25,
    33 => TRIANGLE_G12N33
)

# map a desired quadrature order to the number of nodes
const TRIANGLE_GAUSS_ORDER_TO_NPTS = Dict(
    1 => 1,
    2 => 3,
    3 => 4,
    4 => 6,
    5 => 7,
    6 => 12,
    7 => 13,
    8 => 16,
    9 => 19,
    10 => 25,
    12 => 33
)
const TRIANGLE_GAUSS_NPTS_TO_ORDER = Dict((v,k) for (k,v) in TRIANGLE_GAUSS_ORDER_TO_NPTS)

const TETAHEDRON_GAUSS_QRULES = Dict(
    1 => TETAHEDRON_G1N1,
    4 => TETAHEDRON_G2N4,
    5 => TETAHEDRON_G3N5,
    11 => TETAHEDRON_G4N11,
    14 => TETAHEDRON_G5N14,
    24 => TETAHEDRON_G6N24,
    31 => TETAHEDRON_G7N31,
    43 => TETAHEDRON_G8N43
)

# map a desired quadrature order to the number of nodes
const TETRAHEDRON_GAUSS_ORDER_TO_NPTS = Dict(
    1 => 1,
    2 => 4,
    3 => 5,
    4 => 11,
    5 => 14,
    6 => 24,
    7 => 31,
    8 => 43
)
const TETRAHEDRON_GAUSS_NPTS_TO_ORDER = Dict((v,k) for (k,v) in TETRAHEDRON_GAUSS_ORDER_TO_NPTS)

##
const GAUSS_QRULES = Dict(
    ReferenceTriangle => TRIANGLE_GAUSS_QRULES,
    ReferenceTetrahedron => TETAHEDRON_GAUSS_QRULES
)

"""
    _get_qnodes_and_qweights(R::Type{<:AbstractReferenceShape{D}}, N) where D

Returns the `N`-point symmetric gaussian qnodes and qweights `(x, w)` for integration over `R`.
"""
function _get_gauss_qnodes_and_qweights(R::Type{<:AbstractReferenceShape{D}}, N) where D
    if !haskey(GAUSS_QRULES, R) || !haskey(GAUSS_QRULES[R], N)
        error("quadrature rule not found")
    end
    qrule = GAUSS_QRULES[R][N]
    @assert length(qrule) == N
    # qnodes
    qnodestype = SVector{N,SVector{D,Float64}}
    x = qnodestype([q[1] for q in qrule])
    # qweights
    qweightstype = SVector{N,Float64}
    w = qweightstype([q[2] for q in qrule])
    return x, w
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
# [1] M. R. Teague, “Image analysis via the general theory of moments*,” Journal of the Optical Society of America, vol. 70, no. 8, p. 920, Aug. 1980.
# https://doi.org/10.1364/josa.70.000920
function measure_feature(property::RegionEllipse, t::IndexedTable, labels::AbstractArray, N::Int = 0)
    init = StepRange(typemax(Int), -1, -typemax(Int))
    ℳ₀₀ = zeros(N)
    ℳ₁₀ = zeros(N)
    ℳ₀₁ = zeros(N)
    ℳ₁₁ = zeros(N)
    ℳ₂₀ = zeros(N)
    ℳ₀₂ = zeros(N)
    # TODO Document discretization model that underpins these moment computations.
    for i in CartesianIndices(labels)
        l = labels[i]
        if  l != 0
            y, x = i.I
            ℳ₀₀[l] += 1
            ℳ₁₀[l] += x
            ℳ₀₁[l] += y
            ℳ₁₁[l] += x*y
            ℳ₂₀[l] += x^2 + 1/12
            ℳ₀₂[l] += y^2 + 1/12
        end
    end
    t = pushcol(t, :M₀₀ => ℳ₀₀, :M₁₀ => ℳ₁₀, :M₀₁ => ℳ₀₁, :M₁₁ => ℳ₁₁, :M₂₀ => ℳ₂₀, :M₀₂ => ℳ₀₂)
    fill_properties(property, t)
end

function fill_properties(property::RegionEllipse, t::IndexedTable)
    t = property.centroid ? compute_centroid(t) : t
    t = property.semi_axes ? compute_semi_axes(t) : t
    t = property.orientation ? compute_orientation(t) : t
    t = property.eccentricity ? compute_eccentricity(t) : t
end

function compute_centroid(t::IndexedTable)
    @transform t {centroid = (:M₀₁ / :M₀₀, :M₁₀ / :M₀₀)}
end

function compute_semi_axes(t::IndexedTable)
    semi_axes = select(t, (:M₀₀,  :M₁₀, :M₀₁, :M₁₁, :M₂₀, :M₀₂) => row -> compute_semi_axes(row...))
    t = pushcol(t, :semi_axes => semi_axes)
end

function compute_semi_axes(M₀₀::Real, M₁₀::Real, M₀₁::Real, M₁₁::Real, M₂₀::Real, M₀₂::Real)
    μ′₂₀ = (M₂₀ / M₀₀) - (M₁₀ / M₀₀)^2
    μ′₀₂ = (M₀₂ / M₀₀) - (M₀₁ / M₀₀)^2
    μ′₁₁ = (M₁₁ / M₀₀) - ((M₁₀ / M₀₀) * (M₀₁ / M₀₀))

    # See Equations (7) and (8) in [1].
    l₁ = sqrt((μ′₂₀ + μ′₀₂ + sqrt(4 * μ′₁₁^2 +  (μ′₂₀ - μ′₀₂)^2)) / (1 / 2))
    l₂ = sqrt((μ′₂₀ + μ′₀₂ - sqrt(4 * μ′₁₁^2 +  (μ′₂₀ - μ′₀₂)^2)) / (1 / 2))
    min(l₁, l₂), max(l₁, l₂)
end

function compute_orientation(t::IndexedTable)
    orientation = select(t, (:M₀₀,  :M₁₀, :M₀₁, :M₁₁, :M₂₀, :M₀₂) => row -> compute_orientation(row...))
    t = pushcol(t, :orientation => orientation)
end

function compute_eccentricity(t::IndexedTable)
    eccentricity = select(t, (:M₀₀,  :M₁₀, :M₀₁, :M₁₁, :M₂₀, :M₀₂) => row -> compute_eccentricity(row...))
    t = pushcol(t, :eccentricity => eccentricity)
end

function compute_eccentricity(M₀₀::Real, M₁₀::Real, M₀₁::Real, M₁₁::Real, M₂₀::Real, M₀₂::Real)
    b, a = compute_semi_axes(M₀₀, M₁₀, M₀₁, M₁₁, M₂₀, M₀₂)
    e = sqrt(1 - (b/a)^2)
end

# [1] M. R. Teague, “Image analysis via the general theory of moments*,” Journal of the Optical Society of America, vol. 70, no. 8, p. 920, Aug. 1980.
# https://doi.org/10.1364/josa.70.000920
function compute_orientation(M₀₀::Real, M₁₀::Real, M₀₁::Real, M₁₁::Real, M₂₀::Real, M₀₂::Real)
    μ′₂₀ = M₂₀ / M₀₀ - (M₁₀ / M₀₀)^2
    μ′₀₂ = M₀₂ / M₀₀ - (M₀₁ / M₀₀)^2
    μ′₁₁ = M₁₁ / M₀₀ - (M₁₀ / M₀₀) * (M₀₁ / M₀₀)

    # Ellipse tilt angle for various cases of signs of the second moments [1].
    θ = 0.0
    if μ′₂₀ - μ′₀₂ == 0 && μ′₁₁  == 0
        θ = 0.0
    elseif μ′₂₀ - μ′₀₂ == 0 && μ′₁₁  > 0
        θ = 45.0
    elseif μ′₂₀ - μ′₀₂ == 0 && μ′₁₁  < 0
        θ = -45.0
    elseif μ′₂₀ - μ′₀₂ > 0 && μ′₁₁ == 0
        θ = 0.0
    elseif μ′₂₀ - μ′₀₂ < 0 && μ′₁₁ == 0
        θ = -90.0
    elseif μ′₂₀ - μ′₀₂ > 0 && μ′₁₁ > 0
        # 0 < θ < 45
        ξ = 2*μ′₁₁ / (μ′₂₀ - μ′₀₂)
        θ = (1/2) * atand(ξ)
    elseif μ′₂₀ - μ′₀₂ > 0 && μ′₁₁ < 0
        # -45 < θ < 0
        ξ = 2*μ′₁₁ / (μ′₂₀ - μ′₀₂)
        θ = (1/2) * atand(ξ)
    elseif μ′₂₀ - μ′₀₂ < 0 && μ′₁₁ > 0
        # 45 < θ < 90
        ξ = 2*μ′₁₁ / (μ′₂₀ - μ′₀₂)
        θ = (1/2) * atand(ξ) + 90.0
    elseif μ′₂₀ - μ′₀₂ < 0 && μ′₁₁ < 0
        # -90 < θ < -45
        ξ = 2*μ′₁₁ / (μ′₂₀ - μ′₀₂)
        θ = (1/2) * atand(ξ) - 90.0
    end
    θ
end


text:
using HypercubeTransform

function _split_conditional(lj)
    tc = Base.invokelatest(HypercubeTransform.ascube, lj)
    pr = Soss.prune(lj.model, Soss.observed(lj)...)
    cpr = pr(argvals(lj))
    base = transform(tc, fill(0.5, dimension(tc)))
    med, unflatten = ParameterHandling.flatten(base)
    function lklhd(x::Vector{T}) where {T}
        θ = unflatten(x)
        ℓ::T =  logdensity(lj, θ)::T - logdensity(cpr, θ)::T
        return convert(T,isnan(ℓ) ? -1e10 : ℓ)
    end

    prt = x->first(ParameterHandling.flatten(transform(tc, x)))
    return lklhd, prt, tc, unflatten
end

function create_joint(model,
                      ampobs::Comrade.EHTObservation{F,A};
                      amppriors=(AA=0.1,AP=0.1,AZ=0.1,LM=0.2,JC=0.1,PV=0.1,SM=0.1, SP=0.1)
                      ) where {F, A<:Comrade.EHTVisibilityAmplitudeDatum}
    uamp = Comrade.getdata(ampobs, :u)
    vamp = Comrade.getdata(ampobs, :v)
    bl = Comrade.getdata(ampobs, :baseline)
    s1 = first.(bl)
    s2 = last.(bl)
    st = Tuple(unique(vcat(s1,s2)))
    gpriors = select(amppriors, st)
    spriors = values(gpriors)
    erramp = Comrade.getdata(ampobs, :error)
    amps = Comrade.getdata(ampobs, :amp)

    joint = va(
                image=model,
                gamps=gamps(spriors=spriors,stations=st),
                uamp=μas2rad.(uamp),
                vamp=μas2rad.(vamp),
                s1=s1,
                s2=s2,
                erramp=erramp
                )
    conditioned = (amp = amps,)
    return joint | conditioned
end

"""
    create_ampcpjoint
Takes in a model, ampobs, cpobs file and creates a joint distribution.
Optionally, can fit for the gains. If true the gains are set to unity.
"""
function create_joint(model,
                      ampobs::Comrade.EHTObservation{F,A},
                      cpobs::Comrade.EHTObservation{F,P};
                      amppriors=(AA=0.1,AP=0.1,AZ=0.1,LM=0.2,JC=0.1,PV=0.1,SM=0.1, SP=0.1)
                      ) where {F, A<:Comrade.EHTVisibilityAmplitudeDatum,P<:Comrade.EHTClosurePhaseDatum}
    uamp = Comrade.getdata(ampobs, :u)
    vamp = Comrade.getdata(ampobs, :v)
    bl = Comrade.getdata(ampobs, :baseline)
    s1 = first.(bl)
    s2 = last.(bl)
    st = Tuple(unique(vcat(s1,s2)))
    gpriors = select(amppriors, st)
    spriors = values(gpriors)
    stations = keys(gpriors)
    erramp = Comrade.getdata(ampobs, :error)
    amps = Comrade.getdata(ampobs, :amp)

    u1cp = Comrade.getdata(cpobs, :u1)
    v1cp = Comrade.getdata(cpobs, :v1)
    u2cp = Comrade.getdata(cpobs, :u2)
    v2cp = Comrade.getdata(cpobs, :v2)
    u3cp = Comrade.getdata(cpobs, :u3)
    v3cp = Comrade.getdata(cpobs, :v3)
    errcp = Comrade.getdata(cpobs, :error)
    cps = Comrade.getdata(cpobs, :phase)

    joint = vacp(
                  image=model,
                  gamps=gamps(stations=stations, spriors=spriors,),
                  uamp= μas2rad(uamp),
                  vamp= μas2rad(vamp),
                  s1=s1,
                  s2=s2,
                  erramp=erramp,
                  u1cp = μas2rad(u1cp),
                  v1cp = μas2rad(v1cp),
                  u2cp = μas2rad(u2cp),
                  v2cp = μas2rad(v2cp),
                  u3cp = μas2rad(u3cp),
                  v3cp = μas2rad(v3cp),
                  errcp = errcp
                )
    conditioned = (amp = amps, cphase = cps,)
    return joint | conditioned
end

"""
    create_joint
Takes in a model, visobs file and creates a joint distribution.
Optionally, can fit for the gains. If true the gains are set to unity.
"""
function create_joint(model,
                      visobs::Comrade.EHTObservation{F,A},
                      amppriors,
                      phasepriors
                      ) where {F, A<:Comrade.EHTVisibilityDatum}

    u = Comrade.getdata(visobs, :u)
    v = Comrade.getdata(visobs, :v)
    bl = Comrade.getdata(visobs, :baseline)
    s1 = first.(bl)
    s2 = last.(bl)
    st = Tuple(unique(vcat(s1,s2)))
    gpriors = select(amppriors, st)
    ppriors = select(phasepriors, st)
    sapriors = values(gpriors)
    sppriors = values(ppriors)
    err = Comrade.getdata(visobs, :error)
    visr = Comrade.getdata(visobs, :visr)
    visi = Comrade.getdata(visobs, :visi)

    joint = vis(image = model,
                gamps=gamps(spriors=sapriors, stations=st),
                gphases=gphases(spriors=sppriors, stations=st),
                u=μas2rad(u),
                v=μas2rad(v),
                s1=s1,
                s2=s2,
                err=err
                )
    conditioned = (visr = visr, visi = visi,)

    return joint | conditioned
end

"""
    create_joint
Takes in a model, visobs file and creates a joint distribution.
Optionally, can fit for the gains. If true the gains are set to unity.
"""
function create_joint_wnoise(model,
                      visobs::Comrade.EHTObservation{F,A},
                      amppriors,
                      phasepriors,
                      gainamps,
                      gainphases
                      ) where {F, A<:Comrade.EHTVisibilityDatum}

    u = Comrade.getdata(visobs, :u)
    v = Comrade.getdata(visobs, :v)
    bl = Comrade.getdata(visobs, :baseline)
    s1 = first.(bl)
    s2 = last.(bl)
    st = Tuple(unique(vcat(s1,s2)))
    gpriors = select(amppriors, st)
    ppriors = select(phasepriors, st)
    sapriors = values(gpriors)
    sppriors = values(ppriors)
    err = Comrade.getdata(visobs, :error)
    visr = Comrade.getdata(visobs, :visr)
    visi = Comrade.getdata(visobs, :visi)

    joint = viswnoise(image = model,
                gamps=gamps(spriors=sapriors, stations=st),
                gphases=gphases(spriors=sppriors, stations=st),
                u=μas2rad(u),
                v=μas2rad(v),
                s1=s1,
                s2=s2,
                err=err
                )
    if gainamps && gainphases
        conditioned = (visr = visr, visi = visi,)
    elseif gainamps
        conditioned = (visr = visr, visi = visi,
                        gp = (σ = zeros(length(st)-1),),
                      )
    elseif gainphases
        conditioned = (visr = visr, visi = visi,
                        ga = (σ = ones(length(st)),),
                      )
    end

    return joint | conditioned
end

function create_joint(model,
                     dlca::Comrade.EHTObservation{F,A},
                     dcp::Comrade.EHTObservation{F,P};
                     kwargs...) where {F,
                                       A<:Comrade.EHTLogClosureAmplitudeDatum,
                                       P<:Comrade.EHTClosurePhaseDatum}

    u1a = getdata(dlca, :u1)
    v1a = getdata(dlca, :v1)
    u2a = getdata(dlca, :u2)
    v2a = getdata(dlca, :v2)
    u3a = getdata(dlca, :u3)
    v3a = getdata(dlca, :v3)
    u4a = getdata(dlca, :u4)
    v4a = getdata(dlca, :v4)
    lcamp = getdata(dlca, :amp)
    errcamp = getdata(dlca, :error)

    u1cp = getdata(dcp, :u1)
    v1cp = getdata(dcp, :v1)
    u2cp = getdata(dcp, :u2)
    v2cp = getdata(dcp, :v2)
    u3cp = getdata(dcp, :u3)
    v3cp = getdata(dcp, :v3)
    cp = getdata(dcp, :phase)
    errcp = getdata(dcp, :error)

    m = lcacp(
              image=model,
              u1a=μas2rad(u1a),v1a=μas2rad(v1a),
              u2a=μas2rad(u2a),v2a=μas2rad(v2a),
              u3a=μas2rad(u3a),v3a=μas2rad(v3a),
              u4a=μas2rad(u4a),v4a=μas2rad(v4a),
              errcamp=errcamp,
              u1cp= μas2rad(u1cp),v1cp= μas2rad(v1cp),
              u2cp= μas2rad(u2cp),v2cp= μas2rad(v2cp),
              u3cp= μas2rad(u3cp),v3cp= μas2rad(v3cp),
              errcp=errcp
        )

    conditioned = (lcamp = lcamp, cphase = cp,)
    return m | conditioned
end

abstract type AbstractOptimizer end

"""
    optimize(opt::AbstractOptimizer, cm::Soss.ConditionalModel)
Optimize a Soss conditional model to find the MAP. Usually this
requires bringing in other packages to load. For example to use
BlackBoxOptim you would do

```julia
using ComradeSoss, BlackBoxOptim

opt = BBO()
cm = ... #create log joint
x, logMAP = optimize(opt, cm)
```
"""
function optimize end

abstract type AbstractSampler end
abstract type AbstractMCMC end
abstract type AbstractNested <: AbstractSampler end

Base.@kwdef struct DynestyStatic <: AbstractNested
    nlive::Int = 500
    bound::String = "multi"
    sample::String = "auto"
    walks::Int = 25
    slices::Int = 5
    max_move::Int = 100
end

Base.@kwdef struct DynestyDynamic <: AbstractNested
    bound::String="multi"
    sample::String="auto"
    walks::Int = 25
    slices::Int = 5
    max_move::Int = 100
end

function StatsBase.sample(dy::DynestyDynamic, lj::Soss.ConditionalModel; progress=true, kwargs...)
    lklhd, prt, tc, unflatten = _split_conditional(lj)

    sampler =  dynesty.dynamicsampler.DynamicSampler(lklhd,
                                     prt,
                                     dimension(tc);
                                     bound=dy.bound,
                                     sample=dy.sample,
                                     walks = dy.walks,
                                     slices = dy.slices,
                                     max_move=dy.max_move
                                    )
    sampler.run_nested(;kwargs...)
    res = sampler.results
    samples, weights = res["samples"], exp.(res["logwt"] .- res["logz"][end])
    logz = res["logz"][end]
    logzerr = res["logzerr"][end]
    logl = res["logl"][end]
    stats = (logz=logz, logzerr=logzerr, logl=logl)
    return _create_tv(unflatten, samples, weights), stats
end

"""
    dynesty_sampler(logj; nlive=400, kwargs...)
Takes a problem defined by the logj likelihood object and runs the
dynesty on it using the default options and static sampler.

Returns a chain, state
"""
function StatsBase.sample(dy::DynestyStatic, lj::Soss.ConditionalModel; progress=true, kwargs...)
    lklhd, prt, tc, unflatten = _split_conditional(lj)

    sampler =  dynesty.NestedSampler(lklhd,
                                     prt,
                                     dimension(tc);
                                     nlive=dy.nlive,
                                     bound=dy.bound,
                                     sample=dy.sample,
                                     walks = dy.walks,
                                     slices = dy.slices,
                                     max_move=dy.max_move
                                    )
    sampler.run_nested(;kwargs...)
    res = sampler.results
    samples, weights = res["samples"], exp.(res["logwt"] .- res["logz"][end])
    logz = res["logz"][end]
    logzerr = res["logzerr"][end]
    logl = res["logl"][end]
    stats = (logz=logz, logzerr=logzerr, logl=logl)
    return _create_tv(unflatten, samples, weights), stats
end

function _create_tv(unflatten, samples, weights)
    ts = TupleVector([unflatten(Vector(x)) for x in eachrow(samples)])
    tv = TupleVector(merge(getfield(ts, :data), (weights=weights,)))
    return tv
end


text:
struct BipartiteMatchingInstance{T} <: CombinatorialInstance
    graph::AbstractGraph{T} # Ensured to be bipartite.
    reward::Dict{Edge{T}, Float64}

    partition::Vector{UInt8}
    n_left::Int
    n_right::Int
    vertex_left::Vector{T}
    vertex_right::Vector{T}

    function BipartiteMatchingInstance(
        graph::AbstractGraph{T},
        reward::Dict{Edge{T}, Float64},
    ) where {T}
        # Implicit assumption: all edges in the rewards exist in the graph; all edges in the graph have a reward.

        bmap = bipartite_map(graph)
        if length(bmap) != nv(graph) # Condition from is_bipartite
            error("The input graph is not bipartite!")
        end

        n_left = sum(bmap .== 1)
        n_right = sum(bmap .== 2)

        emap = collect(enumerate(bmap))
        emap1 = collect(filter(v -> v[2] == 1, emap))
        emap2 = collect(filter(v -> v[2] == 2, emap))

        vertex_left = map(first, emap1)
        vertex_right = map(first, emap2)

        @assert length(vertex_left) == n_left
        @assert length(vertex_right) == n_right
        @assert n_left + n_right == nv(graph)

        return new{T}(
            copy(graph),
            copy(reward),
            bmap,
            n_left,
            n_right,
            vertex_left,
            vertex_right,
        )
    end

    function BipartiteMatchingInstance(
        i::BipartiteMatchingInstance{T},
        reward::Dict{Edge{T}, Float64},
    ) where {T}
        # Only for private use. Change the rewards of a given instance.
        return new{T}(
            i.graph,
            reward,
            i.partition,
            i.n_left,
            i.n_right,
            i.vertex_left,
            i.vertex_right,
        )
    end
end

struct BipartiteMatchingSolution{T} <: CombinatorialSolution
    instance::BipartiteMatchingInstance{T}
    solution::Vector{Edge{T}}
    value::Float64
end


text:
using BinDeps
using Compat
using Compat.Sys: iswindows, isapple, isunix

@BinDeps.setup

version = "3370"
baseurl = "ftp://heasarc.gsfc.nasa.gov/software/fitsio/c/"

if isunix()
    archivename = "cfitsio$(version).tar.gz"
elseif iswindows()
    archivename = "cfitsio_MSVC_$(Sys.WORD_SIZE)bit_DLL_$(version).zip"
end

libcfitsio = library_dependency("libcfitsio", aliases=["cfitsio"])
downloadsdir = BinDeps.downloadsdir(libcfitsio)
libdir = BinDeps.libdir(libcfitsio)
srcdir = BinDeps.srcdir(libcfitsio)

if isapple()
    libfilename = "libcfitsio.dylib"
elseif isunix()
    libfilename = "libcfitsio.so"
elseif iswindows()
    libfilename = "cfitsio.dll"
end

# Unix
prefix = joinpath(BinDeps.depsdir(libcfitsio), "usr")
provides(Sources, URI(baseurl*archivename), libcfitsio, unpacked_dir="cfitsio")
provides(BuildProcess,
         (@build_steps begin
             GetSources(libcfitsio)
             @build_steps begin
                 ChangeDirectory(joinpath(srcdir, "cfitsio"))
                 FileRule(joinpath(libdir, libfilename),
                          @build_steps begin
                              `./configure --prefix=$prefix --enable-reentrant`
                              `make shared install`
                          end)
             end
          end), libcfitsio, os = :Unix)

# Windows
provides(BuildProcess,
	 (@build_steps begin
             FileDownloader(baseurl*archivename,
                            joinpath(downloadsdir, archivename))
	     CreateDirectory(srcdir, true)
	     FileUnpacker(joinpath(downloadsdir, archivename), srcdir,
                          joinpath(srcdir, libfilename))
	     CreateDirectory(libdir, true)
	     @build_steps begin
		 ChangeDirectory(srcdir)
	         FileRule(joinpath(libdir, libfilename), @build_steps begin
		          `powershell -Command "cp $(libfilename) $(joinpath(libdir, libfilename))"`
			  end)
             end
	  end), libcfitsio, os = :Windows)

if iswindows()
    push!(BinDeps.defaults, BuildProcess)
end

# OSX
if isapple()
    using Homebrew
    provides(Homebrew.HB, "cfitsio", libcfitsio, os=:Darwin)
end

@BinDeps.install Dict(:libcfitsio => :libcfitsio)

if iswindows()
    pop!(BinDeps.defaults)
end


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Test Nodes" begin
    node = DS.NullNode()
    @test DS.value(node) == nothing
    @test DS.next(node) isa DS.NullNode
    @test DS.prev(node) isa DS.NullNode
    @test repr(node) == ""

    node = DS.ListNode{Int64}(10)
    @test DS.value(node) == 10
    @test DS.next(node) isa DS.NullNode
    @test DS.prev(node) isa DS.NullNode
    @test repr(node) == "10"
end

@testset "Test LinkedList" begin
    # Initial a LinkedList
    T_e = Int64
    T = DS.LinkedList{T_e}
    ll = T()
    @test length(ll) == 0
    @test repr(ll) == "LinkedList{$T_e}([])"

    # Push Test
    for i=1:5
        push!(ll, i)
        @test ll[i] == i
    end
    @test repr(ll) == "LinkedList{$T_e}([1, 2, 3, 4, 5])"

    # Popfirst Test
    for i=1:5
        @test popfirst!(ll) == i
    end
    @test repr(ll) == "LinkedList{$T_e}([])"

    # Pushfirst Test
    for i=6:10
        pushfirst!(ll, i)
        @test ll[1] == i
    end
    @test repr(ll) == "LinkedList{$T_e}([10, 9, 8, 7, 6])"

    # Pop Test
    for i=6:10
        @test pop!(ll) == i
    end
    @test repr(ll) == "LinkedList{$T_e}([])"

    # Setindex Test
    push!(ll, 1)
    ll[1] = 30
    @test ll[1] == 30
    @test repr(ll) == "LinkedList{$T_e}([30])"

    # Insert Test
    ll = T()
    insert!(ll, 30, 1)
    @test ll[1] == 30
    @test repr(ll) == "LinkedList{$T_e}([30])"

    insert!(ll, 10, 1)
    @test ll[1] == 10
    @test repr(ll) == "LinkedList{$T_e}([10, 30])"

    insert!(ll, 20, 2)
    @test ll[2] == 20
    @test repr(ll) == "LinkedList{$T_e}([10, 20, 30])"

    insert!(ll, 40, 4)
    @test ll[4] == 40
    @test repr(ll) == "LinkedList{$T_e}([10, 20, 30, 40])"

    # Delete Test
    @test delete!(ll, 1) == 10
    @test repr(ll) == "LinkedList{$T_e}([20, 30, 40])"

    @test delete!(ll, 2) == 30
    @test repr(ll) == "LinkedList{$T_e}([20, 40])"

    @test delete!(ll, 2) == 40
    @test repr(ll) == "LinkedList{$T_e}([20])"

    @test delete!(ll, 1) == 20
    @test repr(ll) == "LinkedList{$T_e}([])"

    # Throws Test
    ll = T()
    @test_throws BoundsError ll[1]
    insert!(ll, 1, 1)
    @test_throws BoundsError ll[0]
    @test ll[1] == 1
    @test_throws BoundsError ll[2]
    @test_throws BoundsError insert!(ll, 0, 0)
    @test_throws BoundsError delete!(ll, 0)
end


text:
module LoggingProfiler
using IRTools, DataStructures, Colors, Dictionaries
include("logging.jl")
include("whiteblacklist.jl")
include("overdubbing.jl")
include("visualization.jl")
include("luxor.jl")

export @record, @recordfun whitelist!, blacklist!
end

max_stars_repo_name: pevnak/LoggingProfiler

max_issues_repo_name: pevnak/LoggingProfiler

max_forks_repo_name: pevnak/LoggingProfiler


text:
function generate_precompile_execution_file()
 source_code="""
 module sysimagegenerator
    using TestApp
    using PartitionedArrays
    const PArrays = PartitionedArrays
    using MPI

    include("../../mpi/runtests_np4_body.jl")

    prun(all_tests,mpi,(1,1))

    MPI.Finalize()

 end #module
 """
 open("sysimagegenerator.jl","w") do io
   println(io,source_code)
 end
end

using Pkg
using PackageCompiler
generate_precompile_execution_file()
pkgs = Symbol[]
push!(pkgs, :TestApp)

if VERSION >= v"1.4"
    append!(pkgs, [Symbol(v.name) for v in values(Pkg.dependencies()) if v.is_direct_dep],)
else
    append!(pkgs, [Symbol(name) for name in keys(Pkg.installed())])
end

create_sysimage(pkgs,
  sysimage_path=joinpath(@__DIR__,"TestApp.so"),
  precompile_execution_file=joinpath(@__DIR__,"sysimagegenerator.jl"))


text:
using Plots, Einsum, StaticArrays

# Simulation parameters
const Nx, Ny      = 401, 401    # grid points in x and y
const rho0        = 1000    # average density
const Nt          = 600   # number of timesteps
const NL          = 9       # D2Q9 Lattice
const tau         = 0.6    # collision timescale
const omega       = 1/tau # appears in the collision equation, precomputing saves a division.

const cxs         = @SVector[0, 0, 1, 1, 1, 0,-1,-1,-1]
const cys         = @SVector[0, 1, 1, 0,-1,-1,-1, 0, 1]
const opp         = @SVector[1,6,7,8,9,2,3,4,5] # bounce back array, opposite direction. Not used but will be when boundaries are added.
const weights     = @SVector[4/9,1/9,1/36,1/9,1/36,1/9,1/36,1/9,1/36] # sums to 1

##build the struct that will hold the state of the LB system##
struct LatticeState
    rho::Array{Float32, 2}    # macroscale density
    ux::Array{Float32, 2}    # macroscale velocity, x component
    uy::Array{Float32, 2}    # macroscale velocity, y component
    F::Array{Float32, 3}    # particle distribution function
    Feq::Array{Float32, 3}    # equillibrium term for collision
    tmp::Array{Float32, 2}   # place-holder array
end

function LatticeState(Nx,Ny)
    F = ones(Nx,Ny,NL)
    Feq = similar(F)
    ux = zeros(Nx,Ny)
    uy = similar(ux)
    rho = similar(ux)
    tmp = similar(ux)
    sum!(rho, F)

    return LatticeState(rho,ux,uy,F,Feq,tmp)
end

@fastmath @inbounds function apply_drift!(tmp::Array{<:Real},F::Array{<:Real})
    for i in 1:size(F)[3]
        A = @view F[:,:,i]
        circshift!(tmp,A, (cxs[i],0) )
        circshift!(A,tmp, (0,cys[i]) )
    end
end

@fastmath @inbounds function feq_point(rho::Real,ux::Real,uy::Real,cx::Int,cy::Int,weight::Real)
    cu = (cx * ux +  cy * uy)
    return rho*weight*(1. +  3. *cu + 4.5*cu*cu - 1.5*(ux*ux + uy*uy))
end

@fastmath @inbounds function calculate_feq!(Feq::Array{<:Real},rho::Array{<:Real},ux::Array{<:Real},uy::Array{<:Real})
    @einsum Feq[i,j,k] = feq_point(rho[i,j],ux[i,j],uy[i,j],cxs[k],cys[k],weights[k])
end

@fastmath @inbounds function apply_collision!(F::Array{<:Real},Feq::Array{<:Real})
    @. F -= omega .* (F .- Feq)
end

@fastmath @inbounds function calculate_u!(ux::Array{<:Real},uy::Array{<:Real},F::Array{<:Real},rho::Array{<:Real})
    @einsum ux[i,j] = F[i,j,k] * cxs[k] / rho[i,j]
    @einsum uy[i,j] = F[i,j,k] * cys[k] / rho[i,j]
end

@fastmath @inbounds function correct_F!(F::Array{<:Real},rho::Array{<:Real})
    @einsum F[i,j,k] = F[i,j,k] * rho0 / rho[i,j]
end

function iterate_lb!(f::LatticeState, Nt::Int)
    for it in 1:Nt
        apply_drift!(f.tmp,f.F) # Particle streaming step
        sum!(f.rho, f.F) # Calculate fluid density
        calculate_u!(f.ux, f.uy, f.F, f.rho) # Calculate macroscopic velocity
        calculate_feq!(f.Feq, f.rho, f.ux,f.uy) # Calculate Feq then apply the collision step
        apply_collision!(f.F, f.Feq)
    end
end

@fastmath @inbounds function set_initial_condition!(f::LatticeState)
    @. f.rho += [exp(-sqrt((x-Nx/2)^2 + (y-Ny/2)^2)) for x in 1:Nx, y in 1:Ny]; # Initial condition. Modify density with a pulse in the middle
    correct_F!(f.F,f.rho) #correct F to the new density
end

problem = LatticeState(Nx,Ny)
set_initial_condition!(problem)
iterate_lb!(problem,Nt) # Run simulation
Plots.heatmap(problem.rho, c=:redsblues, size=(650,640), clims=(rho0-.2, rho0+0.2), right_margin = 5Plots.mm,aspect_ratio=:equal) # Plot the final density variation
title!("Final density")
savefig("output.png")


text:
export AlgAss, AlgAssElem, AlgGrp, AlgGrpElem, AlgMat, AlgMatElem

abstract type AbsAlgAss{T} <: Ring end

abstract type AbsAlgAssElem{T} <: RingElem end

################################################################################
#
#  AlgAss / AlgAssElem
#
################################################################################

# Associative algebras defined by structure constants (multiplication table)
mutable struct AlgAss{T} <: AbsAlgAss{T}
  base_ring::Ring
  mult_table::Array{T, 3} # e_i*e_j = sum_k mult_table[i, j, k]*e_k
  one::Vector{T}
  has_one::Bool
  iszero::Bool
  iscommutative::Int       # 0: don't know
                           # 1: known to be commutative
                           # 2: known to be not commutative
  basis#::Vector{AlgAssElem{T, AlgAss{T}}
  gens#::Vector{AlgAssElem{T, AlgAss{T}} # "Small" number of algebra generators
  trace_basis_elem::Vector{T}
  issimple::Int

  decomposition#::Vector{Tuple{AlgAss{T}, mor(AlgAss{T}, AlgAss{T})}
  center#Tuple{AlgAss{T}, mor(AlgAss{T}, AlgAss{T})
  maps_to_numberfields
  maximal_order

  # Constructor with default values
  function AlgAss{T}(R::Ring) where {T}
    A = new{T}()
    A.base_ring = R
    A.iscommutative = 0
    A.issimple = 0
    return A
  end

  function AlgAss{T}(R::Ring, mult_table::Array{T, 3}, one::Vector{T}) where {T}
    A = AlgAss{T}(R)
    A.mult_table = mult_table
    A.one = one
    A.has_one = true
    A.iszero = (size(mult_table, 1) == 0)
    return A
  end

  function AlgAss{T}(R::Ring, mult_table::Array{T, 3}) where {T}
    A = AlgAss{T}(R)
    A.mult_table = mult_table
    A.iszero = (size(mult_table, 1) == 0)
    return A
  end
end

mutable struct AlgAssElem{T, S} <: AbsAlgAssElem{T}
  parent::S
  coeffs::Array{T, 1}

  function AlgAssElem{T, S}(A::S) where {T, S}
    z = new{T, S}()
    z.parent = A
    z.coeffs = Array{T, 1}(undef, size(A.mult_table, 1))
    for i = 1:length(z.coeffs)
      z.coeffs[i] = A.base_ring()
    end
    return z
  end

  function AlgAssElem{T, S}(A::AlgAss{T}) where {T, S}
    z = new{T, AlgAss{T}}()
    z.parent = A
    z.coeffs = Array{T, 1}(undef, size(A.mult_table, 1))
    for i = 1:length(z.coeffs)
      z.coeffs[i] = A.base_ring()
    end
    return z
  end

  # This does not make a copy of coeffs
  function AlgAssElem{T, S}(A::S, coeffs::Array{T, 1}) where {T, S}
    z = new{T, AlgAss{T}}()
    z.parent = A
    z.coeffs = coeffs
    return z
  end
end

################################################################################
#
#  AlgGrp / AlgGrpElem
#
################################################################################

# Group rings
mutable struct AlgGrp{T, S, R} <: AbsAlgAss{T}
  base_ring::Ring
  group::S
  group_to_base::Dict{R, Int}
  one::Vector{T}
  basis#::Vector{AlgAssElem{T, AlgAss{T}}
  gens#::Vector{AlgAssElem{T, AlgAss{T}} # "Small" number of algebra generators
  mult_table::Array{Int, 2} # b_i * b_j = b_(mult_table[i, j])
  iscommutative::Int
  trace_basis_elem::Array{T, 1}
  issimple::Int

  decomposition
  center
  maps_to_numberfields
  maximal_order

  function AlgGrp(K::Ring, G::GrpAbFinGen)
    A = AlgGrp(K, G, op = +)
    A.iscommutative = true
    return A
  end

  function AlgGrp(K::Ring, G; op = *)
    A = new{elem_type(K), typeof(G), elem_type(G)}()
    A.iscommutative = 0
    A.issimple = 0
    A.base_ring = K
    A.group = G
    d = Int(order(G))
    A.group_to_base = Dict{elem_type(G), Int}()
    A.mult_table = zeros(Int, d, d)

    for (i, g) in enumerate(G)
      A.group_to_base[deepcopy(g)] = i
    end

    v = Vector{elem_type(K)}(undef, d)
    for i in 1:d
      v[i] = zero(K)
    end
    v[1] = one(K)

    A.one = v

    for (i, g) in enumerate(G)
      for (j, h) in enumerate(G)
        l = op(g, h)
        A.mult_table[i, j] = A.group_to_base[l]
      end
    end

    @assert all(A.mult_table[1, i] == i for i in 1:dim(A))

    return A
  end
end

mutable struct AlgGrpElem{T, S} <: AbsAlgAssElem{T}
  parent::S
  coeffs::Array{T, 1}

  function AlgGrpElem{T, S}(A::S) where {T, S}
    z = new{T, S}()
    z.parent = A
    z.coeffs = Array{T, 1}(undef, size(A.mult_table, 1))
    for i = 1:length(z.coeffs)
      z.coeffs[i] = A.base_ring()
    end
    return z
  end

  function AlgGrpElem{T, S}(A::S, g::U) where {T, S, U}
    return A[A.group_to_base[g]]
  end

  # This does not make a copy of coeffs
  function AlgGrpElem{T, S}(A::S, coeffs::Array{T, 1}) where {T, S}
    z = new{T, S}()
    z.parent = A
    z.coeffs = coeffs
    return z
  end
end

################################################################################
#
#  AbsAlgAssIdl
#
################################################################################

# S is the type of the algebra, T = elem_type(S) and U is the type of matrices
# over the base ring of the algebra
mutable struct AbsAlgAssIdl{S, T, U}
  algebra::S
  basis::Vector{T}
  basis_mat::U

  isleft::Int                      # 0 Not known
                                   # 1 Known to be a left ideal
                                   # 2 Known not to be a left ideal
  isright::Int                     # as for isleft

  iszero::Int

  function AbsAlgAssIdl{S, T, U}(A::S) where {S, T, U}
    I = new{S, T, U}()
    I.algebra = A
    I.isleft = 0
    I.isright = 0
    I.iszero = 0
    return I
  end

  function AbsAlgAssIdl{S, U}(A::S, M::U) where {S, U}
    I = new{S, elem_type(S), U}()
    I.algebra = A
    I.basis_mat = M
    I.isleft = 0
    I.isright = 0
    I.iszero = 0
    return I
  end
end

################################################################################
#
#  AlgAssAbsOrd / AlgAssAbsOrdElem / AlgAssAbsOrdIdl / AlgAssAbsOrdFracIdl
#
################################################################################

# Orders in algebras over the rationals
mutable struct AlgAssAbsOrd{S, T} <: Ring
  algebra::S                       # Algebra containing the order
  dim::Int
  basis#::Vector{AlgAssAbsOrdElem{S, T}}
  basis_alg::Vector{T}             # Basis as array of elements of the algebra
  basis_mat::FakeFmpqMat           # Basis matrix of order wrt basis of the algebra
  basis_mat_inv::FakeFmpqMat       # Inverse of basis matrix
  gen_index::fmpq                  # The det of basis_mat_inv as fmpq
  index::fmpz                      # The det of basis_mat_inv
                                   # (this is the index of the equation order
                                   #  in the given order)
  disc::fmpz                       # Discriminant

  ismaximal::Int                   # 0 Not known
                                   # 1 Known to be maximal
                                   # 2 Known to not be maximal

  #trace_mat::fmpz_mat              # The reduced trace matrix (if known)
  trred_matrix::fmpz_mat

  picard_group#::MapPicardGrp

  function AlgAssAbsOrd{S}(A::S) where {S}
    O = new{S, elem_type(S)}()
    O.algebra = A
    O.dim = dim(A)
    O.ismaximal = 0
    return O
  end

  function AlgAssAbsOrd{S, T}(A::S, basis::Vector{T}) where {S, T}
    # "Default" constructor with default values.
    r = AlgAssAbsOrd{S}(A)
    @assert length(basis) == r.dim
    r.basis_alg = basis
    r.basis_mat = basis_mat(basis)
    return r
  end

  function AlgAssAbsOrd{S}(A::S, basis_mat::FakeFmpqMat) where {S}
    r = AlgAssAbsOrd{S}(A)
    d = dim(A)
    r.basis_mat = basis_mat
    r.basis_alg = Vector{elem_type(S)}(undef, nrows(basis_mat))
    for i in 1:d
      r.basis_alg[i] = elem_from_mat_row(A, basis_mat.num, i, basis_mat.den)
    end
    return r
  end
end

mutable struct AlgAssAbsOrdElem{S, T} <: RingElem
  elem_in_algebra::T
  coordinates::Vector{fmpz}
  has_coord::Bool # needed for mul!
  parent::AlgAssAbsOrd{S, T}

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}) where {S, T}
    z = new{S, T}()
    z.parent = O
    z.elem_in_algebra = algebra(O)()
    z.coordinates = Vector{fmpz}(undef, degree(O))
    z.has_coord = false
    return z
  end

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}, a::T) where {S, T}
    z = new{S, T}()
    z.elem_in_algebra = a
    z.parent = O
    z.coordinates = Vector{fmpz}(undef, degree(O))
    z.has_coord = false
    return z
  end

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}, arr::Vector{fmpz}) where {S, T}
    z = new{S, T}()
    z.elem_in_algebra = dot(O.basis_alg, arr)
    z.coordinates = arr
    z.parent = O
    z.has_coord = true
    return z
  end

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}, a::T, arr::Vector{fmpz}) where {S, T}
    z = new{S, T}()
    z.parent = O
    z.elem_in_algebra = a
    z.coordinates = arr
    z.has_coord = true
    return z
  end
end

mutable struct AlgAssAbsOrdIdl{S, T}
  order::AlgAssAbsOrd                     # Order containing it
  basis::Vector{AlgAssAbsOrdElem{S, T}}   # Basis of the ideal as array of elements of the order
  basis_mat::fmpz_mat                     # Basis matrix of ideal wrt basis of the order
  basis_mat_inv::FakeFmpqMat
  gens::Vector{AlgAssAbsOrdElem{S, T}}    # Generators of the ideal

  isleft::Int                      # 0 Not known
                                   # 1 Known to be a left ideal
                                   # 2 Known not to be a left ideal
  isright::Int                     # as for isleft

  iszero::Int                             # 0: don't know, 1: known to be zero, 2: known to be not zero

  function AlgAssAbsOrdIdl{S, T}(O::AlgAssAbsOrd{S, T}) where {S, T}
    r = new{S, T}()
    r.order = O
    r.isleft = 0
    r.isright = 0
    r.iszero = 0
    return r
  end

  function AlgAssAbsOrdIdl{S, T}(O::AlgAssAbsOrd{S, T}, M::fmpz_mat) where {S, T}
    r = new{S, T}()
    r.order = O
    d = O.dim
    r.basis_mat = M
    r.isleft = 0
    r.isright = 0
    r.iszero = 0
    return r
  end
end

mutable struct AlgAssAbsOrdIdlSet{S, T}
  order::AlgAssAbsOrd{S, T}

  function AlgAssAbsOrdIdlSet{S, T}(O::AlgAssAbsOrd{S, T}) where {S, T}
    z = new{S, T}(O)
    return z
  end
end

function AlgAssAbsOrdIdlSet(O::AlgAssAbsOrd{S, T}) where {S, T}
  return AlgAssAbsOrdIdlSet{S, T}(O)
end

mutable struct AlgAssAbsOrdFracIdl{S, T}
  order::AlgAssAbsOrd{S, T}
  num::AlgAssAbsOrdIdl{S, T}
  den::fmpz
  basis_mat::FakeFmpqMat
  basis_mat_inv::FakeFmpqMat

  function AlgAssAbsOrdFracIdl{S, T}(O::AlgAssAbsOrd{S, T}, a::AlgAssAbsOrdIdl{S, T}, b::fmpz) where {S, T}
    z = new{S, T}()
    z.order = O
    z.basis_mat = FakeFmpqMat(basis_mat(a), deepcopy(b))
    z.num = a
    z.den = b
    return z
  end

  function AlgAssAbsOrdFracIdl{S, T}(O::AlgAssAbsOrd{S, T}, M::FakeFmpqMat) where {S, T}
    z = new{S, T}()
    z.order = O
    z.basis_mat = M
    return z
  end
end

################################################################################
#
#  AlgAssAbsOrdQuoRing / AlgAssAbsOrdQuoRingElem
#
################################################################################

const AlgAssAbsOrdQuoRing{S, T} = AbsOrdQuoRing{AlgAssAbsOrd{S, T}, AlgAssAbsOrdIdl{S, T}} where {S, T}

const AlgAssAbsOrdQuoRingElem{S, T} = AbsOrdQuoRingElem{AlgAssAbsOrd{S, T}, AlgAssAbsOrdIdl{S, T}, AlgAssAbsOrdElem{S, T}} where {S, T}

################################################################################
#
#  AlgMat / AlgMatElem
#
################################################################################

# T == elem_type(base_ring), S == dense_matrix_type(coefficient_ring)
mutable struct AlgMat{T, S} <: AbsAlgAss{T}
  base_ring::Ring
  coefficient_ring::Ring
  one
  has_one::Bool
  basis
  basis_mat # matrix over the base_ring
  dim::Int
  degree::Int
  issimple::Int
  decomposition
  maximal_order
  mult_table::Array{T, 3} # e_i*e_j = sum_k mult_table[i, j, k]*e_k

  function AlgMat{T, S}(R::Ring) where {T, S}
    A = new{T, S}()
    A.base_ring = R
    A.coefficient_ring = R
    A.issimple = 0
    return A
  end

  function AlgMat{T, S}(R1::Ring, R2::Ring) where {T, S}
    A = new{T, S}()
    A.base_ring = R1
    A.coefficient_ring = R2
    A.issimple = 0
    return A
  end
end

mutable struct AlgMatElem{T, S, Mat} <: AbsAlgAssElem{T}
  parent::S
  matrix::Mat # over the coefficient ring of the parent
  coeffs::Vector{T} # over the base ring of the parent
  has_coeffs::Bool

  function AlgMatElem{T, S, Mat}(A::S) where {T, S, Mat}
    z = new{T, S, Mat}()
    z.parent = A
    z.matrix = zero_matrix(base_ring(A), degree(A), degree(A))
    z.has_coeffs = false
    return z
  end

  function AlgMatElem{T, S, Mat}(A::S, M::Mat) where {T, S, Mat}
    z = new{T, S, Mat}()
    z.parent = A
    z.matrix = M
    z.has_coeffs = false
    return z
  end
end


text:
import Base.show, Base.display

function printnode(io::IO, node::Node, maxdepth::Int = 4, depth::Int = 0)
    junction = ifelse(depth > 0, "├", "┌")
    tab = "   │   " ^ depth * "   $(junction)── "
    if node.isleaf
        if node.index != 0
            @printf(io, "%s(%d)\n", tab, node.index)
        else
            @printf(io, "%s*\n", tab)
        end
    elseif depth == maxdepth
        println(io, tab, "...")
    else
        @printf(io, "%s%d, %#-.1f\n", tab, node.index, node.radius)
        #@printf(io, "%s[%d, %d], %d, %#-.1f\n", tab, node.i, node.j, node.index, sqrt(node.radius))
        printnode(io, node.inside, maxdepth, depth + 1)
        printnode(io, node.outside, maxdepth, depth + 1)
    end
end

show(io::IO, node::Node, args...) = printnode(io, node, args...)
display(node::Node) = show(node)

show(io::IO, vptree::VPTree, args...) = show(io, vptree.root, args...)
display(vptree::VPTree) = show(vptree.root)


text:
# promotion methods for CircuitParams
for p in (:Impedance, :Admittance)
    @eval promote_rule(::Type{($p){T}}, ::Type{($p){S}}) where {T<:Real, S<:Real} =
        ($p){promote_type(T, S)}
end

# convert methods for CircuitParams, Array{CircuitParams, 1}, and CircuitData
convert(::Type{Impedance{T}}, z::Impedance{S}) where {T<:Real, S<:Real} =
     Impedance(convert(Complex{T}, z.data))
convert(::Type{Admittance{T}}, y::Admittance{S}) where {T<:Real, S<:Real} =
     Admittance(convert(Complex{T}, y.data))

convert(::Type{Admittance}, z::Impedance{T}) where {T<:Real} = Admittance(1/z.data)
convert(::Type{Impedance}, y::Admittance{T}) where {T<:Real} = Impedance(1/y.data)

# shortcut convert methods
convert(::Type{T}, D::Vector{S}) where {T<:CircuitParams, S<:CircuitParams} =
    [convert(T, d) for d in D]
convert(::Type{T1}, D::CircuitData{S, T2}) where {S<:Real,
    T1<:CircuitParams, T2<:CircuitParams} =
    CircuitData(D.frequency, convert(T1, D.params))


text:
function testNHSolidBodySphere
clear all
close all

% Physical parameters
Param=PhysParameters();

% Grid
nz=10;
Param.nPanel=4;
Param.H=30000;
Param.Grid=CubedGrid(Param.nPanel,OrientFaceSphere,Param);

Param.Grid.nz=nz;
Param.Grid.zP=zeros(nz,1);
Param.Grid.z=zeros(nz+1,1);
Param.Grid.dz=Param.H/nz;
Param.Grid.zP(1)=Param.Grid.dz/2;
for i=2:nz
  Param.Grid.zP(i)=Param.Grid.zP(i-1)+Param.Grid.dz;
end
for i=2:nz+1
  Param.Grid.z(i)=Param.Grid.z(i-1)+Param.Grid.dz;
end

% Model
Param.Deep=false;
Param.HeightLimit=30000.0;
Param.T0E=310.0;
Param.T0P=240.0;
Param.B=2.0;
Param.K=3.0;
Param.LapseRate=0.005;
Param.U0=-0.5;
Param.PertR=1.0/6.0;
Param.Up=0.0; %1.0;
Param.PertExpR=0.1;
Param.PertLon=pi/9.0;
Param.PertLat=2.0*pi/9.0;
Param.PertZ=15000.0;

Param.StrideDamp=6000;
Param.Relax=1.e-3;
Param.Damping=false;
Param.Coriolis=true;
Param.CoriolisType='Sphere';
Param.Buoyancy=true;
Param.Source=false;
Param.Th0=300;
Param.uMax=10;
Param.vMax=0;
Param.NBr=1.e-2;
Param.DeltaT=1;
Param.ExpDist=5;
Param.T0=300;
Param.Equation='Compressible';
Param.TopoS='';
Param.lat0=0;
Param.lon0=pi/2;
Param.ProfVel='SolidBody';
Param.ProfRho='SolidBody';
Param.ProfTheta='SolidBody';
Param.NumV=5;
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
Param.Thermo='';%'Energy';

% Output
Param.Flat=true;
Param.level=3;
Param.fig=1;
Param.SliceXY.Type='XY';
Param.SliceXY.iz=1;

% Discretization

OrdPoly=4;
[CG,Param]=Discretization(OrdPoly,@JacobiSphere3,Param);
LRef=11*1.e5;
dx=2*pi*Param.RadEarth/4/Param.nPanel/OrdPoly;
Param.HyperD=2.e17; %1.e14*(dx/LRef)^3.2;

% Initial conditions

U=zeros(CG.NumG,nz,Param.NumV);
U(:,:,Param.RhoPos)=Project(@fRho,CG,Param);
[U(:,:,Param.uPos),U(:,:,Param.vPos)]=ProjectVec(@fVel,CG,Param);
U(:,:,Param.ThPos)=Project(@fTheta,CG,Param).*U(:,:,Param.RhoPos);
PresStart=Pressure(U(:,:,Param.ThPos),U(:,:,Param.ThPos),U(:,:,Param.ThPos),Param);
Param.fig=PlotCG(U(:,:,Param.ThPos)./U(:,:,Param.RhoPos)...
  ,CG,@TransSphere,@Topo,Param,Param.fig,Param.SliceXY);
if strcmp(Param.Thermo,'Energy')
  U(:,:,Param.ThPos)=PotToEnergy(U,CG,Param);
end

% Integration
CFL=0.125;
dtau=500;
time=0;

IntMethod='Rosenbrock'; %'RungeKutta';
%IntMethod='RungeKutta';
nIter=20000;
%v = VideoWriter ('Galewsky.avi');
%open (v);
Param.RK=RungeKuttaMethod('RK4');
Param.ROS=RosenbrockMethod('ROSRK3');
SimDays=100;
PrintDay=5;
nIter=round(24*3600*SimDays/dtau);
PrintInt=round(24*3600*PrintDay/dtau);
switch IntMethod
  case 'Rosenbrock'
    for i=1:nIter
      i
      sum(sum(U(:,:,Param.wPos)))
      U=Rosenbrock(U,dtau,@FcnNHCurlVec,@Jac,CG,Param);
      time=time+dtau;
      if mod(i,4000)==0
        Pres=Pressure(U(:,:,Param.ThPos),U(:,:,Param.ThPos),U(:,:,Param.ThPos),Param);
        if strcmp(Param.Thermo,'Energy')
          Th=EnergyToPot(U,CG,Param);
        else
          Th=U(:,:,Param.ThPos);
        end
        Param.fig=PlotCG(Th./U(:,:,Param.RhoPos)...
          ,CG,@TransSphere,@Topo,Param,Param.fig,Param.SliceXY);
        Param.fig=PlotCG(U(:,:,Param.uPos),CG,@TransSphere,@Topo,Param,Param.fig,Param.SliceXY);
        Param.fig=PlotCG(U(:,:,Param.vPos),CG,@TransSphere,@Topo,Param,Param.fig,Param.SliceXY);
        W=zeros(size(U,1),nz+1);
        W(:,2:nz+1)=U(:,:,Param.wPos);
        Param.fig=PlotCG(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,@TransSphere...
          ,@Topo,Param,Param.fig,Param.SliceXY);
      end
    end
  case 'RungeKutta'
    for i=1:nIter
      U=RungeKuttaExplicit(U,dtau,@FcnNHCurlVec,CG,Param);

      time=time+dtau;
      if mod(i,1)==0
        Param.fig=PlotCG(U(:,:,Param.ThPos)./U(:,:,Param.RhoPos)...
          ,CG,@TransSphere,@Topo,Param,Param.fig,Param.SliceXY);
        Param.fig=PlotCG(U(:,:,Param.uPos),CG,@TransSphere,@Topo,Param,Param.fig,Param.SliceXZ);
        Param.fig=PlotCG(U(:,:,Param.vPos),CG,@TransSphere,@Topo,Param,Param.fig,Param.SliceXZ);
        W=zeros(size(U,1),nz+1);
        W(:,2:nz+1)=U(:,:,Param.wPos);
        Param.fig=PlotCG(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,@TransSpheres...
          ,@Topo,Param,Param.fig,Param.SliceXZ);
        frame=getframe(gcf);
        writeVideo (v, frame);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
%close(v)
end


text:
@testset "Co-ordinate system conversion" begin
    @testset "Fixed conversions" begin
        ###################################
        ### Testing fixed relationships ###
        ###################################

        lla = LLA(42.3673, -71.0960, 0.)
        lla_ref = LLA(42.36299, -71.09183, 0.)

        # LLA <-> ECEF
        ecef = ECEF(1529073.1560519305, -4465040.019013103, 4275835.339260307)
        @test ECEF(lla, wgs84) ≈ ecef
        @test LLA(ecef, wgs84) ≈ lla

        # LLA <-> ENU
        enu = ENU(-343.493749083977, 478.764855466788, -0.027242885224325164)
        @test ENU(lla, lla_ref, wgs84) ≈ enu
        @test LLA(enu, lla_ref, wgs84) ≈ lla

        # ECEF <-> ENU
        @test ENU(ecef, lla_ref, wgs84) ≈ enu
        @test ECEF(enu, lla_ref, wgs84) ≈ ecef

        # ENU <-> NED
        ned = NED(478.764855466788, -343.493749083977, 0.027242885224325164)
        @test NED(enu) == ned
        @test ENU(ned) == enu
        @test NED((ENU(ned))) == ned
        @test ENU((NED(enu))) == enu

         # LLA <-> NED
         ned = NED(enu.n, enu.e, -enu.u)
         @test NED(lla, lla_ref, wgs84) ≈ ned
         @test LLA(ned, lla_ref, wgs84) ≈ lla

         # ECEF <-> NED
         @test NED(ecef, lla_ref, wgs84) ≈ ned
         @test ECEF(ned, lla_ref, wgs84) ≈ ecef

        # LLA <-> UTM
        (z, h) = (19, true)
        utm = UTM(327412.48528248386, 4.692686244318043e6, 0.0)
        @test UTM(lla, z, h, wgs84) ≈ utm
        @test LLA(utm, z, h, wgs84) ≈ lla

        # ECEF <-> UTM
        @test UTM(ecef, z, h, wgs84) ≈ utm
        @test ECEF(utm, z, h, wgs84) ≈ ecef

        # ENU <-> UTM
        @test UTM(enu, z, h, lla_ref, wgs84) ≈ utm
        @test ENU(utm, z, h, lla_ref, wgs84) ≈ enu

        # NED <-> UTM
        @test UTM(ned, z, h, lla_ref, wgs84) ≈ utm
        @test NED(utm, z, h, lla_ref, wgs84) ≈ ned

        # LLA <-> UTMZ
        utmz = UTMZ(327412.48528248386, 4.692686244318043e6, 0, z, h)
        @test UTMZ(lla, wgs84) ≈ utmz
        @test LLA(utmz, wgs84) ≈ lla

        # ECEF <-> UTMZ
        @test UTMZ(ecef, wgs84) ≈ utmz
        @test ECEF(utmz, wgs84) ≈ ecef

        # ENU <-> UTMZ
        @test UTMZ(enu, lla_ref, wgs84) ≈ utmz
        @test ENU(utmz, lla_ref, wgs84) ≈ enu

        # NED <-> UTMZ
        @test UTMZ(ned, lla_ref, wgs84) ≈ utmz
        @test NED(utmz, lla_ref, wgs84) ≈ ned

        # UTM <-> UTMZ
        @test UTMZ(utm, z, h, wgs84) ≈ utmz
        @test UTM(utmz, z, h, wgs84) ≈ utm

        # Identity
        @test ECEF(ecef, wgs84) == ecef
        @test LLA(lla, wgs84) == lla
        @test ENU(enu, wgs84) == enu
        @test NED(ned, wgs84) == ned
        @test UTM(utm, wgs84) == utm
        @test UTMZ(utmz, wgs84) == utmz
    end

    @testset "Random distance tests" begin
        #############################
        ### Testing random errors ###
        #############################
        """
            randLLA(h_min = -5e6, h_max = 5e6)

        Random LLA coordinate with evenly-distributed angles (probabilty density
        uniform over the surface-area of a sphere)
        """
        function randLLA(h_min = -5e6, h_max = 5e6)
            lon = 180*(2*rand()-1)
            lat = asind(2*rand()-1)
            h = h_min + (h_max-h_min)*rand()
            return LLA(lat,lon,h)
        end

        number_of_utm_distance_tests = 0
        for _ = 1:5000
            lla1 = randLLA()
            lla2 = randLLA()

            d = euclidean_distance(lla1, lla2)

            ecef1 = ECEF(lla1, wgs84)
            ecef2 = ECEF(lla2, wgs84)

            @test euclidean_distance(ecef1, ecef2) ≈ sqrt((ecef1.x - ecef2.x)^2 + (ecef1.y - ecef2.y)^2 + (ecef1.z - ecef2.z)^2)
            @test euclidean_distance(ecef1, ecef2) ≈ d

            utmz1 = UTMZ(lla1, wgs84)
            utmz2 = UTMZ(lla2, wgs84)

            @test euclidean_distance(utmz1, utmz2) ≈ d

            if utmz1.zone == utmz2.zone && utmz1.isnorth == utmz2.isnorth
                number_of_utm_distance_tests += 1

                utm1 = UTM(utmz1)
                utm2 = UTM(utmz2)

                zone = utmz1.zone
                isnorth = utmz1.isnorth

                @test euclidean_distance(utm1,  utm2,  zone, isnorth) ≈ d
                @test euclidean_distance(utm1,  ecef2, zone, isnorth) ≈ d
                @test euclidean_distance(ecef1, utm2,  zone, isnorth) ≈ d
            end

            enu000 = ENU(0.0, 0.0, 0.0)
            enu = ENU(ecef1, ecef2, wgs84)

            @test euclidean_distance(enu, enu000) ≈ d

            ned000 = NED(0.0, 0.0, 0.0)
            ned = NED(ecef1, ecef2, wgs84)

            @test euclidean_distance(ned, ned000) ≈ d
        end
        @test number_of_utm_distance_tests > 0
    end
end # @testset


text:
export stack, unstack, create_dummies

"""
`stack(t, by = pkeynames(t); select = excludecols(t, by), variable = :variable, value = :value)`

Reshape a table from the wide to the long format. Columns in `by` are kept as indexing columns.
Columns in `select` are stacked. In addition to the id columns, two additional columns labeled `variable` and `value`
are added, containg the column identifier and the stacked columns.

## Examples

```jldoctest stack
julia> t = table(1:4, [1, 4, 9, 16], [1, 8, 27, 64], names = [:x, :xsquare, :xcube], pkey = :x);

julia> stack(t)
Table with 8 rows, 3 columns:
x  variable  value
──────────────────
1  :xsquare  1
1  :xcube    1
2  :xsquare  4
2  :xcube    8
3  :xsquare  9
3  :xcube    27
4  :xsquare  16
4  :xcube    64
```
"""
function stack(t::D, by = pkeynames(t); select = isa(t, NDSparse) ? valuenames(t) : excludecols(t, by),
    variable = :variable, value = :value) where {D<:Dataset}

    (by != pkeynames(t)) && return stack(reindex(t, by, select); variable = :variable, value = :value)

    valuecols = columns(t, select)
    valuecol = [valuecol[i] for i in 1:length(t) for valuecol in valuecols]

    labels = fieldnames(valuecols)
    labelcol = [label for i in 1:length(t) for label in labels]

    bycols = map(arg -> repeat(arg, inner = length(valuecols)), columns(t, by))
    convert(collectiontype(D), Columns(bycols), Columns(labelcol, valuecol, names = [variable, value]))
end

function unstack(::Type{D}, ::Type{T}, key, val, cols::AbstractVector{S}) where {D <:Dataset, T, S}
    dest_val = Columns((DataValues.DataValueArray{T}(length(val)) for i in cols)...; names = cols)
    for (i, el) in enumerate(val)
        for j in el
            k, v = j
            isnull(columns(dest_val, S(k))[i]) || error("Repeated values with same label are not allowed")
            columns(dest_val, S(k))[i] = v
        end
    end
    convert(collectiontype(D), key, dest_val)
end

"""
`unstack(t, by = pkeynames(t); variable = :variable, value = :value)`

Reshape a table from the long to the wide format. Columns in `by` are kept as indexing columns.
Keyword arguments `variable` and `value` denote which column contains the column identifier and
which the corresponding values.

## Examples

```jldoctest unstack
julia> t = table(1:4, [1, 4, 9, 16], [1, 8, 27, 64], names = [:x, :xsquare, :xcube], pkey = :x);

julia> long = stack(t)
Table with 8 rows, 3 columns:
x  variable  value
──────────────────
1  :xsquare  1
1  :xcube    1
2  :xsquare  4
2  :xcube    8
3  :xsquare  9
3  :xcube    27
4  :xsquare  16
4  :xcube    64

julia> unstack(long)
Table with 4 rows, 3 columns:
x  xsquare  xcube
─────────────────
1  1        1
2  4        8
3  9        27
4  16       64
```
"""
function unstack(t::D, by = pkeynames(t); variable = :variable, value = :value) where {D<:Dataset}
    tgrp = groupby((value => identity,), t, by, select = (variable, value))
    S = eltype(colnames(t))
    cols = S.(union(columns(t, variable)))
    T = eltype(columns(t, value))
    unstack(D, T isa Type{<:DataValue} ? eltype(T) : T, pkeys(tgrp), columns(tgrp, value), cols)
end

"""
`create_dummies(t, cols; categories = Dict())`

Create a column corresponding to every categorical value present in the columns specified. `cols` is supposed be
an array of column names or a single column name (`Symbol`) on which the dummy variables are to be created.
The columns would have a 1 for the row with the value corresponding with the `val` as specified in the column name and 0 otherwise.
`categories` is a dictionary from column name to a vector of specific unique values.

## Examples

```jldoctest create_dummies
julia> t = table([1, 4, 9, 16], ["val1", "val1", "val2", "val2"], names = [:x, :y], pkey = :x);

julia> create_dummies(t, :y)
Table with 4 rows, 3 columns:
x   y_val1   y_val2
──────────────────
1   1        0
4   1        0
9   0        1
16  0        1

julia> t = table([1, 4, 9, 16], ["val1", "val1", "val2", "val2"], names = [:x, :y], pkey = :x);

julia> create_dummies(t, :y, categories = Dict(:y => ["val1"]))
Table with 4 rows, 3 columns:
x   y_val1
───────────
1      1
4      1
9      0
16     0
```
"""
function create_dummies(t, cols; categories=Dict())
    arr = []
    for i in cols
        if haskey(categories, i)
            uniq = categories[i]
        else
            uniq = collect(keys(reduce(CountMap(), t, select = i).value))
        end

        for j in uniq
            push!(arr, Symbol("$(i)_$(j)") => i => f(j))
        end
    end
    t = setcol(t, arr)
    t = popcol(t, cols)
    t
end
create_dummies(t, col::Symbol; categories = Dict()) = create_dummies(t, [col], categories = categories)
f(j) = x -> x == j ? Int8(1) : Int8(0)


text:
@ghdef mutable struct CheckRun
    name::Union{String,Nothing}
    details_url::Union{String,Nothing}
    external_id::Union{String,Nothing}
    started_at::Union{DateTime, Nothing}
    status::Union{String,Nothing}
    conclusion::Union{String, Nothing}
    completed_at::Union{DateTime, Nothing}
    output::Union{Output, Nothing}
    actions::Union{Vector{Action}, Nothing}

    # output
    id::Union{Int, Nothing}
    head_sha::Union{String, Nothing}
    app::Union{App, Nothing}
    pull_requests::Union{Vector{PullRequest}, Nothing}
    check_suite::Union{CheckSuite, Nothing}
    html_url::Union{HTTP.URI, Nothing}
end
namefield(cr::CheckRun) = cr.id

"""
    GitHub.check_runs([api,] repo::Repo, suite::CheckSuite; options...)

# Parameters
 - `check_name::String`: Returns check runs with the specified `name`.
 - `status::String`: Returns check runs with the specified `status`. Can be one of `"queued"`, `"in_progress"`, or `"completed"`.
 - `filter::String`: Filters check runs by their `completed_at` timestamp. Can be one of `"latest"` (returning the most recent check runs) or `"all"`. Default: `"latest"`

# External links
- https://developer.github.com/v3/checks/runs/#list-check-runs-in-a-check-suite
"""
@api_default function check_runs(api::GitHubAPI, repo::Repo, suite::CheckSuite; headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/check-suites/$(name(suite))/check-runs";
        headers=headers, options...)
    map(CheckRun, results["check_runs"]), page_data, results["total_count"]
end

"""
    GitHub.check_runs([api,] repo::Repo, ref; options...)

List the `CheckRun`s for `ref`

# Parameters
 - `check_name::String`: Returns check runs with the specified `name`.
 - `status::String`: Returns check runs with the specified `status`. Can be one of `"queued"`, `"in_progress"`, or `"completed"`.
 - `filter::String`: Filters check runs by their `completed_at` timestamp. Can be one of `"latest"` (returning the most recent check runs) or `"all"`. Default: `"latest"`

# External links
- https://developer.github.com/v3/checks/runs/#list-check-runs-for-a-specific-ref
"""
@api_default function check_runs(api::GitHubAPI, repo::Repo, ref; headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    results, page_data, total_count = gh_get_paged_json(api, "/repos/$(name(repo))/commits/$(ref)/check-runs", "check_runs";
                                           headers=headers, options...)
    map(CheckRun, results), page_data, total_count
end

"""
    GitHub.create_check_run([api,] repo::Repo; options...)

- https://developer.github.com/v3/checks/runs/#create-a-check-run
"""
@api_default function create_check_run(api::GitHubAPI, repo::Repo; headers = Dict(), kwargs...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    result = gh_post_json(api, "/repos/$(name(repo))/check-runs"; headers=headers, kwargs...)
    return CheckRun(result)
end

"""
    GitHub.update_check_run([api,] repo::Repo, cr::CheckRun; options...)

- https://developer.github.com/v3/checks/runs/#update-a-check-run
"""
@api_default function update_check_run(api::GitHubAPI, repo::Repo, cr::CheckRun; kwargs...)
    update_check_run(api, repo, cr.id; kwargs...)
end
@api_default function update_check_run(api::GitHubAPI, repo::Repo, id::Int; headers = Dict(), kwargs...)
    # TODO: deprecate?
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    result = gh_patch_json(api, "/repos/$(name(repo))/check-runs/$(id)"; headers=headers, kwargs...)
    return CheckRun(result)
end


text:
import BlockArrays: BlockIndex, BlockIndexRange, globalrange, nblocks, global2blockindex, blockindex2global

@testset "Blocks" begin
    @test Int(Block(2)) === Integer(Block(2)) === Number(Block(2)) === 2
    @test Block((Block(3), Block(4))) === Block(3,4)
end

#=
[1,1  1,2] | [1,3  1,4  1,5]
--------------------------
[2,1  2,2] | [2,3  2,4  2,5]
[3,1  3,2] | [3,3  3,4  3,5]
----------------------------
[4,1  4,2] | [4,3  4,4  4,5]
[5,1  5,2] | [5,3  5,4  5,5]
[6,1  6,2] | [6,3  6,4  6,5]
=#

@testset "BlockSizes / BlockIndices" begin
    block_size = BlockArrays.BlockSizes([1,2,3], [2, 3])

    @test nblocks(block_size) == (3,2)
    @test nblocks(block_size, 1) == 3
    @test nblocks(block_size, 2) == 2

    @test @inferred(globalrange(block_size, (1,1))) == (1:1, 1:2)
    @test @inferred(globalrange(block_size, (1,2))) == (1:1, 3:5)
    @test @inferred(globalrange(block_size, (2,1))) == (2:3, 1:2)
    @test @inferred(globalrange(block_size, (2,2))) == (2:3, 3:5)

    # Test for allocations inside a function to avoid noise due to global
    # variable references
    wrapped_allocations = (bs, i) -> @allocated(globalrange(bs, i))
    @test wrapped_allocations(block_size, (1, 1)) == 0

    @test @inferred(global2blockindex(block_size, (3, 1))) == BlockIndex((2,1), (2,1))
    @test @inferred(global2blockindex(block_size, (1, 4))) == BlockIndex((1,2), (1,2))
    @test @inferred(global2blockindex(block_size, (4, 5))) == BlockIndex((3,2), (1,3))

    wrapped_allocations = (bs, i) -> @allocated(global2blockindex(bs, i))
    @test wrapped_allocations(block_size, (3, 1)) == 0

    @test @inferred(blockindex2global(block_size, BlockIndex((2,1), (2,1)))) == (3, 1)
    @test @inferred(blockindex2global(block_size, BlockIndex((1,2), (1,2)))) == (1, 4)
    @test @inferred(blockindex2global(block_size, BlockIndex((3,2), (1,3)))) == (4, 5)

    wrapped_allocations = (bs, i) -> @allocated(blockindex2global(bs, i))
    @test wrapped_allocations(block_size, BlockIndex((2,1), (2,1))) == 0

    @test block_size == BlockArrays.BlockSizes(1:3, 2:3)

    buf = IOBuffer()
    print(buf, block_size)
    @test String(take!(buf)) == "[1, 2, 3] × [2, 3]"

    @test BlockArrays.searchlinear([1,2,3], 5) == 3
    @test Base.dataids(block_size) == Base.dataids(block_size)

    @test Block(1)[1] == BlockIndex((1,),(1,))
    @test Block(1)[1:2] == BlockIndexRange(Block(1),(1:2,))
    @test Block(1,1)[1,1] == BlockIndex((1,1),(1,1))
    @test Block(1,1)[1:2,1:2] == BlockIndexRange(Block(1,1),(1:2,1:2))

    A = BlockVector([1,2,3],[1,2])
    @test A[Block(2)[2]] == 3
    @test A[Block(2)[1:2]] == [2,3]
    @test A[getindex.(Block.(1:2), 1)] == [1,2]

    @test_throws BlockBoundsError A[Block(3)]
    @test_throws BlockBoundsError A[Block(3)[1]]
    @test_throws BoundsError A[Block(3)[1:1]] # this is likely an error
    @test_throws BoundsError A[Block(2)[3]]
    @test_throws BoundsError A[Block(2)[3:3]]
end

@testset "sortedin" begin
    v = [1,3,4]
    @test BlockArrays.sortedin(1,v)
    @test !BlockArrays.sortedin(2,v)
    @test !BlockArrays.sortedin(0,v)
    @test !BlockArrays.sortedin(5,v)
end


text:
using DiffEqBayesStan, OrdinaryDiffEq, ParameterizedFunctions,
      RecursiveArrayTools, Distributions, MCMCChains, Test

println("\n\nFour parameter case\n")

f1 = @ode_def begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d

u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())

t = collect(range(1,stop=10,length=50))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(1.5,0.01),0.5,2),truncated(Normal(1.0,0.01),0.5,1.5),
          truncated(Normal(3.0,0.01),0.5,4),truncated(Normal(1.0,0.01),0.5,2)]

bayesian_result = stan_inference(prob1,t,data,priors;
    num_chains=1,num_samples=1000,num_warmups=1000,
    vars =(DiffEqBayesStan.StanODEData(),InverseGamma(4,1)))

@test mean(get(bayesian_result.chains,:theta_1)[1]) ≈ 1.5 atol=1e-1
@test mean(get(bayesian_result.chains,:theta_2)[1]) ≈ 1.0 atol=1e-1
@test mean(get(bayesian_result.chains,:theta_3)[1]) ≈ 3.0 atol=1e-1
@test mean(get(bayesian_result.chains,:theta_4)[1]) ≈ 1.0 atol=1e-1


text:
using Random, Test
using StatsBase
using Distributions
using DependentBootstrap
using DataFrames, TimeSeries

#Quick and dirty function for simulating AR(1) data with deterministic start point for random number generation
function temp_ar(seedint::Int)
    #srand(seedInt)
    Random.seed!(seedint)
    N = 100
    e = randn(N)
    x = NaN*ones(Float64,N)
    x[1] = 0.0
    for n = 2:N
        x[n] = 0.8 * x[n-1] + e[n]
    end
    return x
end

#Get the AR(1) data
x = temp_ar(1234);
xmat = hcat(x, temp_ar(5678));
xvv = [ xmat[:, k] for k = 1:size(xmat,2) ]

#bootstrap methods
bootmethodvec = [:iid, :stationary, :moving, :circular, :nooverlap];
bootmethodtypevec = [DependentBootstrap.BootIID(), DependentBootstrap.BootStationary(), DependentBootstrap.BootMoving(),
                     DependentBootstrap.BootCircular(), DependentBootstrap.BootNoOverlap()]
blocklengthvec = [0.0, 5.0];
blocklengthmethodvec = [:ppw2009]
blocklengthmethodtypevec = [DependentBootstrap.BLPPW2009()]

#Test constructor
@testset "BootInput constructor test" begin
    for kbm = 1:length(bootmethodvec)
        for kbl = 1:length(blocklengthvec)
            for kblm = 1:length(blocklengthmethodvec)
                bi = BootInput(x, 100, blocklengthvec[kbl], 200, bootmethodtypevec[kbm],
                               blocklengthmethodtypevec[kblm], var, std, 300, mean)
                @test bi.numresample == 200
                @test bi.bootmethod == bootmethodtypevec[kbm]
                @test bi.blocklengthmethod == blocklengthmethodtypevec[kblm]
                @test bi.flevel1 == var
                @test bi.flevel2 == std
                @test bi.numobsperresample == 300
                @test bi.fblocklengthcombine == mean
                bi = BootInput(x, numresample=200, bootmethod=bootmethodvec[kbm], blocklength=blocklengthvec[kbl],
                               blocklengthmethod=blocklengthmethodvec[kblm], flevel1=var, flevel2=std,
                               numobsperresample=300, fblocklengthcombine=mean)
                @test bi.numresample == 200
                @test bi.bootmethod == bootmethodtypevec[kbm]
                @test bi.blocklengthmethod == blocklengthmethodtypevec[kblm]
                @test bi.flevel1 == var
                @test bi.flevel2 == std
                @test bi.numobsperresample == 300
                @test bi.fblocklengthcombine == mean
            end
        end
    end
end

bootmethodvec = Symbol[:iid, :stationary, :moving]
correctblocklength1 = Vector{Float64}[Float64[1.0, 7.002404488495543, 8.015752150100226]];
correctblocklength2 = Vector{Float64}[Float64[1.0, 6.837001835299588, 7.826413377230708]];

#Test block length selection procedure
@testset "Block length selection test" begin
    for kbm = 1:length(bootmethodvec)
        for kblm = 1:length(blocklengthmethodvec)
            bi = BootInput(xmat[:, 1], numresample=200, bootmethod=bootmethodvec[kbm], blocklengthmethod=blocklengthmethodvec[kblm])
            @test isapprox(bi.blocklength, correctblocklength1[kblm][kbm])
            bi = BootInput(xmat[:, 2], numresample=200, bootmethod=bootmethodvec[kbm], blocklengthmethod=blocklengthmethodvec[kblm])
            @test isapprox(bi.blocklength, correctblocklength2[kblm][kbm])
            bi = BootInput(xmat, numresample=200, bootmethod=bootmethodvec[kbm], blocklengthmethod=blocklengthmethodvec[kblm], fblocklengthcombine=mean)
            @test isapprox(bi.blocklength, mean([correctblocklength1[kblm][kbm],correctblocklength2[kblm][kbm]]))
        end
    end
end

#Test univariate bootstrap method
bootmethodvec = [:iid, :stationary, :moving, :circular, :nooverlap];
correctbootunivariate = Float64[0.02231785457420185,0.0805280736777265,0.07755558515710691,0.07248357346968781,0.0718412998084052]
correctbootunivariatebl1 = 0.021374650187840242*ones(Float64, length(bootmethodvec))
@testset "Univariate bootstrap test" begin
    for kbm = 1:length(bootmethodvec)
        Random.seed!(1234)
        y = dboot(x, numresample=1000, bootmethod=bootmethodvec[kbm], blocklength=5, flevel1=mean, flevel2=var)
        @test isapprox(y, correctbootunivariate[kbm])
    end
    for kbm = 1:length(bootmethodvec)
        Random.seed!(1234)
        y = dboot(x, numresample=500, bootmethod=bootmethodvec[kbm], blocklength=1, flevel1=mean, flevel2=var)
        @test isapprox(y, correctbootunivariatebl1[kbm])
    end
end

#Test multivariate bootstrap
bootmethodvec = [:iid, :stationary, :moving, :circular, :nooverlap];
correctbootmultmatrix = Float64[0.19194008990823777,0.40101052326568176,0.5186276897440075,0.44038974430819566,0.3960765922525029]
correctbootmultvv = Float64[0.1461957546031049,0.3546585480604615,0.316914521186327,0.34363235435327766,0.361937416083102]
@testset "Multivariate bootstrap test" begin
    for kbm = 1:length(bootmethodvec)
        Random.seed!(1234)
        y = dbootvar(xmat, numresample=1000, bootmethod=bootmethodvec[kbm], blocklength=5, flevel1=x->minimum(x[:,1].*x[:,2]))
        @test isapprox(y, correctbootmultmatrix[kbm])
        Random.seed!(1234)
        y = dbootconf(xvv, numresample=1000, bootmethod=bootmethodvec[kbm], blocklength=5, flevel1=flevel1=x->mean([mean(x[1]),mean(x[2])]))
        @test isapprox(y[2], correctbootmultvv[kbm])
    end
end

#Test dbootdata and dbootinds
@testset "Test dbootdata and dbootinds" begin
    a = dbootinds(x, numresample=20, blocklength=5, bootmethod=:stationary)
    @test length(a) == 20
    a = dbootdata(x, numresample=20, blocklength=5, bootmethod=:stationary)
    @test length(a) == 20
end

#Test exotic dataset types
@testset "Exotic data types test" begin
    Random.seed!(1234)
    xdf = DataFrame(xmat)
    #y = dboot(xdf, bootmethod=:stationary, blocklength=5, numresample=1000, flevel1=x->mean(DataFrames.columns(x)[1]))
    y = dboot(xdf, bootmethod=:stationary, blocklength=5, numresample=1000, flevel1=(x->mean(x[:,1])))
    @test isapprox(y, 0.0805280736777265)
    dtvec = [ Date(2000)+Day(n) for n = 1:size(xmat,1) ]
    xta1 = TimeSeries.TimeArray(dtvec, x)
    xta2 = TimeSeries.TimeArray(dtvec, xmat)
    Random.seed!(1234)
    y1 = dboot(xta1, bootmethod=:stationary, blocklength=5, numresample=1000, flevel1=x->mean(values(x)))
    Random.seed!(1234)
    y2 = dboot(xta2, bootmethod=:stationary, blocklength=5, numresample=1000, flevel1=x->mean(values(x)))
    @test isapprox(y1, 0.0805280736777265)
    @test isapprox(y2, 0.05077003035163576)
end


text:
module LoadQML

export load_qml_app

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
  error("$depsfile not found, package QML did not build properly")
end
include(depsfile)

load_qml_app(filename::AbstractString) = ccall((:load_qml_app, _l_loadqml), Void, (Cstring,), filename)

end # module


text:
using Photometry.Aperture: circular_overlap_core,
                           circular_overlap_single_exact,
                           circular_overlap_single_subpixel,
                           area_arc,
                           area_triangle,
                           inside_ellipse,
                           inside_triangle,
                           inside_rectangle,
                           circle_line,
                           circle_segment,
                           circle_segment_single2,
                           triangle_unitcircle_overlap,
                           elliptical_overlap_exact,
                           elliptical_overlap_single_subpixel,
                           rectangular_overlap_exact,
                           rectangular_overlap_single_subpixel

@testset "overlap - circular" begin

    @testset "circular overlap" for r in 0:10:100, use_subpixel in [true, false]

        ap = CircularAperture(0, 0, r)
        if use_subpixel
            ap = Subpixel(ap, 10)
        end
        @test all(p -> 0 ≤ p ≤ 1, ap)
    end

    @testset "overlap core (r = $r)" for r in 1:10
        @test circular_overlap_core(0, 0, 0, 0, r) ≈ 0
        @test circular_overlap_single_exact(0, 0, r, r, r) ≈ r^2 * π / 4
        @test circular_overlap_core(0, 0, r, r, r) ≈ r^2 * π / 4
    end

    @testset "overlap single exact" begin
        @test circular_overlap_single_exact(0, 0, 0, 0, 1) ≈ 0
        @test circular_overlap_single_exact(0, 0, 20, 20, 15) ≈ 176.714586764429
    end

    @testset "overlap subpixel" begin
        @test circular_overlap_single_subpixel(0, 0, 0, 0, 1, 1) ≈ 1

        @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 2) ≈ 0.25
        @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 5) ≈ 0.44
        @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 10) ≈ 0.43
        @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 100) ≈ 0.4423
    end

    @testset "type stability" begin
        r = 5
        @inferred circular_overlap_core(0, 0, r, r, r)
        @inferred circular_overlap_single_exact(0, 0, r, r, r)
        @inferred circular_overlap_single_subpixel(0, 0, 20, 20, 15, 5)
    end

end # circles

@testset "overlap - elliptical" begin

    @testset "elliptical overlap" for r in 0:10:100, use_subpixel in [true, false]
        ap = EllipticalAperture(0, 0, r, r)
        if use_subpixel
            ap = Subpixel(ap, 10)
        end
        @test all(p -> 0 ≤ p ≤ 1, ap)
    end

    @testset "overlap subpixel (elliptical apperture)" begin
        @test elliptical_overlap_single_subpixel(0, 0, 0, 0, 1, 1, 0, 1) ≈ 1

        @test elliptical_overlap_single_subpixel(0, 0, 20, 20, 1 / 225, 1 / 225, 0, 2) ≈ 0.25
        @test elliptical_overlap_single_subpixel(0, 0, 20, 20, 1 / 225, 1 / 225, 0, 5) ≈ 0.44
        @test elliptical_overlap_single_subpixel(0, 0, 20, 20, 1 / 225, 1 / 225, 0, 10) ≈ 0.43
        @test elliptical_overlap_single_subpixel(0, 0, 20, 20, 1 / 225, 1 / 225, 0, 100) ≈ 0.4423
    end

    @testset "circle line" begin
        for line in readlines(joinpath(@__DIR__, "..", "data", "circle_line.csv"))
            tokens = split(line, ',')
            numbers = parse.(Float64, tokens)
            x1, y1, x2, y2 = numbers[1:4]
            point1, point2 = circle_line(x1, y1, x2, y2)
            @test point1[1] ≈ numbers[5] atol = 1e-6
            @test point1[2] ≈ numbers[6] atol = 1e-6
            @test point2[1] ≈ numbers[7] atol = 1e-6
            @test point2[2] ≈ numbers[8] atol = 1e-6
        end
    end

    @testset "circle segment" begin
        for line in readlines(joinpath(@__DIR__, "..", "data", "circle_segment.csv"))
            tokens = split(line, ',')
            numbers = parse.(Float64, tokens)
            x1, y1, x2, y2 = numbers[1:4]
            point1, point2 = circle_segment(x1, y1, x2, y2)
            @test point1[1] ≈ numbers[5] atol = 1e-6
            @test point1[2] ≈ numbers[6] atol = 1e-6
            @test point2[1] ≈ numbers[7] atol = 1e-6
            @test point2[2] ≈ numbers[8] atol = 1e-6
        end
    end

    @testset "circle segment single" begin
        for line in readlines(joinpath(@__DIR__, "..", "data", "circle_segment_single.csv"))
            tokens = split(line, ',')
            numbers = parse.(Float64, tokens)
            x1, y1, x2, y2 = numbers[1:4]
            point1 = circle_segment_single2(x1, y1, x2, y2)
            @test point1[1] ≈ numbers[5] atol = 1e-6
            @test point1[2] ≈ numbers[6] atol = 1e-6
        end
    end

    @testset "triangle unitcircle overlap" begin
        for line in readlines(joinpath(@__DIR__, "..", "data", "triangle_unitcircle_overlap.csv"))
            tokens = split(line, ',')
            numbers = parse.(Float64, tokens)
            x1, y1, x2, y2, x3, y3 = numbers[1:6]
            @test triangle_unitcircle_overlap(x1, y1, x2, y2, x3, y3) ≈ numbers[end] atol = 1e-6
        end
    end
    @testset "elliptical overlap" begin
        @test elliptical_overlap_exact(0.5, 2.5, 1.5, 3.5, 3.0, 3.0, 0) ≈ 0.311725 atol = 1e-6
        @test elliptical_overlap_exact(0, 2, 1, 3, 3.0, 3.0, 0) ≈ 0.943480 atol = 1e-6
    end
end # overlap elliptical

@testset "overlap - rectangular" begin
    @testset "exact" begin
        @test rectangular_overlap_exact(0, 0, 1, 1, 1, 1, 0) ≈ 0.25
    end

    @testset "type stability" begin
        # @inferred rectangular_overlap_exact(0, 0, 1, 1, 2, 2, 0) # TODO area in LazySets.jl not type stable
        @inferred rectangular_overlap_single_subpixel(0, 0, 1, 1, 2, 2, 0, 5)
    end
end # overlap rectangular

@testset "overlap - utils" begin

    @testset "inside triangle" begin
        t1 = [0, 0, 0, 1, 1, 0]
        @test inside_triangle(0, 0, t1...)
        @test inside_triangle(0.1, 0.1, t1...)
        @test inside_triangle(0.25, 0.25, t1...)
    end

    @testset "area triangle" begin
        @test area_triangle(0, 0, 0, 0, 0, 0) ≈ 0
        @test area_triangle(0, 0, 1, 0, 0, 1) ≈ 0.5
        @test area_triangle(1, 1, 1, 4, 6, 1) ≈ 7.5

    end

    @testset "area arc (r = $r)" for r in 1:10
        @test area_arc(0, 0, 0, 0, r) ≈ 0
        @test area_arc(0, r, r, 0, r) ≈ r^2 / 2 * (π / 2 - 1)
    end

    @testset "inside ellipse" begin
        @test inside_ellipse(0, 0, 0, 0, 1, 1, 1)
        @test !inside_ellipse(10, 10, 5, 5, 1, 1, 1)
        @test !inside_ellipse(5, 3, 0, 0, 1 / 16, 1 / 32, 0)
        @test inside_ellipse(0, 0, 0, 0, 5, 6.2, 0)
        @test inside_ellipse(1, 2, 0, 0, 1 / 37, 1 / 36, -1 / 80)
    end

    @testset "inside rectangle" begin
        @test inside_rectangle(0, 0, 3, 4, 0)
        @test !inside_rectangle(10, 10, 3, 4, 0)
    end

    @testset "type stability" begin
        @inferred area_arc(0, 0, 0, 0, 10)
        @inferred area_triangle(0, 0, 1, 0, 0, 1)
        @inferred inside_triangle(0, 0, 0, 0, 0, 1, 1, 0)
        @inferred inside_ellipse(0, 0, 5, 5, 1, 1, 1)
        @inferred inside_rectangle(0, 0, 3, 4, 0)
    end

end # overlap utils


text:
#using Pkg
#Pkg.add("BenchmarkTools")
using BenchmarkTools
using LinearAlgebra
A = rand(100,100)
@btime norm(A)
function timing_from_a_function()
B = rand(10,10)
@btime norm($B)
end
timing_from_a_function()
A = rand(10_000,10_000)
B = @benchmark(norm($A))


text:
using Statistics

function quantile2(mat::AbstractArray, probs::AbstractVector; byrow::Bool = true)
    (nrow, ncol) = size(mat);
    nprobs = length(probs);

    if byrow
        res = fill(0.0, nrow, nprobs);
        for i = 1:nrow
            res[i, :] = quantile(mat[i, :], probs);
        end
    else
        res = fill(0.0, ncol, nprobs);
        for j = 1:ncol
            res[j, :] = quantile(mat[:, j], probs);
        end
    end
    res;
end

function quantile2_nan(mat::AbstractArray, probs::AbstractVector; byrow::Bool = true)
    (nrow, ncol) = size(mat);
    nprobs = length(probs);

    if byrow
        res = fill(0.0, nrow, nprobs);
        for i = 1:nrow
            res[i, :] = quantile(skipmissing(mat[i, :]), probs);
        end
    else
        res = fill(0.0, ncol, nprobs);
        for j = 1:ncol
            res[j, :] = quantile(skipmissing(mat[:, j]), probs);
        end
    end
    res;
end

max_stars_repo_name: CUG-hydro/heatwave

max_issues_repo_name: CUG-hydro/heatwave

max_forks_repo_name: CUG-hydro/heatwave


text:
function myBubbleSort(x)
if length(x) < 1
return "Your original is not an array!"
else
println(" original: $x")
s=x
index = [1:length(s)...]
rank = 0 * index
for i in 2:length(s), j in 1:length(s)-1
if s[j] > s[j+1]
s[j], s[j+1] = s[j+1], s[j]
index[j], index[j+1] = index[j+1], index[j]
end
end
for i in 1:length(index)
rank[index[i]] = i
end
end
println(" sorted:   $s\n index:    $index\n rank:     $rank")
end


text:
# Additional tests of `find_package_path`. Many of these have the
# purpose to cover error cases, making them somewhat contrived.

# Not a developed package.
@test_throws ErrorException find_package_path("AutoHashEquals")

# Find package by module and path.
using Multibreak
package_path = find_package_path(Multibreak)
@test find_package_path(package_path) == package_path

# Find package by name.
if Base.Sys.islinux()
    @test find_package_path("Multibreak") == package_path
else
    # Workaround for CI path weirdnesses on Windows and Mac.
    @test splitpath(package_path)[end-2:end] == splitpath(find_package_path("Multibreak"))[end-2:end]
end

# Not a package path.
corrupt_path = joinpath(package_path, "no_such_dir")
@test_throws ErrorException find_package_path(corrupt_path)

# Unknown package.
@test_throws ErrorException find_package_path("ZeroethTest")

# Current active environment is not a package. Return the current directory.
@test find_package_path(nothing) == pwd()

# Activate a package directory and find it with `find_package_path`.
Pkg.activate("Multibreak")
if Base.Sys.islinux()
    @test find_package_path(nothing) == package_path
else
    # Workaround for CI path weirdnesses on Windows and Mac.
    @test splitpath(package_path)[end-2:end] == splitpath(find_package_path(nothing))[end-2:end]
end
Pkg.activate()

# Early Julia 1.x versions insist on installing the General registry
# when doing `activate`. It is expected not to be there in the
# LocalRegistry tests, so remove it again if it has been installed.
if isdir(joinpath(first(DEPOT_PATH), "registries", "General"))
    Pkg.Registry.rm("General")
end


text:
function electre(decisionMat::DataFrame, weights::Array{Float64,1})::ElectreResult

    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    normalizedMat = normalize(decisionMat)
    weightednormalizedMat = w * normalizedMat

    fitnessTable = []
    nonfitnessTable = []

    for i in 1:nalternatives
        for j in 1:nalternatives
            if i != j
                betterlist = []
                worselist = []
                for h in 1:ncriteria
                    if weightednormalizedMat[i, h] >= weightednormalizedMat[j, h]
                        push!(betterlist, h)
                    else
                        push!(worselist, h)
                    end
                end
                fitnesstableelement = Dict(:i => i, :j => j, :set => betterlist)
                push!(fitnessTable, fitnesstableelement)

                nonfitnesstableelement = Dict(:i => i, :j => j, :set => worselist)
                push!(nonfitnessTable, nonfitnesstableelement)
            end
        end
    end

    fitnessmatrix = zeros(Float64, nalternatives, nalternatives)
    nonfitnessmatrix = zeros(Float64, nalternatives, nalternatives)

    for elements in fitnessTable
        i = elements[:i]
        j = elements[:j]
        elemlist = elements[:set]

        CC = sum(w[elemlist])
        fitnessmatrix[i, j] = CC
    end

    nonfitnessmatrix = zeros(Float64, nalternatives, nalternatives)
    for elements in nonfitnessTable
        i = elements[:i]
        j = elements[:j]
        elemlist = elements[:set]

        r_ik       = weightednormalizedMat[i, elemlist]
        r_jk       = weightednormalizedMat[j, elemlist]
        r_ik_full  = weightednormalizedMat[i, :]
        r_jk_full  = weightednormalizedMat[j, :]

        nom = maximum(abs.(r_ik - r_jk))
        dom = maximum(abs.(r_ik_full - r_jk_full))
        nonfitnessmatrix[i, j] = nom / dom
    end

    C = zeros(Float64, nalternatives)
    D = zeros(Float64, nalternatives)
    for i in 1:nalternatives
        C[i] = sum(fitnessmatrix[i, :]) - sum(fitnessmatrix[:, i])
        D[i] = sum(nonfitnessmatrix[i, :]) - sum(nonfitnessmatrix[:, i])
    end

    best_C_index = sortperm(C) |> last
    best_D_index = sortperm(D) |> first
    best = nothing

    if best_C_index == best_D_index
        best = (best_C_index,)
    else
        best = (best_C_index, best_D_index)
    end

    result = ElectreResult(
        decisionMat,
        w,
        weightednormalizedMat,
        fitnessTable,
        nonfitnessTable,
        fitnessmatrix,
        nonfitnessmatrix,
        C,
        D,
        best
    )

    return result
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM