text:
# cat.jl
# test

using LinearMaps: LinearMap
using LinearMapsAA: LinearMapAA, LinearMapAM, LinearMapAO, LinearMapAX
using LinearMapsAA: redim, undim
using LinearAlgebra: I
using Test: @test, @testset

# test

#=
    this approach using eval() works only in the global scope!

    N = 6
    forw = x -> [cumsum(x); 0] # non-square to stress test
    back = y -> reverse(cumsum(reverse(y[1:N])))
    prop = (name="cumsum", extra=1)
    A = LinearMapAA(forw, back, (N+1, N), prop)

    list1 = [
        :([A I]), :([I A]), :([2*A 3*A]),
        :([A; I]), :([I; A]), :([2*A; 3*A]),
        :([A A I]), :([A I A]), :([2*A 3*A 4*A]),
        :([A; A; I]), :([A; I; A]), :([2*A; 3*A; 4*A]),
        :([I A I]), :([I A A]),
        :([I; A; I]), :([I; A; A]),
    #   :([I I A]), :([I; I; A]), # unsupported
        :([A A; A A]), :([A I; 2A 3I]), :([I A; 2I 3A]),
    #   :([I A; A I]), :([A I; I A]), # weird sizes
        :([A I A; 2A 3I 4A]), :([I A I; 2I 3A 4I]),
    #   :([A I A; I A A]), :([I A 2A; 3A 4I 5A]), # weird
    #   :([I I A; I A I]), # unsupported
        :([A A I; 2A 3A 4I]),
        :([A I I; 2A 3I 4I]),
    ]

    M = Matrix(A)
    list2 = [
        :([M I]), :([I M]), :([2*M 3*M]),
        :([M; I]), :([I; M]), :([2*M; 3*M]),
        :([M M I]), :([M I M]), :([2*M 3*M 4*M]),
        :([M; M; I]), :([M; I; M]), :([2*M; 3*M; 4*M]),
        :([I M I]), :([I M M]),
        :([I; M; I]), :([I; M; M]),
    #   :([I I M]), :([I; I; M]), # unsupported
        :([M M; M M]), :([M I; 2M 3I]), :([I M; 2I 3M]),
    #   :([I M; M I]), :([M I; I M]), # weird sizes
        :([M I M; 2M 3I 4M]), :([I M I; 2I 3M 4I]),
    #   :([M I M; I M M]), :([I M 2M; 3M 4I 5M]), # weird
    #   :([I I M; I M I]), # unsupported
        :([M M I; 2M 3M 4I]),
        :([M I I; 2M 3I 4I]),
    ]

    length(list1) != length(list2) && throw("bug")

    for ii in 1:length(list1)
        e1 = list1[ii]
        b1 = eval(e1)
        e2 = list2[ii]
        b2 = eval(e2)
        @test b1 isa LinearMapAX
    end
=#

"""
    LinearMapAA_test_cat(A::LinearMapAM)
test hcat vcat hvcat
"""
function LinearMapAA_test_cat(A::LinearMapAM)
    Lm = LinearMap{eltype(A)}(x -> A*x, y -> A'*y, size(A,1), size(A,2))
    M = Matrix(A)

#    LinearMaps supports *cat of LM and UniformScaling only in v2.6.1
#    but see: https://github.com/Jutho/LinearMaps.jl/pull/71
#    B0 = [M Lm] # fails!

#=
    # cannot get cat with AbstractMatrix to work
    M1 = reshape(1:35, N+1, N-1)
    H2 = [A M1]
    @test H2 isa LinearMapAX
    @test Matrix(H2) == [Matrix(A) H2]
    H1 = [M1 A]
    @test H1 isa LinearMapAX
    @test Matrix(H1) == [M1 Matrix(A)]

    M2 = reshape(1:(3*N), 3, N)
    V1 = [M2; A]
    @test V1 isa LinearMapAX
    @test Matrix(V1) == [M2; Matrix(A)]
    V2 = [A; M2]
    @test V2 isa LinearMapAX
    @test Matrix(V2) == [Matrix(A); M2]
=#

    fun0 = A -> [
        [A I], [I A], [2A 3A],
        [A; I], [I; A], [2A; 3A],
        [A A I], [A I A], [2A 3A 4A],
        [A; A; I], [A; I; A], [2A; 3A; 4A],
        [I A I], [I A A],
        [I; A; I], [I; A; A],
    #   [I I A], [I; I; A], # unsupported
        [A A; A A], [A I; 2A 3I], [I A; 2I 3A],
    #   [I A; A I], [A I; I A], # weird sizes
        [A I A; 2A 3I 4A], [I A I; 2I 3A 4I],
    #   [A I A; I A A], [I A 2A; 3A 4I 5A], # weird
    #   [I I A; I A I], # unsupported
        [A A I; 2A 3A 4I],
        [A I I; 2A 3I 4I],
    #   [A Lm], # need one LinearMap test for codecov (see below)
    ]

    list1 = fun0(A)
    list2 = fun0(M)

    if true # need one LinearMap test for codecov
        push!(list1, [A Lm])
        push!(list2, [M M]) # trick because [M Lm] fails
    end

    for ii in 1:length(list1)
        b1 = list1[ii]
        b2 = list2[ii]
        @test b1 isa LinearMapAX
        @test b2 isa AbstractMatrix
        @test Matrix(b1) == b2
        @test Matrix(b1') == Matrix(b1)'
    end

    true
end

"""
    LinearMapAA_test_cat(A::LinearMapAO)
test hcat vcat hvcat for AO
"""
function LinearMapAA_test_cat(A::LinearMapAO)
    M = Matrix(A)

    @testset "cat AO" begin
        H = [A A]
        V = [A; A]
        B = [A 2A; 3A 4A]
        @test H isa LinearMapAO
        @test V isa LinearMapAO
        @test B isa LinearMapAM # !!
        @test Matrix(H) == [M M]
        @test Matrix(V) == [M; M]
        @test Matrix(B) == [M 2M; 3M 4M]
    end

    @testset "cat OI" begin
        @test [A I] isa LinearMapAM
        @test [A; I] isa LinearMapAM
        @test [A A; I I] isa LinearMapAM
    end

    @testset "cat AO mixed" begin
        B = redim(A, idim=(1,A._idim...)) # force incompatible dim
        @test [A B] isa LinearMapAM
        @test [A; B] isa LinearMapAM
        Z = [A undim(A) I A._lmap] # Matrix(B)]
        @test Z isa LinearMapAM
        @test Z.hcat == "OAIL"
    end

    true
end

N = 6; M = 8 # non-square to stress test
forw = x -> [cumsum(x); 0; 0]
back = y -> reverse(cumsum(reverse(y[1:N])))
A = LinearMapAA(forw, back, (M, N))
@test LinearMapAA_test_cat(A)

forw = x -> [cumsum(x ; dims=2) [0] [0]]
back = y -> reverse(cumsum(reverse(y[[1],1:N] ; dims=2) ; dims=2) ; dims=2)
O = LinearMapAA(forw, back, (M, N) ; idim=(1,N), odim=(1,M))
@test LinearMapAA_test_cat(O)


text:
module GenerateTests

using DocumenterTools
using Test
using Pkg
using Example

@testset "Generate" begin
    function check_docdir(path; format)
        @test isdir(path)
        @test isfile(joinpath(path, "mkdocs.yml")) == (format === :markdown)
        @test isfile(joinpath(path, ".gitignore"))
        @test isfile(joinpath(path, "Project.toml"))
        @test isfile(joinpath(path, "make.jl"))
        @test isdir(joinpath(path, "src"))
        @test isfile(joinpath(path, "src", "index.md"))
    end
    mktempdir() do tmp; cd(tmp) do
        # generate from module argument
        Pkg.generate("Pkg1")
        Pkg.develop(PackageSpec(path = "Pkg1"))
        @eval using Pkg1
        DocumenterTools.generate(Pkg1; format = :html)
        check_docdir(joinpath("Pkg1", "docs"); format = :html)
        @test_throws ArgumentError DocumenterTools.generate(Pkg1)
        Pkg.rm(PackageSpec("Pkg1"))

        # generate from path
        Pkg.generate("Pkg2")
        DocumenterTools.generate(joinpath("Pkg2", "docs"); format = :markdown)
        check_docdir(joinpath("Pkg2", "docs"); format = :markdown)
        @test_throws ArgumentError DocumenterTools.generate(joinpath("Pkg2", "docs"))

        # generate where name can't be determined
        mkdir("Pkg3")
        @test_throws ArgumentError DocumenterTools.generate(joinpath("Pkg3", "docs"))
        mkdir(joinpath("Pkg3", "src"))
        @test_throws ArgumentError DocumenterTools.generate(joinpath("Pkg3", "docs"))
        write(joinpath("Pkg3", "src", "Pkg3.jl"), "module Pkg3\nend\n")
        write(joinpath("Pkg3", "src", "Pkg4.jl"), "module Pkg4\nend\n")
        @test_throws ArgumentError DocumenterTools.generate(joinpath("Pkg3", "docs"))
        DocumenterTools.generate(joinpath("Pkg3", "docs"); name = "Pkg3", format = :pdf)
        check_docdir(joinpath("Pkg3", "docs"), format = :pdf)

        # throw for a package that is installed and not deved
        @test_throws ArgumentError DocumenterTools.generate(Example)

        # throw for a submodule
        @test_throws ArgumentError DocumenterTools.generate(Broadcast)
    end; end
end

end

ext: jl


text:
@testset "Check miscellaneous utility functions" begin
    @test DSGE.abbrev_symbol(:aaaaa, 3) == "aaa"
    v = [1,2,3]
    DSGE.sorted_list_insert!(v, 4)
    @test v == [1,2,3,4]
    @test DSGE.quarter_range(DSGE.quartertodate("2000-Q1"), DSGE.quartertodate("2000-Q2")) == [Date("2000-03-31",
                                                                                "yyyy-mm-dd"),
                                                                                Date("2000-06-30",
                                                                                "yyyy-mm-dd")]
end

@testset "Check detexify" begin
    @test DSGE.detexify(:α) == :alpha
    @test DSGE.detexify(:β) == :beta
    @test DSGE.detexify(:γ) == :gamma
    @test DSGE.detexify(:δ) == :delta
    @test DSGE.detexify(:ϵ) == :epsilon
    @test DSGE.detexify(:ε) == :epsilon
    @test DSGE.detexify(:ζ) == :zeta
    @test DSGE.detexify(:η) == :eta
    @test DSGE.detexify(:θ) == :theta
    @test DSGE.detexify(:ι) == :iota
    @test DSGE.detexify(:κ) == :kappa
    @test DSGE.detexify(:λ) == :lambda
    @test DSGE.detexify(:μ) == :mu
    @test DSGE.detexify(:ν) == :nu
    @test DSGE.detexify(:ξ) == :xi
    @test DSGE.detexify(:π) == :pi
    @test DSGE.detexify(:ρ) == :rho
    @test DSGE.detexify(:ϱ) == :rho
    @test DSGE.detexify(:σ) == :sigma
    @test DSGE.detexify(:ς) == :sigma
    @test DSGE.detexify(:τ) == :tau
    @test DSGE.detexify(:υ) == :upsilon
    @test DSGE.detexify(:ϕ) == :phi
    @test DSGE.detexify(:φ) == :varphi
    @test DSGE.detexify(:χ) == :chi
    @test DSGE.detexify(:ψ) == :psi
    @test DSGE.detexify(:ω) == :omega

    @test DSGE.detexify(:Α) == :Alpha
    @test DSGE.detexify(:Β) == :Beta
    @test DSGE.detexify(:Γ) == :Gamma
    @test DSGE.detexify(:Δ) == :Delta
    @test DSGE.detexify(:Ε) == :Epsilon
    @test DSGE.detexify(:Ζ) == :Zeta
    @test DSGE.detexify(:Η) == :Eta
    @test DSGE.detexify(:Θ) == :Theta
    @test DSGE.detexify(:Ι) == :Iota
    @test DSGE.detexify(:Κ) == :Kappa
    @test DSGE.detexify(:Λ) == :Lambda
    @test DSGE.detexify(:Μ) == :Mu
    @test DSGE.detexify(:Ν) == :Nu
    @test DSGE.detexify(:Ξ) == :Xi
    @test DSGE.detexify(:Π) == :Pi
    @test DSGE.detexify(:Ρ) == :Rho
    @test DSGE.detexify(:Σ) == :Sigma
    @test DSGE.detexify(:Τ) == :Tau
    @test DSGE.detexify(:Υ) == :Upsilon
    @test DSGE.detexify(:Φ) == :Phi
    @test DSGE.detexify(:Χ) == :Chi
    @test DSGE.detexify(:Ψ) == :Psi
    @test DSGE.detexify(:Ω) == :Omega

    @test DSGE.detexify("α") == "alpha"
    @test DSGE.detexify("β") == "beta"
    @test DSGE.detexify("γ") == "gamma"
    @test DSGE.detexify("δ") == "delta"
    @test DSGE.detexify("ϵ") == "epsilon"
    @test DSGE.detexify("ε") == "epsilon"
    @test DSGE.detexify("ζ") == "zeta"
    @test DSGE.detexify("η") == "eta"
    @test DSGE.detexify("θ") == "theta"
    @test DSGE.detexify("ι") == "iota"
    @test DSGE.detexify("κ") == "kappa"
    @test DSGE.detexify("λ") == "lambda"
    @test DSGE.detexify("μ") == "mu"
    @test DSGE.detexify("ν") == "nu"
    @test DSGE.detexify("ξ") == "xi"
    @test DSGE.detexify("π") == "pi"
    @test DSGE.detexify("ρ") == "rho"
    @test DSGE.detexify("ϱ") == "rho"
    @test DSGE.detexify("σ") == "sigma"
    @test DSGE.detexify("ς") == "sigma"
    @test DSGE.detexify("τ") == "tau"
    @test DSGE.detexify("υ") == "upsilon"
    @test DSGE.detexify("ϕ") == "phi"
    @test DSGE.detexify("φ") == "varphi"
    @test DSGE.detexify("χ") == "chi"
    @test DSGE.detexify("ψ") == "psi"
    @test DSGE.detexify("ω") == "omega"

    @test DSGE.detexify("Α") == "Alpha"
    @test DSGE.detexify("Β") == "Beta"
    @test DSGE.detexify("Γ") == "Gamma"
    @test DSGE.detexify("Δ") == "Delta"
    @test DSGE.detexify("Ε") == "Epsilon"
    @test DSGE.detexify("Ζ") == "Zeta"
    @test DSGE.detexify("Η") == "Eta"
    @test DSGE.detexify("Θ") == "Theta"
    @test DSGE.detexify("Ι") == "Iota"
    @test DSGE.detexify("Κ") == "Kappa"
    @test DSGE.detexify("Λ") == "Lambda"
    @test DSGE.detexify("Μ") == "Mu"
    @test DSGE.detexify("Ν") == "Nu"
    @test DSGE.detexify("Ξ") == "Xi"
    @test DSGE.detexify("Π") == "Pi"
    @test DSGE.detexify("Ρ") == "Rho"
    @test DSGE.detexify("Σ") == "Sigma"
    @test DSGE.detexify("Τ") == "Tau"
    @test DSGE.detexify("Υ") == "Upsilon"
    @test DSGE.detexify("Φ") == "Phi"
    @test DSGE.detexify("Χ") == "Chi"
    @test DSGE.detexify("Ψ") == "Psi"
    @test DSGE.detexify("Ω") == "Omega"

    @test DSGE.detexify("′") == "'"
end

@testset "Check speye" begin
    @test DSGE.speye(2) == SparseMatrixCSC{Float64}(I, 2, 2)
    @test DSGE.speye(Complex, 2) == SparseMatrixCSC{Complex}(I, 2, 2)
end

@info "The following warning is expected."
@testset "Check macro for approximate equality of matrices" begin
    @test !(DSGE.test_matrix_eq2([0. 0.; 1. 2.] + [0. 0.; 1e-9 0.], [0. 0.; 1. 2.], "a", "a", 1e-10, 1e-10))
    @test DSGE.test_matrix_eq2([0. 0.; 1. 2.] + [0. 0.; 1e-9 0.], [0. 0.; 1. 2.], "a", "a", 1e-10, 2e-7)
end

@testset "Check overloading of operators for Complex numbers" begin
    @test Complex(1.,1.) > Complex(.99999, 1e10)
    @test 1. > Complex(.99999, 1e10)
    @test Complex(.99999, 1e10) < 1.
    @test min(Complex(1., 1.), Complex(.99999, 1e10)) == 0.99999
    @test min(1., Complex(.99999, 1e10)) == 0.99999
    @test min(Complex(1., 1.), 0.99999) == 0.99999
    @test max(Complex(1., 1.), Complex(.99999, 1e10)) == 1.
    @test max(1., Complex(.99999, 1e10)) == 1.
    @test max(Complex(1., 1.), 0.99999) == 1.
end

nothing


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

using Distributed
using BenchmarkTools
using Random
using Test
rmprocs(workers())
addprocs(4)
nworkers()

@everywhere using SharedArrays

# +
# shared array
function compute_rho( steps)
    nx = 64
    np = 1000000
    xmin, xmax = -6, 6
    Lx = xmax - xmin
    rng = MersenneTwister(123)
    xp = randn(rng, np);
    wp = similar(xp)
    fill!(wp, 1.0 / np)
    rho = SharedArray(zeros(Float64, nx))

    for step = 1:steps

        fill!(rho, 0.0)
        @sync @distributed (+) for i in 1:length(rho)
            x_norm = (xp[i]-xmin) / Lx
            ip = trunc(Int,  x_norm * nx)+1
            rho[ip] += wp[i]
        end

    end

    rho

end

@time compute_rho(100);
# -

@test sum(rho) ≈ 1.0


text:
#=
Helper functions that works around the fact, that there is no generic
Table interface for this functionality. Once this is in e.g. Tables.jl,
it should be removed from GeometryBasics!
=#

"""
    attributes(hasmeta)
Returns all attributes of meta as a Dict{Symbol, Any}.
Needs to be overloaded, and returns empty dict for non overloaded types!
Gets overloaded by default for all Meta types.
"""
function attributes(hasmeta)
    return Dict{Symbol, Any}()
end

function attributes(hasmeta::StructArray)
    return Dict{Symbol, Any}((name => getproperty(hasmeta, name) for name in propertynames(hasmeta)))
end

"""
    getcolumns(t, colnames::Symbol...)

Gets a column from any Array like (Table/AbstractArray).
For AbstractVectors, a column will be the field names of the element type.
"""
function getcolumns(tablelike, colnames::Symbol...)
    return getproperty.((tablelike,), colnames)
end

getcolumn(t, colname::Symbol) = getcolumns(t, colname)[1]

"""
    MetaType(::Type{T})

Returns the Meta Type corresponding to `T`
E.g:
```julia
MetaType(Point) == PointMeta
```
"""
MetaType(::Type{T}) where T = error("No Meta Type for $T")

"""
    MetaFree(::Type{T})

Returns the original type containing no metadata for `T`
E.g:
```julia
MetaFree(PointMeta) == Point
```
"""
MetaFree(::Type{T}) where T = error("No meta free Type for $T")

"""
    meta(x::MetaObject)

Returns the metadata of `x`
"""
meta(x::T) where T = error("$T has no meta!")

metafree(x::T) where T = x

macro meta_type(name, mainfield, supertype, params...)
    MetaName = Symbol("$(name)Meta")
    field = QuoteNode(mainfield)
    NoParams = Symbol("$(MetaName)NoParams")

    params_sym = map(params) do param
        param isa Symbol && return param
        param isa Expr && param.head == :(<:) && return param.args[1]
        error("Unsupported type parameter: $(param)")
    end

    expr = quote
        struct $MetaName{$(params...), Typ <: $supertype{$(params_sym...)}, Names, Types} <: $supertype{$(params_sym...)}
            main::Typ
            meta::NamedTuple{Names, Types}
        end

        const $NoParams{Typ, Names, Types} = $MetaName{$(params_sym...), Typ, Names, Types} where {$(params_sym...)}

        function Base.getproperty(x::$MetaName{$(params_sym...), Typ, Names, Types},
                                  field::Symbol) where {$(params...), Typ, Names, Types}
            field === $field && return getfield(x, :main)
            field === :main && return getfield(x, :main)
            Base.sym_in(field, Names) && return getfield(getfield(x, :meta), field)
            error("Field $field not part of Element")
        end

        function GeometryBasics.MetaType(XX::Type{<: $supertype{$(params_sym...)} where {$(params...)}})
            return $MetaName
        end

        function GeometryBasics.MetaType(
                ST::Type{<: $supertype{$(params_sym...)}},
                ::Type{NamedTuple{Names, Types}}) where {$(params...), Names, Types}
            return $MetaName{$(params_sym...), ST, Names, Types}
        end

        GeometryBasics.MetaFree(::Type{<: $MetaName{Typ}}) where Typ = Typ
        GeometryBasics.MetaFree(::Type{<: $MetaName}) = $name
        GeometryBasics.metafree(x::$MetaName) = getfield(x, :main)
        GeometryBasics.metafree(x::AbstractVector{<: $MetaName}) = getproperty(x, $field)
        GeometryBasics.meta(x::$MetaName) = getfield(x, :meta)
        GeometryBasics.meta(x::AbstractVector{<: $MetaName}) = getproperty(x, :meta)

        function GeometryBasics.meta(main::$supertype{$(params_sym...)}; meta...) where {$(params...)}
            isempty(meta) && return elements # no meta to add!
            return $MetaName(main; meta...)
        end

        function GeometryBasics.meta(elements::AbstractVector{XX}; meta...) where XX <: $supertype{$(params_sym...)} where {$(params...)}
            isempty(meta) && return elements # no meta to add!
            n = length(elements)
            for (k, v) in meta
                if v isa AbstractVector
                    mn = length(v)
                    mn != n && error("Metadata array needs to have same length as data.
                    Found $(n) data items, and $mn metadata items")
                else
                    error("Metadata needs to be an array with the same length as data items. Found: $(typeof(v))")
                end
            end
            nt = values(meta)
            # get the first element to get the per element named tuple type
            ElementNT = typeof(map(first, nt))

            return StructArray{MetaType(XX, ElementNT)}(($(mainfield) = elements, nt...))
        end

        function GeometryBasics.attributes(hasmeta::$MetaName)
            return Dict{Symbol, Any}((name => getproperty(hasmeta, name) for name in propertynames(hasmeta)))
        end

        function (MT::Type{<: $MetaName})(args...; meta...)
            nt = values(meta)
            obj = MetaFree(MT)(args...)
            return MT(obj, nt)
        end

        function (MT::Type{<: $MetaName})(main::$(name); meta...)
            nt = values(meta)
            return MT(main, nt)
        end

        function Base.propertynames(::$MetaName{$(params_sym...), Typ, Names, Types}) where {$(params...), Typ, Names, Types}
            return ($field, Names...)
        end

        function StructArrays.staticschema(::Type{$MetaName{$(params_sym...), Typ, Names, Types}}) where {$(params...), Typ, Names, Types}
            NamedTuple{($field, Names...), Base.tuple_type_cons(Typ, Types)}
        end

        function StructArrays.createinstance(
                ::Type{$MetaName{$(params_sym...), Typ, Names, Types}},
                metafree, args...
            ) where {$(params...), Typ, Names, Types}
            $MetaName(metafree, NamedTuple{Names, Types}(args))
        end
    end
    return esc(expr)
end

@meta_type(Point, position, AbstractPoint, Dim, T)
Base.getindex(x::PointMeta, idx::Int) = getindex(metafree(x), idx)

@meta_type(NgonFace, ngon, AbstractNgonFace, N, T)
Base.getindex(x::NgonFaceMeta, idx::Int) = getindex(metafree(x), idx)

@meta_type(SimplexFace, simplex, AbstractSimplexFace, N, T)
Base.getindex(x::SimplexFaceMeta, idx::Int) = getindex(metafree(x), idx)

@meta_type(Polygon, polygon, AbstractPolygon, N, T)

@meta_type(LineString, lines, AbstractVector, P <: Line)
Base.getindex(x::LineStringMeta, idx::Int) = getindex(metafree(x), idx)
Base.size(x::LineStringMeta) = size(metafree(x))

@meta_type(MultiPoint, points, AbstractVector, P <: AbstractPoint)
Base.getindex(x::MultiPointMeta, idx::Int) = getindex(metafree(x), idx)
Base.size(x::MultiPointMeta) = size(metafree(x))

@meta_type(MultiLineString, linestrings, AbstractVector, P <: LineString)
Base.getindex(x::MultiLineStringMeta, idx::Int) = getindex(metafree(x), idx)
Base.size(x::MultiLineStringMeta) = size(metafree(x))

@meta_type(MultiPolygon, polygons, AbstractVector, P <: Polygon)
Base.getindex(x::MultiPolygonMeta, idx::Int) = getindex(metafree(x), idx)
Base.size(x::MultiPolygonMeta) = size(metafree(x))

@meta_type(Mesh, mesh, AbstractMesh, Element <: Polytope)
Base.getindex(x::MeshMeta, idx::Int) = getindex(metafree(x), idx)
Base.size(x::MeshMeta) = size(metafree(x))

"""

    MetaT(geometry, meta::NamedTuple)
    MetaT(geometry; meta...)

Returns a `MetaT` that holds a geometry and its metadata

`MetaT` acts the same as `Meta` method.
The difference lies in the fact that it is designed to handle
geometries and metadata of different/heterogeneous types.

eg: While a Point MetaGeometry is a `PointMeta`, the MetaT representation is `MetaT{Point}`
The downside being it's not subtyped to `AbstractPoint` like a `PointMeta` is.

Example:
```julia
julia> MetaT(Point(1, 2), city = "Mumbai")
MetaT{Point{2,Int64},(:city,),Tuple{String}}([1, 2], (city = "Mumbai",))
```
"""
struct MetaT{T, Names, Types}
    main::T
    meta::NamedTuple{Names, Types}
end

MetaT(x; kwargs...) = MetaT(x, values(kwargs))

"""

    metafree(x::MetaT)
    metafree(x::Array{MetaT})

Free the MetaT from metadata
i.e. returns the geometry/array of geometries
"""
function metafree(x::MetaT)
    getfield(x, :main)
end
metafree(x::AbstractVector{<: MetaT}) = map(metafree, x)

"""

    meta(x::MetaT)
    meta(x::Array{MetaT})

Returns the metadata of a `MetaT`
"""
function meta(x::MetaT)
    getfield(x, :meta)
end
meta(x::AbstractVector{<: MetaT}) = map(meta, x)

# helper methods
function Base.getproperty(x::MetaT, field::Symbol)
    if field == :main
        metafree(x)
    elseif field == :meta
        meta(x)
    else
        getproperty(meta(x), field)
    end
end

Base.propertynames(x::MetaT) = (:main, propertynames(meta(x))...)
getnamestypes(::Type{MetaT{T, Names, Types}}) where {T, Names, Types} = (T, Names, Types)

# explicitly give the "schema" of the object to StructArrays
function StructArrays.staticschema(::Type{F}) where {F<:MetaT}
    T, names, types = getnamestypes(F)
    NamedTuple{(:main, names...), Base.tuple_type_cons(T, types)}
end

# generate an instance of MetaT type
function StructArrays.createinstance(::Type{F}, x, args...) where {F<:MetaT}
     T , names, types = getnamestypes(F)
     MetaT(x, NamedTuple{names, types}(args))
end

"""
Puts an iterable of MetaT's into a StructArray
"""
function meta_table(iter)
    cols = Tables.columntable(iter)
    meta_table(first(cols), Base.tail(cols))
end

function meta_table(main, meta::NamedTuple{names, types}) where {names, types}
    F = MetaT{eltype(main), names, StructArrays.eltypes(types)}
    return StructArray{F}(; main=main, meta...)
end

Base.getindex(x::MetaT, idx::Int) = getindex(metafree(x), idx)
Base.size(x::MetaT) = size(metafree(x))


text:
#=

  "TRIBE: Trait and Range Interspecific Biogeographic Evolution" submodule package

=#

module TRIBE

using Random: randexp
using DelimitedFiles: readdlm, writedlm
using ProgressMeter: Progress, next!
using Optim: minimizer, optimize, Options
using LinearAlgebra: BLAS.axpy!, BLAS.gemm!, eigvecs, eigvals, diagm
using RCall: reval, @rput

# other submodules dependencies
using ..Utils

include("tribe/tribe_utils.jl")
include("tribe/cont_DA_prop.jl")
include("tribe/disc_DA_prop.jl")
include("tribe/tribe_data_initializer.jl")
include("tribe/area_lineage_averages.jl")
include("tribe/tribe_loglik_functions.jl")
include("tribe/mcmc_utils.jl")
include("tribe/mcmc.jl")
include("tribe/parameter_updates.jl")
include("tribe/mcmc_burn.jl")
include("tribe/under_prior.jl")
include("tribe/proposal_functions.jl")
include("tribe/tribe_wrapper.jl")
include("tribe/tribe_simulation_utils.jl")

end # module TRIBE

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
abstract type HMCIntegrator end

@with_kw mutable struct LeapfrogIntegrator <: HMCIntegrator
    step_size::Float64 = 0.0
end

@with_kw mutable struct JitteredLeapfrogIntegrator <: HMCIntegrator
    step_size::Float64 = 0.0
    jitter_rate::Float64 = 1.0
end

@with_kw mutable struct TemperedLeapfrogIntegrator <: HMCIntegrator
    step_size::Float64 = 0.0
    tempering_rate::Float64 = 1.05
end

function AHMCIntegrator(integrator::LeapfrogIntegrator)
    return AdvancedHMC.Leapfrog(integrator.step_size)
end

function AHMCIntegrator(integrator::JitteredLeapfrogIntegrator)
    return AdvancedHMC.JitteredLeapfrog(integrator.step_size, integrator.jitter_rate)
end

function AHMCIntegrator(integrator::TemperedLeapfrogIntegrator)
    return AdvancedHMC.TemperedLeapfrog(integrator.step_size, integrator.tempering_rate)
end


text:
function compute_world_X_depth()
    # These values are rather undocumented and come from the source of the
    # glpclview example in libfreenect.
    fx = 594.21
    fy = 591.04
    a = -0.0030711
    b = 3.3309495
    cx = 339.5
    cy = 242.7
    W = [1 / fx 0       0  -cx / fx
         0      -1 / fy 0  cy / fy
         0      0       0  -1
         0      0       a  b]
    P = [ 0.  0. -1.  0.
         -1.  0.  0.  0.
          0.  1.  0.  0.
          0.  0.  0.  1.]

    return P * W
end
const world_X_depth = compute_world_X_depth()

"""
$(TYPEDSIGNATURES)

Helper for getting a pointcloud from the camera at `index`.

The resulting point cloud is relative to the camera, with X forward, Y left,
and Z up.

This function uses a fixed homography matrix - you may have better results
for your own hardware by calibrating your Kinect.
"""
function sync_get_pointcloud(index::Integer)
    depth, timestamp = sync_get_depth(index, depth_11bit)
    cloud = zeros((3, 480, 640))
    for i ∈ 1:480, j ∈ 1:640
        uvw = Float64[j - 1, i - 1, depth[i, j], 1.0]
        xyzw = world_X_depth * uvw
        cloud[:, i, j] .= xyzw[1:3] ./ xyzw[4]
    end

    return cloud, timestamp
end
export sync_get_pointcloud


text:
export defaultRecoParams

function defaultRecoParams()
  params = Dict{Symbol,Any}()
  params[:reco] = "direct"
  params[:reconSize] = (32,32)
  params[:sparseTrafoName] = "Wavelet"
  params[:regularization] = "L1"
  params[:λ] = 0.0
  params[:normalizeReg] = false
  params[:solver] = "admm"
  params[:ρ] = 5.e-2
  params[:iterations] = 30

  return params
end

function setupDirectReco(acqData::AcquisitionData, recoParams::Dict)
  reconSize = recoParams[:reconSize]
  weights = samplingDensity(acqData,recoParams[:reconSize])
  # field map
  cmap = get(recoParams, :cmap, ComplexF64[])

  return reconSize, weights, cmap
end

# """
# Auxilary struct that holds parameters relevant for image reconstruction

# # Fields
# * `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
# * `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
# * `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
# * `reg::Regularization`                 - Regularization to be used
# * `normalize::Bool`                     - adjust regularization parameter according to the size of k-space data
# * `solvername::String`                  - name of the solver to use
# * `senseMaps::Array{ComplexF64}`        - coil sensitivities
# * `correctionMap::Array{ComplexF64}`    - fieldmap for the correction of off-resonance effects
# * `method::String="nfft"`               - method to use for time-segmentation when correctio field inhomogeneities
# """
# mutable struct RecoParameters{N}
#   reconSize::NTuple{N,Int64}
#   weights::Vector{Vector{ComplexF64}}
#   sparseTrafo::AbstractLinearOperator
#   reg::Vector{Regularization}
#   normalize::Bool
#   encodingOps::
#   solvername::String
#   senseMaps::Array{ComplexF64}
# end

"""
    setupIterativeReco(acqData::AcquisitionData, recoParams::Dict)

builds relevant parameters and operators from the entries in `recoParams`

# relevant parameters
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `reg::Regularization`                 - Regularization to be used
* `normalize::Bool`                     - adjust regularization parameter according to the size of k-space data
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
* `correctionMap::Array{ComplexF64}`    - fieldmap for the correction of off-resonance effects
* `method::String="nfft"`               - method to use for time-segmentation when correctio field inhomogeneities

`sparseTrafo` and `reg` can also be speficied using their names in form of a string.
"""
function setupIterativeReco(acqData::AcquisitionData, recoParams::Dict)

  red3d = dims(trajectory(acqData,1))==2 && length(recoParams[:reconSize])==3
  if red3d  # acqData is 3d data converted to 2d
    reconSize = (recoParams[:reconSize][2], recoParams[:reconSize][3])
    recoParams[:reconSize] = reconSize
  else  # regular case
    reconSize = recoParams[:reconSize]
  end

  # density weights
  densityWeighting = get(recoParams,:densityWeighting,true)
  if densityWeighting
    weights = samplingDensity(acqData,reconSize)
  else
    numContr = numContrasts(acqData)
    weights = Array{Vector{ComplexF64}}(undef,numContr)
    for contr=1:numContr
      numNodes = size(acqData.kdata[contr],1)
      weights[contr] = [1.0/sqrt(prod(reconSize)) for node=1:numNodes]
    end
  end

  # sparsifying transform
  if haskey(recoParams,:sparseTrafo) && typeof(recoParams[:sparseTrafo]) != String
    sparseTrafo = recoParams[:sparseTrafo]
  elseif haskey(recoParams,:sparseTrafo)
    sparseTrafoName = get(recoParams, :sparseTrafo, "nothing")
    sparseTrafo = SparseOp(sparseTrafoName, reconSize; recoParams...)
  else
    sparseTrafo=nothing
  end

  # bare regularization (without sparsifying transform)
  regName = get(recoParams, :regularization, "L1")
  λ = get(recoParams,:λ,0.0)
  reg = Regularization(regName, λ; shape=reconSize, recoParams...)

  # normalize regularizer ?
  normalize = get(recoParams, :normalizeReg, false)

  encOps = get(recoParams, :encodingOps, nothing)

  # solvername
  solvername = get(recoParams, :solver, "fista")

  # sensitivity maps
  senseMaps = get(recoParams, :senseMaps, ComplexF64[])
  if red3d && !isempty(senseMaps) # make sure the dimensions match the trajectory dimensions
    senseMaps = permutedims(senseMaps,[2,3,1,4])
  end

  return reconSize, weights, sparseTrafo, vec(reg), normalize, encOps, solvername, senseMaps
end

function getEncodingOperatorParams(; kargs...)
  encKeys = [:correctionMap, :method, :toeplitz, :oversamplingFactor, :kernelSize, :K, :K_tol]
  return Dict([key=>kargs[key] for key in intersect(keys(kargs),encKeys)])
end


text:
@everywhere using SeisDownload
using Dates, LightXML

#==================================================#
# Input Parameters
MAX_MEM_PER_CPU = 2.0 # [GB] maximum allocated memory for one cpu
DownloadType = "Earthquake" # Choise of "Noise" or "Earthquake"

network     = ["*"] #Currently only one network is available when you want to specify.
station     = ["*"]
location    = ["*"]
channel     = ["BH?"]

#Specify region [lat_min, lat_max, lon_min, lon_max, dep_min, dep_max], with lat, lon in decimal degrees (°) and depth in km with + = down.
#dep_min and dep_max is optional.
IsLocationBox = true
locationbox   = [30.9332, 41.4757, 129.9166, 145.3711]

method      = "FDSN" # Method to download data.
datasource  = "IRIS"

SecondsBeforePick = 10 # [s] start downloading `SecondsBeforePick`[s] before picked time
SecondsAfterPick  = 10 * 60 # [s] end downloading `SecondsBeforePick`[s] after picked time

# Time info for Noise case
catalog  = "./fdsnws-event_JP.xml" #QuakeML format

IsResponseRemove = true #whether instrumental response is removed or not
pre_filt    = (0.001, 0.002, 10.0, 20.0) #prefilter tuple used obspy remove_response: taper between f1 and f2, f3 and f4 with obspy

fodir       = "./dataset"
foname      = "JPquake" # data is saved at ./dataset/$foname.jld2
#==================================================#

# Read QuakeML
# This can be customized by users

xdoc = parse_file(catalog)
xroot = root(xdoc)  # an instance of XMLElement
# print its name
ces = collect(child_elements(xroot))
ev = ces[1]["event"]
Ne = length(ev)

event = [] # event information (starttime, lat, lon, mag...)

for i = 1:Ne

    Npicks = 1
    #println("index of quake ", i, " Number of picks  ",Npicks)

    #event information
    publicid    = attribute(ev[i]["origin"][1], "publicID")
    catalog     = attribute(ev[i]["origin"][1], "catalog")
    time        = content(ev[i]["origin"][1]["time"][1])[2:end-1]
    latitude    = content(ev[i]["origin"][1]["latitude"][1])[2:end-1]
    longitude   = content(ev[i]["origin"][1]["longitude"][1])[2:end-1]
    depth       = content(ev[i]["origin"][1]["depth"][1])[2:end-1]
    magnitude   = content(ev[i]["magnitude"][1]["mag"][1])[2:end-1]
    type        = content(ev[i]["magnitude"][1]["type"][1])
    origin      = Dict("publicid" => publicid, "catalog" => catalog, "time" => time,
    "latitude"=>latitude, "longitude"=>longitude, "depth"=>depth, "magnitude"=>magnitude, "type"=>type)

    # loop through all phase picked
    # create a dictionary for each quake with channel information and time.
    ptime=Vector{String}(undef,Npicks)
    net=Vector{String}(undef,Npicks)
    sta=Vector{String}(undef,Npicks)
    cha=Vector{String}(undef,Npicks)
    loc=Vector{String}(undef,Npicks)
    src=Vector{String}(undef,Npicks)

    pickphase_dict = []

    for j=1:Npicks
        ptime[j] = content(ev[i]["origin"][1]["time"][1])[2:end-1]

        net[j] = network[1] # network code
        sta[j] = station[1] # station code
        loc[j] = location[1] # station code
        cha[j] = channel[1] # channel code
        src[j] = datasource # data center

        starttime = DateTime(ptime[j])
        starttime = starttime - Dates.Second(SecondsBeforePick)
        endtime  = starttime + Dates.Second(SecondsAfterPick)
        # prints for sanity check

        pdict_temp = Dict("starttime" => starttime, "endtime" => endtime,
         "net" => net[j], "sta" => sta[j], "loc" => loc[j], "cha" => cha[j], "src" => src[j])

        push!(pickphase_dict, pdict_temp)

    end

    event_temp = Dict("origin" => origin, "pickphase" => pickphase_dict)
    push!(event, event_temp)
end

free(xdoc)

mkpath(fodir)
fopath=joinpath(fodir, foname*".jld2")

#if lat-log box or not
IsLocationBox ? reg=locationbox : reg=[]
IsResponseRemove ? pre_filt = pre_filt : pre_filt = []

InputDictionary = Dict([
      "DownloadType"    => DownloadType,
      "method"          => method,
      "event"           => event,
      "IsLocationBox"   => IsLocationBox,
      "reg"             => reg,
      "IsResponseRemove"=> IsResponseRemove,
      "pre_filt"        => pre_filt,
      "fopath"          => fopath

    ])

# mass request with input Dictionary
seisdownload(InputDictionary, MAX_MEM_PER_CPU=float(MAX_MEM_PER_CPU))


text:
const _openers = b"([{"
const _closers = b")]}"
"""
    chkparens(cmd) -> 0 | -1 | pos

Return the position of a misbalanced parenthesis or 0 if the given string
is balanced, and -1 if there are more open parentheses than closed.
"""
function chkparens(cmd::Array{UInt8})
    # Allow a language prefix such as k) or J).
    isempty(cmd) && return 0
    if length(cmd) > 1 && cmd[2] == 0x29  # ')'
        cmd = cmd[3:end]
    end
    stack = UInt8[]
    for (pos, x) in enumerate(cmd)
        i = findnext(_openers, x, 1)
        if i > 0
            push!(stack, x)
            continue
        end
        i = findnext(_closers, x, 1)
        if i > 0 && (isempty(stack) || pop!(stack) != _openers[i])
            return pos
        end
    end
    isempty(stack) ? 0 : -1
end
chkparens(cmd::AbstractString) = chkparens(transcode(UInt8, cmd))


text:
module Calculus2
    export
        @gradexpr,
        @hessexpr,
        abs2dual,
        absdual,
        conjdual,
        deparse,
        Dual,
        dual,
        Dual128,
        dual128,
        Dual64,
        dual64,
        dual_show,
        DualPair,
        epsilon,
        grad!,
        grad,
        gradexpr,
        gradfunc!,
        gradfunc,
        hess!,
        hess,
        hessexpr,
        hessfunc!,
        hessfunc,
        integrate,
        isdual,
        isgrad!,
        isgrad,
        ishess!,
        ishess,
        ∇!,
        ∇,
        ∫

    include("utils/maxdiff.jl")

    include("symbolic/symbolic.jl")
    include("symbolic/gradexpr.jl")
    include("symbolic/hessexpr.jl")
    include("symbolic/deparse.jl")

    include("finite/epsilon.jl")
    include("finite/grad.jl")
    include("finite/jacobian.jl")
    include("finite/hess.jl")

    include("automatic/dualnumbers.jl")
    include("automatic/grad.jl")

    include("integrate/simpsons.jl")
    include("integrate/montecarlo.jl")

    include("api/grad.jl")
    include("api/gradfunc.jl")
    include("api/gradexpr.jl")
    include("api/hess.jl")
    include("api/hessfunc.jl")
    include("api/hessexpr.jl")
    include("api/checker.jl")
    include("api/integrate.jl")
    include("api/unicode.jl")
end


text:
using ControlSystems
using PredictiveControl
using Test

# This is a controllable and unstable discrete-time system
Aunstab = [2.0 1.0;
           0.0 1.0]

# This is a controllable and stable discrete-time system
Astab = [0.9 1.0
         0.0 0.9]

B = [0.0;
     1.0]
C = [1.0 0.0;
     0.0 1.0]
D = [0.0;
     0.0]
Ts = 0.1

sysStab   = StateSpace( Astab, B, C, D, Ts )
sysUnstab = StateSpace( Aunstab, B, C, D, Ts )

Q = [1.0 1.0;
     1.0 1.0]
R = [1.0]
N = 2

em = Matrix{Float64}(undef, 0, 0)

# Ensure the dimension checks on the weight matrices work
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, em,  R,  Q )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N,  Q, em,  Q )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N,  Q,  R, em )

# Define an incorrect controller
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, K = [ 1.0 ] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, K = [ 1.0; 2.0 ] )

# Define bad constraints
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, F = [ 1.0 ],     E = [ 1.0 ],      g = [ 1.0 ] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, F = [ 1.0 2.0 ], E = [ 1.0 1.0 ],  g = [ 1.0 ] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, F = [ 1.0 2.0 ], E = [ 1.0 ],      g = [ 1.0; 1.0 ] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, F = [ 1.0 2.0 ], E = [ 1.0; 1.0 ], g = [ 1.0 ] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, F = [ 1.0 2.0 ], E = [ 1.0; 1.0 ], g = [ 1.0 ] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, F = [ 1.0 2.0; 1.0 2.0], E = [ 1.0 ], g = [ 1.0 ] )

# Define bad bounds
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, xₗ = [1.0; 2.0; 3.0] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, xᵤ = [1.0; 2.0; 3.0] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, uₗ = [1.0; 2.0] )
@test_throws DimensionMismatch ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, uᵤ = [1.0; 2.0] )

# Test bounds constraints
let clqr = ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, xᵤ=1.0 )
    @test length( clqr.xᵤ ) == 2
    @test clqr.xᵤ[1] == 1.0
    @test clqr.xᵤ[2] == 1.0
end

let clqr = ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, xₗ=1.0 )
    @test length( clqr.xₗ ) == 2
    @test clqr.xₗ[1] == 1.0
    @test clqr.xₗ[2] == 1.0
end

let clqr = ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, uᵤ=1.0 )
    @test length( clqr.uᵤ ) == 1
    @test clqr.uᵤ[1] == 1.0
end

let clqr = ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, uₗ=1.0 )
    @test length( clqr.uₗ ) == 1
    @test clqr.uₗ[1] == 1.0
end

let clqr = ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, xᵤ=[Inf; 1.0] )
    @test length( clqr.xᵤ ) == 2
    @test ismissing( clqr.xᵤ[1] )
    @test clqr.xᵤ[2] == 1.0
end

let clqr = ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :Q, xₗ=[-Inf; 1.0] )
    @test length( clqr.xₗ ) == 2
    @test ismissing( clqr.xₗ[1] )
    @test clqr.xₗ[2] == 1.0
end

# Use P=Q and no prestabilization controller
let clqr = ConstrainedTimeInvariantLQR( sysUnstab, N, Q, R, :Q )
    @test iszero( clqr.K )

    @test clqr.Qₖ == clqr.Q
    @test clqr.Q  == clqr.P

    @test !isstable( PredictiveControl.getsystem( clqr, prestabilized = false ) )
    @test !isstable( PredictiveControl.getsystem( clqr ) )
end

# Use P=Q and an LQR prestabilization controller
let clqr = ConstrainedTimeInvariantLQR( sysUnstab, N, Q, R, :Q, K = :dlqr )
    @test clqr.Qₖ != clqr.Q
    @test clqr.Qₖ == clqr.P

    @test !isstable( PredictiveControl.getsystem( clqr, prestabilized = false ) )
    @test isstable( PredictiveControl.getsystem( clqr ) )
end

# Use P=(the solution to the DARE) and the LQR controller for prestabilization
let clqr = ConstrainedTimeInvariantLQR( sysUnstab, N, Q, R, :dare, K = :dlqr )
    Aₖ = Aunstab - B*clqr.K

    @test clqr.Qₖ != clqr.Q
    @test clqr.Qₖ != clqr.P
    @test clqr.P  ≈ Aₖ'*clqr.P*Aₖ + clqr.Qₖ

    @test !isstable( PredictiveControl.getsystem( clqr, prestabilized = false ) )
    @test isstable( PredictiveControl.getsystem( clqr ) )
end

# Use P=(the solution to the lyap) and no stabilization
let clqr = ConstrainedTimeInvariantLQR( sysStab, N, Q, R, :dlyap )
    @test iszero( clqr.K )

    @test clqr.Qₖ == clqr.Q
    @test clqr.Q  != clqr.P
    @test clqr.P  ≈ Astab'*clqr.P*Astab + clqr.Q
end


text:
using CoveragePackage
using Test

@testset "CoveragePackage.jl" begin
    e = CoveragePackage.ExtractKeyAsField(1, "2")
    e2 = CoveragePackage.ExtractKeyAsField(1, "2")
	@test hash(e) === hash(e2)
	@test e == e2
end

max_stars_repo_name: racinmat/CoveragePackage

max_issues_repo_name: racinmat/CoveragePackage

max_forks_repo_name: racinmat/CoveragePackage


text:
const TAB_RIGHT = Dict( "text-align"=>"right")
const TAB_LEFT = Dict( "text-align"=>"left")
const TAB_CENTRE = Dict( "text-align"=>"center")
const TAB_RED = Dict( "color"=>"#BB3311")
const TAB_GREEN = Dict( "color"=>"#33BB11")
const TAB_BOLD = Dict( "font-weight"=>"Bold")
const TAB_TOTAL = Dict( "font-weight"=>"Bold", "background-color"=>"#ccccee" )
const RESERVED = Dict( "background-color"=>"#dddddd")

function u( d ... )
    Dict( union(d...))
end

function draw_lorenz( pre::Vector, post::Vector )
    np = size(pre)[1]
    step = 1/np
    xr = 0:step:1
    # insert a zero element
    v1 = insert!(copy(pre),1,0)
    v2 = insert!(copy(post),1,0)

    layout = Layout(
        title="Lorenz Curve",
        xaxis_title="Population Share",
        yaxis_title="Net Income Share",
        xaxis_range=[0, 1],
        yaxis_range=[0, 1],
        legend=attr(x=0.01, y=0.95),
        width=350,
        height=350)

    pre = scatter(
        x=xr,
        y=v1,
        mode="line",
        name="Pre" )

    post = scatter(
        x=xr,
        y=v2,
        mode="line",
        name="Post" )

    diag = scatter(y=[0,1], x=[0,1], showlegend=false, name="")
    # return [pre, post, diag]
    #return Dict([:data=>[pre, post, diag],:layout=>layout])
    return PlotlyJS.Plot( [pre, post, diag], layout)
end

function drawDeciles( pre::Vector, post :: Vector )
    v = pre-post;
    println( "v=$v")
    layout = Layout(
        title="Gain/Loss by decile",
        xaxis_title="Decile",
        yaxis_title="£pw",
        width=350,
        height=350)
    return PlotlyJS.Plot( bar( x=1:11, y=v), layout )
end

"""
There's a pre-computed histogram in summary, but
in turns out to be less hassle to just feed plotly
the raw data & let it do its own histogram.
"""
function drawMRS( pre :: DataFrame, post :: DataFrame )
    layout = Layout(
        title="Marginal Tax Rates",
        xaxis_title="METR(%)",
        yaxis_title="People",
        width=350,
        height=350)
    preb = histogram( pre, x=:metr, y=0:5:200, histnorm="probability density")
    postb = histogram( post, x=:metr, y=0:5:200, histnorm="probability density")
    return PlotlyJS.Plot( [preb,postb], layout )
end

"""
This doesn't work well.
"""
function drawMRS( pre :: Histogram, post :: Histogram )
    layout = Layout(
        title="Marginal Tax Rates",
        xaxis_title="METR(%)",
        yaxis_title="Number",
        width=350,
        height=350)
    preb = bar( x=pre.edges, y=pre.weights)
    postb = bar( x=post.edges, y=post.weights)
    return PlotlyJS.Plot( [preb, postb], layout )
end

function thing_table( names::Vector{String}, v1::Vector, v2::Vector, up_is_good::Vector{Int} )
    table_header =
        html_thead(
            html_tr([html_th(""),
            html_th("Before",style=TAB_RIGHT),
            html_th("After",style=TAB_RIGHT),
            html_th("Change",style=TAB_RIGHT)])
        )
    rows = []
    n = size(names)[1]
    diff = v2 - v1
    for i in 1:n
        colour = Dict()
        if (up_is_good[i] !== 0) && (! (diff[i] ≈ 0))
            if diff[i] > 0
                colour = up_is_good[i] == 1 ? TAB_GREEN : TAB_RED
             else
                colour = up_is_good[i] == 1 ? TAB_RED : TAB_GREEN
            end # neg diff
        end # non zero diff
        ds = diff[i] ≈ 0 ? "-" : fp(diff[i])
        row = html_tr([html_td(names[i]),
            html_td(f2(v1[i]),style=TAB_RIGHT),
            html_td(f2(v2[i]),style=TAB_RIGHT),
            html_td( ds ,style=u(TAB_RIGHT,colour))
            ])
        push!( rows, row )
    end
    table_body = html_tbody(rows)
    return dbc_table([table_header,table_body], bordered = false)
end
const MR_UP_GOOD = [1,0,0,0,0,0,0,-1,-1]

const COST_UP_GOOD = [1,1,1,1,-1,-1,-1,-1,-1,-1,-1]

function costs_table( incs1 :: DataFrame, incs2 :: DataFrame )
    v1 = extract_incs( incs1, COST_ITEMS ) ./ 1_000_000
    v2 = extract_incs( incs2, COST_ITEMS ) ./ 1_000_000
    thing_table( COST_LABELS, v1, v2, COST_UP_GOOD )
end

function frame_to_dash_table(
    df :: DataFrame;
    up_is_good :: Vector{Int},
    prec :: Int = 2,
    caption :: String = "",
    totals_col :: Int = -1 )
    table_header =
        html_thead(
            html_tr([html_th(""),
                html_th("Before",style=TAB_RIGHT),
                html_th("After",style=TAB_RIGHT),
                html_th("Change",style=TAB_RIGHT)
            ]) # tr
        )
    rows = []
    i = 0
    for r in eachrow( df )
        i += 1
        colour = Dict()
        if (up_is_good[i] !== 0) && (! (r.Change ≈ 0))
            if r.Change > 0
                colour = up_is_good[i] == 1 ? TAB_GREEN : TAB_RED
             else
                colour = up_is_good[i] == 1 ? TAB_RED : TAB_GREEN
            end # neg diff
        end # non zero diff
        ds = r.Change ≈ 0 ? "-" : format(r.Change, commas=true, precision=prec )
        if ds != "-" && r.Change > 0
            ds = "+$(ds)"
        end
        row_style = i == totals_col ? TAB_TOTAL : Dict()
        row = html_tr( [
            html_th( r.Item, style=TAB_LEFT),
            html_td(format(r.Before, commas=true, precision=prec),style=TAB_RIGHT),
            html_td(format(r.After, commas=true, precision=prec),style=TAB_RIGHT),
            html_td( ds, style=u(TAB_RIGHT,colour))], style=row_style)
        push!( rows, row )
    end
    table_body = html_tbody(rows)
    table_caption = caption != "" ? html_caption( caption ) : nothing
    return dbc_table([table_header,table_caption,table_body], bordered = false)
end

function costs_table( incs1 :: DataFrame, incs2 :: DataFrame )
    df = costs_dataframe( incs1, incs2 )
    return frame_to_dash_table( df, prec=0, up_is_good=COST_UP_GOOD,
        caption="Tax Liabilities and Benefit Entitlements, £m pa, 2021/22" )
    # thing_table( COST_LABELS, v1, v2, COST_UP_GOOD )
end

function overall_cost( incs1:: DataFrame, incs2:: DataFrame )
    n1 = incs1[1,:net_cost]
    n2 = incs2[1,:net_cost]

    eni1 = incs1[1,:employers_ni]
    eni2 = incs2[1,:employers_ni]
    d = (n1-eni1) - (n2-eni2)
    if d ≈ 0
        return
    end
    d /= 1_000_000
    colour = "info"
    extra = ""
    change_str = "In total, your changes cost less than £1m"
    change_val = ""
    if abs(d) > 1
        change_val = f0(abs(d))
        if d > 0
            colour = "success"
            change_str = "In total, your changes raise £"
            extra = "m."
        else
            colour = "danger"
            change_str = "In total, your changes cost £"
            extra = "m."
        end
    end
    return dbc_alert( [change_str, html_strong( change_val ), extra ], color = colour )
end

function mr_table( mr1, mr2 )
    df = mr_dataframe( mr1.hist, mr2.hist, mr1.mean, mr2.mean )
    n = size(df)[1]
    table = frame_to_dash_table(
        df,
        prec=0,
        up_is_good=MR_UP_GOOD,
        caption="Working age individuals with Marginal Effective Tax Rates
                (METRs) in the given range. METR is the percentage of the next £1 you earn that is taken away in taxes or
                reduced means-tested benefits.",
        totals_col = n )
    return table

    # thing_table( MR_LABELS, mr1.weights, mr2.weights, MR_UP_GOOD)
end

function ineq_table( ineq1 :: InequalityMeasures, ineq2 :: InequalityMeasures )
    df = ineq_dataframe( ineq1, ineq2 )
    up_is_good = fill( -1, 6 )
    return frame_to_dash_table(
        df,
        prec=2,
        up_is_good=up_is_good,
        caption="Standard Inequality Measures, using Before Housing Costs Equivalised Net Income." )
end

function pov_table(
    pov1 :: PovertyMeasures,
    pov2 :: PovertyMeasures,
    ch1  :: GroupPoverty,
    ch2  :: GroupPoverty )
    df = pov_dataframe( pov1, pov2, ch1, ch2 )
    up_is_good = fill( -1, 7 )
    return frame_to_dash_table(
        df,
        prec=2,
        up_is_good=up_is_good,
        caption="Standard Poverty Measures, using Before Housing Costs Equivalised Net Income." )
end

function gain_lose_table( gl :: NamedTuple )
    losepct = md_format(100*gl.losers/gl.popn)
    gainpct = md_format(100*gl.gainers/gl.popn)
    ncpct = md_format(100*gl.nc/gl.popn)
    table_header =
        html_thead(
            html_tr([html_th(""), html_th(""),html_th("%",style=TAB_RIGHT)])
        )
    row1 = html_tr([html_th("Gainers"), html_td(f0(gl.gainers),style=TAB_RIGHT),html_td(gainpct,style=TAB_RIGHT) ])
    row2 = html_tr([html_th("Losers"), html_td(f0(gl.losers),style=TAB_RIGHT),html_td(losepct,style=TAB_RIGHT)])
    row3 = html_tr([html_th("Unchanged"), html_td(f0(gl.nc),style=TAB_RIGHT),html_td(ncpct,style=TAB_RIGHT)])
    table_body = html_tbody([row1, row2, row3])
    table_caption = html_caption( "Individuals living in households where net income has risen, fallen, or stayed the same respectively.")
    table = dbc_table([table_header,table_caption,table_body], bordered = false)
    return table
end

function rb_table( sys :: TaxBenefitSystem )
    table_header_1 =
        html_thead(
            html_tr([
                html_th("Rates(%)"),
                html_th("Threshold(£pa)")
	        ])
        )
    nr = size( sys.it.non_savings_rates)[1]
    rows = []

    rates = copy( sys.it.non_savings_rates ) .* 100
    bands = copy( sys.it.non_savings_thresholds ) .* WEEKS_PER_YEAR
    for i in 1:nr
        bs = (i < nr) ? f0(bands[i]) : "Remainder"
        row = html_tr([
            html_td(f2(rates[i]),style=TAB_RIGHT),
            html_td(bs,style=TAB_RIGHT)])
        push!( rows, row )
    end
    table_body = html_tbody(rows)
    table = dbc_table([table_header_1,table_body], bordered = false)
    println( typeof(table))
    return table
end

# as an actual table
function make_output_table_t( results::NamedTuple, sys::TaxBenefitSystem )
    header = []

    hrow_1 = html_tr([
        html_td(html_h4("Gainers and Losers"),style=TAB_CENTRE,colSpan=2),
		html_td(html_h4("Costs and Revenues (£m pa)"),style=TAB_CENTRE),
        html_td(html_h4("Incentives"),style=TAB_CENTRE)
	])

    row1 = html_tr([
        html_td(
			gain_lose_table( results.gain_lose)),
		html_td(dcc_graph(figure=drawDeciles(
			results.summary.deciles[1][:,3],
			BASE_RESULTS.summary.deciles[1][:,3]))),
        html_td(
             costs_table(
                BASE_RESULTS.summary.income_summary[1],
                results.summary.income_summary[1])),
        html_td(
            [mr_table(
                    BASE_RESULTS.summary.metrs[1],
                    results.summary.metrs[1] )
            ],style=TAB_CENTRE )
        ])

    hrow_2 = html_tr([
        html_td(html_h4("Poverty"),style=TAB_CENTRE),
		html_td(html_h4("Inequality"),style=TAB_CENTRE,colSpan=2)
	])

	row2 = html_tr([
        html_td(
			pov_table(
				BASE_RESULTS.summary.poverty[1],
				results.summary.poverty[1],
                BASE_RESULTS.results.child_poverty[1],
                results.summary.child_poverty[1])),
        html_td( ineq_table(
            BASE_RESULTS.summary.inequality[1],
            results.summary.inequality[1])),
        html_td( dcc_graph(figure=draw_lorenz(
                BASE_RESULTS.summary.deciles[1][:,2],
                results.summary.deciles[1][:,2]))),
        html_td() # spacer
    ]) # TR

    table_body = html_tbody([hrow_1, row1, hrow_2, row2 ]) #, row2
    table = dbc_table([ table_body], bordered = false)

    return table
end

# as bootstrap rows
function make_output_table(
    results :: AllOutput,
    sys     :: TaxBenefitSystem )
    row1 = dbc_row([
        dbc_col([
            html_h4("Gainers and Losers"),
            dbc_row([
                dbc_col(gain_lose_table( results.gain_lose)),
                dbc_col(
                    dcc_graph(
                        figure=drawDeciles(
                            results.summary.deciles[1][:,3],
                            BASE_RESULTS.summary.deciles[1][:,3])
                        )
                    )
            ]) # inner row for chart & table
        ]), # row1, col1
        dbc_col([
            html_h4(" Revenues and Costs"),
            costs_table(
                BASE_RESULTS.summary.income_summary[1],
                results.summary.income_summary[1]),
            overall_cost(
                BASE_RESULTS.summary.income_summary[1],
                results.summary.income_summary[1])
        ]), # row1 col2
        dbc_col([
                html_h4("Incentives - Marginal Effective Tax Rates.")
                mr_table(
                    BASE_RESULTS.summary.metrs[1],
                    results.summary.metrs[1] )
        ]) # row1 col3
    ]) # row 1

	row2 = dbc_row([
        dbc_col([
            html_h4("Poverty"),
			pov_table(
				BASE_RESULTS.summary.poverty[1],
				results.summary.poverty[1],
				BASE_RESULTS.summary.child_poverty[1],
				results.summary.child_poverty[1]
                )
        ]), # row2 col 1
        dbc_col( [
            html_h4("Inequality"),
            ineq_table(
                BASE_RESULTS.summary.inequality[1],
                results.summary.inequality[1])
        ]), # row2 col 2
        dbc_col( [
                dcc_graph(
                    figure=draw_lorenz(
                        BASE_RESULTS.summary.deciles[1][:,2],
                        results.summary.deciles[1][:,2])
                )
        ]) # row 2 col 3
    ]) # row 2
    return html_div([ row1, row2 ])
end

function it_fieldset()
    it = html_fieldset([
        html_legend( "Income Tax"),
        dbc_row([
            dbc_col(
                dbc_label("Basic Rate (%)"; html_for="basic_rate")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "basic_rate",
                    min = 0,
                    max = 100,
                    size = "4",
                    value = 20.0,
                    step = 0.5 )
            ) # col
        ]),
        dbc_row([
            dbc_col(
                dbc_label("Higher Rate (%)"; html_for="higher_rate")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "higher_rate",
                    min = 0,
                    max = 100,
                    size = "4",
                    value = 41.0,
                    step = 0.5 )
            ) # col
        ]), # row
        dbc_row([
            dbc_col(
                dbc_label("Top Rate (%)"; html_for="top_rate")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "top_rate",
                    min = 0,
                    max = 100,
                    size = "4",
                    value = 46.0,
                    step = 0.5 )
            ) # col
        ]), # row
        dbc_row([
            dbc_col(
                dbc_label("Personal Allowance £pa"; html_for="pa")
            ), # col
            dbc_col(
                dbc_input(
                    type="number",
                    id = "pa",
                    min = 0,
                    max = 50_000.0,
                    size = "4",
                    value = 12_570,
                    step = 1,
                    style=RESERVED )
            ) # col
        ]), # row
    ])
    return it
end

function ni_fieldset()
    ni = html_fieldset([
        html_legend( "National Insurance"),
        dbc_row([
            dbc_col(
                dbc_label("Employee's Rate(%)"; html_for="ni_prim")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "ni_prim",
                    min = 0,
                    max = 50,
                    size = "4",
                    value = 12.0,
                    step = 0.1,
                    style=RESERVED )
            ) # col
        ]),

        dbc_row([
            dbc_col(
                dbc_label("Employer's Rate(%)"; html_for="ni_sec")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "ni_sec",
                    min = 0,
                    max = 50,
                    size = "4",
                    value = 13.8,
                    step = 0.1,
                    style=RESERVED )
            ) # col
        ]), # row

    ])
    return ni
end

function ben_fieldset()
    bens = html_fieldset([
        html_legend( "Benefits"),
        dbc_row([
            dbc_col(
                dbc_label("Child Benefit £pw (1st child)"; html_for="cb")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "cb",
                    min = 0,
                    max = 100,
                    size = "4",
                    value = 21.15,
                    step = 0.05,
                    style=RESERVED )
            ) # col
        ]), # row
        dbc_row([
            dbc_col(
                dbc_label("New State Pension £pw"; html_for="pen")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "pen",
                    min = 0,
                    max = 500,
                    size = "4",
                    value = 179.60,
                    step = 0.10,
                    style=RESERVED )
            ) # col
        ]),
        dbc_row([
            dbc_col(
                dbc_label("Universal Credit Taper (%)"; html_for="uctaper")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "uctaper",
                    min = 0,
                    max = 100,
                    size = "4",
                    value = 55.0,
                    step = 0.5,
                    style=RESERVED )
            ) # col
        ]),
        dbc_row([
            dbc_col(
                dbc_label("Universal Credit: single 25+ adult £pm"; html_for="ucs")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "ucs",
                    min = 0,
                    max = 1000,
                    size = "4",
                    value = 324.84,
                    step = 0.01,
                    style=RESERVED )
            ) # col
        ]), # row
        dbc_row([
            dbc_col(
                dbc_label("Working Tax Credit: basic amount £pa"; html_for="wtcb")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "wtcb",
                    min = 0,
                    max = 10_000,
                    size = "6",
                    value = 2_005,
                    step = 0.50,
                    style=RESERVED )
            ) # col
        ]), # row
        dbc_row([
            dbc_col(
                dbc_label("Scottish Child Payment; £pw"; html_for="scp")
            ),
            dbc_col(
                dbc_input(
                    type="number",
                    id = "scp",
                    min = 0,
                    max = 100,
                    size = "4",
                    value = 10.0,
                    step = 0.05 )
            ) # col
        ]),
        dbc_row([
            dbc_col(
                dbc_label("Scottish Child Payment; maximum age"; html_for="scp")
            ),
            dbc_col(
                dcc_slider(
					id = "scp_age",
					min = 1,
					max = 15,
					marks = Dict([Symbol("$v") => Symbol("$v") for v in 0:15]),
					value = 5,
					step = 1
				)
            ) # col
        ]),
         # row
    ])  # fieldset
    return bens
end

function submit_button()
    dbc_button(
        id = "submit-button",
        class_name="primary",
        color = "primary",
        name = "Run",
        value = "Run",
        children = "submit"
    )
end


text:
module BiharmonicTests

using Test
using Gridap

# Analytical manufactured solution
u(x) = x[1]*(x[1]-1)*x[2]*(x[2]-1)
f(x) = Δ(Δ(u))(x)
g(x) = Δ(u)(x)
@test f(VectorValue(0.5,0.5)) == 8.0
@test g(VectorValue(0.5,0.5)) == -1

# Domain
domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

# FE space
order = 2
V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order),dirichlet_tags="boundary")
U = TrialFESpace(V,u)

# Triangulation
Ω = Triangulation(model)
Γ = BoundaryTriangulation(model)
Λ = SkeletonTriangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)
dΛ = Measure(Λ,degree)
nΓ = get_normal_vector(Γ)
nΛ = get_normal_vector(Λ)

# Weak form
const h = (domain[2]-domain[1]) / partition[1]
const γ = 1
a(u,v) = ∫( Δ(u)*Δ(v) )dΩ +
         ∫( - mean(Δ(u))*jump(∇(v)⋅nΛ) - jump(∇(u)⋅nΛ)*mean(Δ(v)) + γ/h*jump(∇(u)⋅nΛ)*jump(∇(v)⋅nΛ) )dΛ
l(v) = ∫( v*f )dΩ + ∫( g*(∇(v)⋅nΓ) )dΓ
op = AffineFEOperator(a,l,U,V)

uₕ = solve(op)

# Error
e = u - uₕ
l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))
el2 = l2(e)
eh1 = h1(e)
tol = 1.0e-10
@test el2 < tol
@test eh1 < tol

end


text:
# Collection of the actual benchmarks that are run.

"""
Run the set of benchmarks for sequential and random array access for arrays that fit within
the DRAM cache in 2LM Mode.
"""
function run_2lm_fit_in_dram()
    @assert Threads.nthreads() == 24

    # Define closures to pass to to the server.
    measurements = [
        # Only really interested in the TLB misses
        #
        # Restrict counters to cores 25:48 - the cores we're actually using.
        () -> core_events([
            PCM.EventDescription(0x08, 0x0E, :dtlb_load_miss),
            PCM.EventDescription(0x49, 0x0E, :dtlb_store_miss),
        ], 25:48),
        # Chunk up `uncore` events into groups of 4.
        #
        # First collect data on tag statistics
        () -> uncore_events((
            tag_hit         = tagchk_hit(),
            tag_miss_clean  = tagchk_miss_clean(),
            tag_miss_dirty  = tagchk_miss_dirty(),
            dram_rq         = dram_rq(),
        )),
        # Bandwidth statistics
        () -> uncore_events((
            dram_reads  = cas_count_rd(),
            dram_writes = cas_count_wr(),
            pmm_reads   = pmm_read_cmd(),
            pmm_writes  = pmm_write_cmd(),
        )),
        # Queue Depths
        () -> uncore_events((
            unc_clocks  = unc_clocks(),
            pmm_rq      = pmm_rq(),
            pmm_wq      = pmm_wq(),
            dram_wq     = dram_wq(),
        ))
    ]

    params = KernelParams(
        29,     # Around 2G per threads (2^29 * 4 B). Total array size ~51 GB.
        1.0,    # Sample time: 1 second
        "2lm",  # Run this in 2LM
        10,     # 10 inner iterations - each iteration is pretty quick.
        measurements,
        [16],       # Only use AvX 512
        [true, false],    # Use standard loads/stores
        "mmap-1gb"
    )

    runkernels(
        params,
        sz -> hugepage_mmap(Float32, sz, Pagesize1G()),
    )
end

"""
Run the set of benchmarks for sequential and random array access for arrays that do not fit
within the DRAM cache in 2LM Mode.
"""
function run_2lm_exceeds_dram()
    @assert Threads.nthreads() == 24

    # Define closures to pass to to the server.
    measurements = [
        # Only really interested in the TLB misses
        #
        # Restrict counters to cores 25:48 - the cores we're actually using.
        () -> core_events([
            PCM.EventDescription(0x08, 0x0E, :dtlb_load_miss),
            PCM.EventDescription(0x49, 0x0E, :dtlb_store_miss),
        ], 25:48),
        # Chunk up `uncore` events into groups of 4.
        #
        # First collect data on tag statistics
        () -> uncore_events((
            tag_hit         = tagchk_hit(),
            tag_miss_clean  = tagchk_miss_clean(),
            tag_miss_dirty  = tagchk_miss_dirty(),
            dram_rq         = dram_rq(),
        )),
        # Bandwidth statistics
        () -> uncore_events((
            dram_reads  = cas_count_rd(),
            dram_writes = cas_count_wr(),
            pmm_reads   = pmm_read_cmd(),
            pmm_writes  = pmm_write_cmd(),
        )),
        # Queue Depths
        () -> uncore_events((
            unc_clocks  = unc_clocks(),
            pmm_rq      = pmm_rq(),
            pmm_wq      = pmm_wq(),
            dram_wq     = dram_wq(),
        ))
    ]

    params = KernelParams(
        32,     # Around 17G per threads (2^32 * 4 B). Total array size ~413 GB.
        1.0,    # Sample time: 1 second
        "2lm",  # Run this in 2LM
        2,      # 2 iterations - these take quite a while.
        measurements,
        [16],       # Only use AvX 512
        [false, true],
        "mmap-1gb"
    )

    runkernels(
        params,
        sz -> hugepage_mmap(Float32, sz, Pagesize1G()),
    )
end

function pmm_direct_test()
    # Get the number of threads.
    # Pick the array size based on the number of workers.
    nthreads = Threads.nthreads()

    # Minimum total array size of 30 GB
    min_total_size = 30E9 / sizeof(Float32)

    # Calculate a nice power of 2 that will get us there
    array_size = convert(Int, log2(nextpow(2, ceil(Int, min_total_size / nthreads))))

    measurements = [
        # Bandwidth statistics
        () -> uncore_events((
            dram_reads  = cas_count_rd(),
            dram_writes = cas_count_wr(),
            pmm_reads   = pmm_read_cmd(),
            pmm_writes  = pmm_write_cmd(),
        )),
        # Queue Depths
        () -> uncore_events((
            unc_clocks  = unc_clocks(),
            pmm_rq      = pmm_rq(),
            pmm_wq      = pmm_wq(),
            dram_wq     = dram_wq(),
        )),
    ]

    params = KernelParams(
        array_size,
        1.0,
        "1lm",
        5,
        measurements,
        [16],           # only use AvX 512
        [true, false],  # Use nontemporal and standard loads and stores
        "pmm",
    )

    # Run these experiments with a PersistentArray instead of a normal array.
    runkernels(
        params,
        sz -> PersistentArray{Float32}(undef, sz),
    )
end

function dram_direct_test(;test = false)
    @assert Threads.nthreads() == 24

    measurements = [
        # Bandwidth statistics
        () -> uncore_events((
            dram_reads  = cas_count_rd(),
            dram_writes = cas_count_wr(),
            pmm_reads   = pmm_read_cmd(),
            pmm_writes  = pmm_write_cmd(),
        )),
        # Queue Depths
        () -> uncore_events((
            unc_clocks  = unc_clocks(),
            pmm_rq      = pmm_rq(),
            dram_wq     = dram_wq(),
            dram_rq     = dram_rq(),
        )),
    ]

    # Set to a small size if we're doing a debug run.
    _size = test ? 20 : 29

    params = KernelParams(
        _size,
        1.0,
        "1lm",
        10,
        measurements,
        [16],           # only use AvX 512
        [true, false],  # Use nontemporal and standard loads and stores
        "dram",
    )

    # Run these experiments with a PersistentArray instead of a normal array.
    runkernels(
        params,
        sz -> hugepage_mmap(Float32, sz, Pagesize1G()),
    )
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
"""
    cross(fs::AbstractVector{<:GP})

Creates a multi-output GP from an `AbstractVector` of `GP`s.
"""
function LinearAlgebra.cross(fs::AbstractVector{<:AbstractGP})
    consistency_checks(fs)
    return CompositeGP((cross, fs), first(fs).gpc)
end

function consistency_checks(fs)
    @assert length(fs) >= 1
    @assert all([f.gpc == first(fs).gpc for f in fs])
end
ChainRulesCore.@non_differentiable consistency_checks(::Any)

const cross_args{T<:AbstractVector{<:AbstractGP}} = Tuple{typeof(cross), T}

function mean((_, fs)::cross_args, x::BlockData)
    blks = map((f, blk)->mean(f, blk), fs, blocks(x))
    return _collect(mortar(blks))
end

function cov((_, fs)::cross_args, x::BlockData)
    Cs = reshape(map((f, blk)->cov(f, (cross, fs), blk, x), fs, blocks(x)), :, 1)
    return _collect(mortar(reshape(Cs, :, 1)))
end

function var((_, fs)::cross_args, x::BlockData)
    cs = map(var, fs, blocks(x))
    return _collect(mortar(cs))
end

function cov((_, fs)::cross_args, x::BlockData, x′::BlockData)
    Cs = reshape(map((f, blk)->cov(f, (cross, fs), blk, x′), fs, blocks(x)), :, 1)
    return _collect(mortar(reshape(Cs, :, 1)))
end

function var((_, fs)::cross_args, x::BlockData, x′::BlockData)
    cs = map(var, fs, blocks(x), blocks(x′))
    return _collect(mortar(cs))
end

function cov((_, fs)::cross_args, f′::AbstractGP, x::BlockData, x′::AV)
    Cs = reshape(map((f, x)->cov(f, f′, x, x′), fs, blocks(x)), :, 1)
    return _collect(mortar(Cs))
end
function cov(f::AbstractGP, (_, fs)::cross_args, x::AV, x′::BlockData)
    Cs = reshape(map((f′, x′)->cov(f, f′, x, x′), fs, blocks(x′)), 1, :)
    return _collect(mortar(Cs))
end

function var(args::cross_args, f′::AbstractGP, x::BlockData, x′::AV)
    return diag(cov(args, f′, x, x′))
end
function var(f::AbstractGP, args::cross_args, x::AV, x′::BlockData)
    return diag(cov(f, args, x, x′))
end

#
# Build a single FiniteGP from a collection of FiniteGPs.
#

function finites_to_block(fs::AV{<:FiniteGP})
    return FiniteGP(
        cross(map(f->f.f, fs)),
        BlockData(map(f->f.x, fs)),
        make_block_noise(map(f->f.Σy, fs)),
    )
end

make_block_noise(Σys::Vector{<:Diagonal}) = Diagonal(Vector(mortar(diag.(Σys))))
make_block_noise(Σys::Vector) = block_diagonal(Σys)


text:
countprimes(stop::ConvT; tuplet::Int = 1) = countprimes(conv128(stop),tuplet)
#countprimes(stop::ConvT) = countprimes(conv128(stop)) ## no sense
function countprimes(stop; tuplet::Int = 1)
    if tuplet == 1
        _countprimes(stop)
    elseif tuplet == 2
        countprimes2(stop)
    elseif tuplet == 3
        countprimes3(stop)
    elseif tuplet == 4
        countprimes4(stop)
    elseif tuplet == 5
        countprimes5(stop)
    elseif tuplet == 6
        countprimes6(stop)
    else
        error("tuplet must be between 1 and 6")
    end
end

countprimes(start::ConvT, stop::ConvT; tuplet::Int = 1) = countprimes(conv128(start),conv128(stop),tuplet)
#countprimes(start::ConvT, stop::ConvT) = countprimes(conv128(start),conv128(stop))
function countprimes(start,stop; tuplet::Int = 1, alg::Symbol = :tabsieve )
    if tuplet == 1
        if alg == :next || stop > stoplimit
            countprimesb(start,stop)
        elseif alg == :nexta
            countprimesc(start,stop)
        elseif alg == :tabsieve
            _countprimes(start,stop)
        elseif alg == :sieve
            ntcountprimes(start,stop)
        else
            error("algorithm must be one of :sieve, :next, :nexta, :tabsieve")
        end
    elseif tuplet == 2
        countprimes2(start,stop)
    elseif tuplet == 3
        countprimes3(start,stop)
    elseif tuplet == 4
        countprimes4(start,stop)
    elseif tuplet == 5
        countprimes5(start,stop)
    elseif tuplet == 6
        countprimes6(start,stop)
    else
        error("tuplet must be between 1 and 6")
    end
end

scountprimes(stop::ConvT; tuplet::Int = 1) = scountprimes(conv128(stop),tuplet)
#scountprimes(stop::ConvT) = scountprimes(conv128(stop))
function scountprimes(stop; tuplet::Int = 1)
    if tuplet == 1
        _scountprimes(stop)
    elseif tuplet == 2
        scountprimes2(stop)
    elseif tuplet == 3
        scountprimes3(stop)
    elseif tuplet == 4
        scountprimes4(stop)
    elseif tuplet == 5
        scountprimes5(stop)
    elseif tuplet == 6
        scountprimes6(stop)
    else
        error("tuplet must be between 1 and 6")
    end
end

scountprimes(start::ConvT, stop::ConvT; tuplet::Int = 1) = scountprimes(conv128(start),conv128(stop),tuplet)
#scountprimes(start::ConvT, stop::ConvT) = scountprimes(conv128(start),conv128(stop))
function scountprimes(start,stop; tuplet::Int = 1)
    if tuplet == 1
        _scountprimes(start,stop)
    elseif tuplet == 2
        scountprimes2(start,stop)
    elseif tuplet == 3
        scountprimes3(start,stop)
    elseif tuplet == 4
        scountprimes4(start,stop)
    elseif tuplet == 5
        scountprimes5(start,stop)
    elseif tuplet == 6
        scountprimes6(start,stop)
    else
        error("tuplet must be between 1 and 6")
    end
end

printprimes(stop::ConvT; tuplet::Int = 1) = printprimes(conv128(stop),tuplet)
#printprimes(stop::ConvT) = printprimes(conv128(stop))
function printprimes(stop; tuplet::Int = 1)
    if tuplet == 1
        _printprimes(stop)
    elseif tuplet == 2
        printprimes2(stop)
    elseif tuplet == 3
        printprimes3(stop)
    elseif tuplet == 4
        printprimes4(stop)
    elseif tuplet == 5
        printprimes5(stop)
    elseif tuplet == 6
        printprimes6(stop)
    else
        error("tuplet must be between 1 and 6")
    end
end

printprimes(start::ConvT, stop::ConvT; tuplet::Int = 1) = printprimes(conv128(start),conv128(stop),tuplet)
#printprimes(start::ConvT, stop::ConvT) = printprimes(conv128(start),conv128(stop))
function printprimes(start,stop; tuplet::Int = 1)
    if tuplet == 1
        _printprimes(start,stop)
    elseif tuplet == 2
        printprimes2(start,stop)
    elseif tuplet == 3
        printprimes3(start,stop)
    elseif tuplet == 4
        printprimes4(start,stop)
    elseif tuplet == 5
        printprimes5(start,stop)
    elseif tuplet == 6
        printprimes6(start,stop)
    else
        error("tuplet must be between 1 and 6")
    end
end