text:
using Test
using HypothesisTests
using LittleGibbs

include("./testmodels.jl")

function test_gdemo_m(N, α₀, θ₀)
    # Testing consistency of conditional distribution with joint distribution,
    # see: Grosse, R. B. & Duvenaud, D. K. Testing MCMC Code. arXiv:1412.5218 [cs, stat] (2014).
    # p(m | λ, x)/p(m′ | λ, x) =!= p(m, λ, x)/p(m′, λ, x)

    (m, λ), x = gdemo_forward(N, α₀, θ₀)
    cond_m = gdemo_cond_m(α₀, θ₀, x, λ)
    m′ = m + randn()
    @test isapprox(logpdf(cond_m, m′) - logpdf(cond_m, m),
                   gdemo_logjoint(α₀, θ₀, x, m′, λ) - gdemo_logjoint(α₀, θ₀, x, m, λ),
                   rtol=0.1)
end

function test_gdemo_λ(N, α₀, θ₀)
    # Testing consistency of conditional distribution with joint distribution,
    # see: Grosse, R. B. & Duvenaud, D. K. Testing MCMC Code. arXiv:1412.5218 [cs, stat] (2014).
    # p(λ | m, x)/p(λ′ | m, x) =!= p(λ, m, x)/p(λ′, m, x)

    (m, λ), x = gdemo_forward(N, α₀, θ₀)
    cond_λ = gdemo_cond_λ(α₀, θ₀, x, m)
    λ′ = max(0, λ + randn())
    @test isapprox(logpdf(cond_λ, λ′) - logpdf(cond_λ, λ),
                   gdemo_logjoint(α₀, θ₀, x, m, λ′) - gdemo_logjoint(α₀, θ₀, x, m, λ),
                   rtol=0.2)
end

function test_gdemo(observations, α₀, θ₀, m_true, λ_true; m_init=0.0, λ_init=1.0)
    # Unit tests for conditionals
    test_gdemo_m(100, α₀, θ₀)
    test_gdemo_λ(100, α₀, θ₀)

    # Tests for expectations
    conditionals = gdemo(observations; α₀=α₀, θ₀=θ₀)

    model = BlockModel((m=conditionals.m, λ=(m,) -> Normal(λ_true, 0)))
    sampler = Gibbs((m=m_init, λ=λ_init))
    chain = sample(model, sampler, 10_000)
    chain_m = map(t -> t.params.m, chain[1_000:end])
    @test isapprox(mean(chain_m), m_true, atol=0.2, rtol=0.0)

    model = BlockModel((m=(λ,) -> Normal(m_true, 0), λ=conditionals.λ))
    sampler = Gibbs((m=m_init, λ=λ_init))
    chain = sample(model, sampler, 10_000)
    chain_λ = map(t -> t.params.λ, chain[1_000:end])
    @test isapprox(mean(chain_λ), λ_true, atol=0.2, rtol=0.0)

    model = BlockModel(conditionals)
    sampler = Gibbs((m=m_init, λ=λ_init))
    chain = sample(model, sampler, 10_000)
    chain_m = map(t -> t.params.m, chain[1_000:end])
    chain_λ = map(t -> t.params.λ, chain[1_000:end])
    @test isapprox(mean(chain_m), m_true, atol=0.2, rtol=0.0)
    @test isapprox(mean(chain_λ), λ_true, atol=0.2, rtol=0.0)
end

@testset "gdemo" begin
    N = 1000

    α₀, θ₀ = 2.0, inv(3.0)
    λ_true = rand(Gamma(α₀, θ₀))
    σ_true = √(1 / λ_true)
    m_true = rand(Normal(0, σ_true))
    test_gdemo(rand(Normal(m_true, σ_true), N), α₀, θ₀, m_true, λ_true; λ_init = α₀ * θ₀)

    α₀, θ₀ = 1.5, 1.5
    λ_true = rand(Gamma(α₀, θ₀))
    σ_true = √(1 / λ_true)
    m_true = rand(Normal(0, σ_true))
    test_gdemo(rand(Normal(m_true, σ_true), N), α₀, θ₀, m_true, λ_true; λ_init = α₀ * θ₀)

    # x = y = -2:0.05:2
    # @show exact_posterior(-1.5, -1.5)
    # display(contour(x, y, exact_posterior))
end

@testset "mixture" begin
    π = [0.5, 0.5]
    K = length(π)
    m = 0.5
    λ = 2.0
    σ = 0.1
    observations = [σ .* randn(10); 1 .+ σ .* randn(10)]
    N = length(observations)

    conditionals = mixture(π, K, m, λ, σ, observations)
    model = BlockModel(conditionals)
    sampler = Gibbs((z=rand((0, 1), N), μ=randn(2)))

    chain = sample(model, sampler, 11_000)
    chain_z = map(t -> t.params[:z], chain[1001:end])
    chain_μ = map(t -> t.params[:μ], chain[1001:end])

    ẑ = mean(chain_z)
    μ̂ = mean(chain_μ)
    z_true = [fill(1, 10); fill(2, 10)]
    μ_true = [0, 1]
    @info (ẑ=ẑ, μ̂=μ̂)
    @test isapprox(ẑ, z_true, atol=0.2, rtol=0.0) || isapprox(ẑ, z_true[[11:20; 1:10]], atol=0.2, rtol=0.0)
    @test isapprox(sort(μ̂), μ_true, atol=0.2, rtol=0.0)
end

# @testset "ising" begin
#     N = 20
#     β = 0.5
#     conditionals = isingdemo(β, N)
#     model = BlockModel(conditionals)
#     sampler = Gibbs(merge(NamedTuple(), [Symbol(:x, i) => 0 for i = 1:N]))
# end


text:
using LightGraphs

# assumes that genotype and phenotype graph have been built
function get_house_of_cards_null_model(graph_data::GGraphData)
    null_model = deepcopy(graph_data)

    weights = [e.weight for e in edges(null_model.phenotype_graph)]

    shuffle!(weights)

    for (i, edge) in edges(null_model.phenotype_graph) |> enumerate
        u, v = src(edge), dst(edge)

        add_edge!(null_model.phenotype_graph, u, v, weights[i])
    end

    return null_model
end

# assumes that genotype and phenotype graph have been built
function get_complete_house_of_cards_null_model(graph_data::GGraphData)
    null_model = deepcopy(graph_data)

    if isfile(CALCULATED_FOLDER * "phenotype_similarities")
        similarities = deserialize(CALCULATED_FOLDER * "phenotype_similarities")
    else
        similarities = _sample_phenotype_similarities(graph_data::GGraphData, 250)
        serialize(CALCULATED_FOLDER * "phenotype_similarities", similarities)
    end

    all_edges = edges(null_model.genotype_graph)
    sources = src.(all_edges)
    destinations = dst.(all_edges)
    weights::Vector{Real} = rand(similarities, length(sources))

    null_model.phenotype_graph = SimpleWeightedDiGraph(sources, destinations, weights)

    @info "$(nv(null_model.phenotype_graph)) $(ne(null_model.phenotype_graph))"

    return null_model
end

function _sample_phenotype_similarities(graph_data::GGraphData, n)
    similarities = []

    simple_genotype_graph = build_simple_graph(graph_data.genotype_graph)
    all_edges = unique(edges(simple_genotype_graph))

    for i in 1:n
        edge = rand(all_edges)
        u, v = src(edge), dst(edge)

        if has_edge(graph_data.phenotype_graph, u, v)
            similarity = graph_data.phenotype_graph.weights[v, u]
        else
            genotype_u, genotype_v = graph_data.genotype_vertex_mapping(u), graph_data.genotype_vertex_mapping(v)
            similarity = get_phenotype_similarity(genotype_u, genotype_v, 0.05, 10, 20)
        end

        push!(similarities, similarity)

        @info "$i / $n sampled ($similarity)."
    end

    return similarities
end

ext: jl

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using NCDatasets
using Test

fname = tempname()

for define_standard_name in [false,true]

    NCDataset(fname,"c") do ds
        # Dimensions

        ds.dim["xi_rho"] = 6
        ds.dim["xi_u"] = 5
        ds.dim["xi_v"] = 6
        ds.dim["xi_psi"] = 5
        ds.dim["eta_rho"] = 6
        ds.dim["eta_u"] = 6
        ds.dim["eta_v"] = 5
        ds.dim["eta_psi"] = 5
        ds.dim["N"] = 30
        ds.dim["s_rho"] = 30
        ds.dim["s_w"] = 31
        ds.dim["tracer"] = 63
        ds.dim["ocean_time"] = Inf # unlimited dimension

        # Declare variables

        nclon_rho = defVar(ds,"lon_rho", Float64, ("xi_rho", "eta_rho"))
        nclon_rho.attrib["long_name"] = "longitude of RHO-points"
        nclon_rho.attrib["units"] = "degree_east"
        nclon_rho.attrib["field"] = "lon_rho, scalar"
        if define_standard_name
            nclon_rho.attrib["standard_name"] = "longitude"
        end

        nclat_rho = defVar(ds,"lat_rho", Float64, ("xi_rho", "eta_rho"))
        nclat_rho.attrib["long_name"] = "latitude of RHO-points"
        nclat_rho.attrib["units"] = "degree_north"
        nclat_rho.attrib["field"] = "lat_rho, scalar"
        if define_standard_name
            nclon_rho.attrib["standard_name"] = "latitude"
        end

        nclon_u = defVar(ds,"lon_u", Float64, ("xi_u", "eta_u"))
        nclon_u.attrib["long_name"] = "longitude of U-points"
        nclon_u.attrib["units"] = "degree_east"
        nclon_u.attrib["field"] = "lon_u, scalar"

        nclat_u = defVar(ds,"lat_u", Float64, ("xi_u", "eta_u"))
        nclat_u.attrib["long_name"] = "latitude of U-points"
        nclat_u.attrib["units"] = "degree_north"
        nclat_u.attrib["field"] = "lat_u, scalar"

        nczeta = defVar(ds,"zeta", Float32, ("xi_rho", "eta_rho", "ocean_time"))
        nczeta.attrib["long_name"] = "free-surface"
        nczeta.attrib["units"] = "meter"
        nczeta.attrib["time"] = "ocean_time"
        nczeta.attrib["coordinates"] = "lat_rho lon_rho"
        nczeta.attrib["field"] = "free-surface, scalar, series"

        ncubar = defVar(ds,"ubar", Float32, ("xi_u", "eta_u", "ocean_time"))
        ncubar.attrib["long_name"] = "vertically integrated u-momentum component"
        ncubar.attrib["units"] = "meter second-1"
        ncubar.attrib["time"] = "ocean_time"
        ncubar.attrib["coordinates"] = "lat_u lon_u"
        ncubar.attrib["field"] = "ubar-velocity, scalar, series"
    end

    local ds
    ds = NCDataset(fname)
    @test name(coord(ds["zeta"],"longitude")) == "lon_rho"
    @test name(coord(ds["ubar"],"longitude")) == "lon_u"
    @test coord(ds["ubar"],"foobar") == nothing
    close(ds)
end


text:
int_rules_8_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.6*Gamma*functions =#
    @apply_utils Antiderivative(Gamma(~n, ~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * Gamma(~n, ~a + ~b * ~x)) / ~b - Gamma(~n + 1, ~a + ~b * ~x) / ~b <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative(Gamma(0, ~(b') * ~x) / ~x, ~x) => (~b * ~x * HypergeometricPFQ([1, 1, 1], [2, 2, 2], -(~b) * ~x) - EulerGamma * log(~x)) - (1 / 2) * log(~b * ~x) ^ 2 <-- FreeQ(~b, ~x)
    #= Antiderivative(Gamma(1,(~b')*(~x))/(~x),~x) := Antiderivative(1/(x*ℯ^(b*x)),x) <--  FreeQ(b,x) =#
    @apply_utils Antiderivative(Gamma(~n, ~(b') * ~x) / ~x, ~x) => -(Gamma(~n - 1, ~b * ~x)) + (~n - 1) * Antiderivative(Gamma(~n - 1, ~b * ~x) / ~x, ~x) <-- FreeQ(~b, ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative(Gamma(~n, ~(b') * ~x) / ~x, ~x) => Gamma(~n, ~b * ~x) / ~n + (1 / ~n) * Antiderivative(Gamma(~n + 1, ~b * ~x) / ~x, ~x) <-- FreeQ(~b, ~x) && ILtQ(~n, 0)
    @apply_utils Antiderivative(Gamma(~n, ~(b') * ~x) / ~x, ~x) => Gamma(~n) * log(~x) - ((~b * ~x) ^ ~n / (~n) ^ 2) * HypergeometricPFQ([~n, ~n], [1 + ~n, 1 + ~n], -(~b) * ~x) <-- FreeQ([~b, ~n], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * Gamma(~n, ~(b') * ~x), ~x) => ((~d * ~x) ^ (~m + 1) * Gamma(~n, ~b * ~x)) / (~d * (~m + 1)) - ((~d * ~x) ^ ~m * Gamma(~m + ~n + 1, ~b * ~x)) / (~b * (~m + 1) * (~b * ~x) ^ ~m) <-- FreeQ([~b, ~d, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~c + ~(d') * ~x) ^ ~(m') * Gamma(~n, ~a + ~(b') * ~x), ~x) => (1 / ~b) * Subst(Antiderivative(((~d * ~x) / ~b) ^ ~m * Gamma(~n, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && EqQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative(Gamma(~n, ~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x), ~x) => Antiderivative((~a + ~b * ~x) ^ (~n - 1) / ((~c + ~d * ~x) * ℯ ^ (~a + ~b * ~x)), ~x) + (~n - 1) * Antiderivative(Gamma(~n - 1, ~a + ~b * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * Gamma(~n, ~(a') + ~(b') * ~x), ~x) => Block([$UseGamma = True], ((~c + ~d * ~x) ^ (~m + 1) * Gamma(~n, ~a + ~b * ~x)) / (~d * (~m + 1)) + (~b / (~d * (~m + 1))) * Antiderivative(((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * ~x) ^ (~n - 1)) / ℯ ^ (~a + ~b * ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && ((IGtQ(~m, 0) || (IGtQ(~n, 0) || IntegersQ(~m, ~n))) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * Gamma(~n, ~(a') + ~(b') * ~x), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * Gamma(~n, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative(LogGamma(~(a') + ~(b') * ~x), ~x) => PolyGamma(-2, ~a + ~b * ~x) / ~b <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * LogGamma(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ ~m * PolyGamma(-2, ~a + ~b * ~x)) / ~b - ((~d * ~m) / ~b) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * PolyGamma(-2, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * LogGamma(~(a') + ~(b') * ~x), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * LogGamma(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x)
    @apply_utils Antiderivative(PolyGamma(~n, ~(a') + ~(b') * ~x), ~x) => PolyGamma(~n - 1, ~a + ~b * ~x) / ~b <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * PolyGamma(~n, ~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ ~m * PolyGamma(~n - 1, ~a + ~b * ~x)) / ~b - ((~d * ~m) / ~b) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * PolyGamma(~n - 1, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && GtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * PolyGamma(~n, ~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * PolyGamma(~n, ~a + ~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) * PolyGamma(~n + 1, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && LtQ(~m, -1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * PolyGamma(~n, ~(a') + ~(b') * ~x), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * PolyGamma(~n, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative(Gamma(~(a') + ~(b') * ~x) ^ ~(n') * PolyGamma(0, ~(a') + ~(b') * ~x), ~x) => Gamma(~a + ~b * ~x) ^ ~n / (~b * ~n) <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative(factorial(~(a') + ~(b') * ~x) ^ ~(n') * PolyGamma(0, ~(c') + ~(b') * ~x), ~x) => factorial(~a + ~b * ~x) ^ ~n / (~b * ~n) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && EqQ(~c, ~a + 1)
    @apply_utils Antiderivative(Gamma(~p, ~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * Gamma(~p, ~d * (~a + ~b * log(~c * (~x) ^ ~n))) + ~b * ~d * ~n * ℯ ^ (-(~a) * ~d) * Antiderivative((~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ (~p - 1) / (~c * (~x) ^ ~n) ^ (~b * ~d), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative(Gamma(~p, ~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))) / ~x, ~x) => (1 / ~n) * Subst(Gamma(~p, ~d * (~a + ~b * ~x)), ~x, log(~c * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * Gamma(~p, ~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * Gamma(~p, ~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) + ((~b * ~d * ~n * ℯ ^ (-(~a) * ~d) * (~e * ~x) ^ (~b * ~d * ~n)) / ((~m + 1) * (~c * (~x) ^ ~n) ^ (~b * ~d))) * Antiderivative((~e * ~x) ^ (~m - ~b * ~d * ~n) * (~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(Gamma(~p, ~(f') * (~(a') + ~(b') * log(~(c') * (~d + ~(e') * ~x) ^ ~(n')))), ~x) => (1 / ~e) * Subst(Antiderivative(Gamma(~p, ~f * (~a + ~b * log(~c * (~x) ^ ~n))), ~x), ~x, ~d + ~e * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x)
    @apply_utils Antiderivative((~g + ~(h') * ~x) ^ ~(m') * Gamma(~p, ~(f') * (~(a') + ~(b') * log(~(c') * (~d + ~(e') * ~x) ^ ~(n')))), ~x) => (1 / ~e) * Subst(Antiderivative(((~g * ~x) / ~d) ^ ~m * Gamma(~p, ~f * (~a + ~b * log(~c * (~x) ^ ~n))), ~x), ~x, ~d + ~e * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p], ~x) && EqQ(~e * ~g - ~d * ~h, 0)
end

max_stars_repo_path: src/rules/8 Special functions/8.6 Gamma functions.jl


text:
using Documenter
using VectorSphericalHarmonics

DocMeta.setdocmeta!(VectorSphericalHarmonics, :DocTestSetup, :(using VectorSphericalHarmonics); recursive=true)

makedocs(;
    modules=[VectorSphericalHarmonics],
    authors="Jishnu Bhattacharya",
    repo="https://github.com/jishnub/VectorSphericalHarmonics.jl/blob/{commit}{path}#L{line}",
    sitename="VectorSphericalHarmonics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jishnub.github.io/VectorSphericalHarmonics.jl",
        assets=String[],
    ),
    pages=[
        "Reference" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jishnub/VectorSphericalHarmonics.jl",
)


text:
abstract type AbstractReMat{T} <: AbstractMatrix{T} end

"""
    ReMat{T,S} <: AbstractMatrix{T}

A section of a model matrix generated by a random-effects term.

# Fields
- `trm`: the grouping factor as a `StatsModels.CategoricalTerm`
- `refs`: indices into the levels of the grouping factor as a `Vector{Int32}`
- `levels`: the levels of the grouping factor
- `z`: transpose of the model matrix generated by the left-hand side of the term
- `wtz`: a weighted copy of `z` (`z` and `wtz` are the same object for unweighted cases)
- `λ`: a `LowerTriangular` matrix of size `S×S`
- `inds`: a `Vector{Int}` of linear indices of the potential nonzeros in `λ`
- `adjA`: the adjoint of the matrix as a `SparseMatrixCSC{T}`
"""
mutable struct ReMat{T,S} <: AbstractReMat{T}
    trm
    refs::Vector{Int32}
    levels
    cnames::Vector{String}
    z::Matrix{T}
    wtz::Matrix{T}
    λ::LowerTriangular{T,Matrix{T}}
    inds::Vector{Int}
    adjA::SparseMatrixCSC{T,Int32}
    scratch::Matrix{T}
end

"""
    amalgamate(reterms::Vector{AbstractReMat})

Combine multiple ReMat with the same grouping variable into a single object.
"""
amalgamate(reterms::Vector{AbstractReMat{T}}) where {T} = _amalgamate(reterms,T)

function _amalgamate(reterms::Vector, T::Type)
    factordict = Dict{Symbol, Vector{Int}}()
    for (i, rt) in enumerate(reterms)
        push!(get!(factordict, fname(rt), Int[]), i)
    end
    length(factordict) == length(reterms) && return reterms
    value = AbstractReMat{T}[]
    for (f, inds) in factordict
        if isone(length(inds))
            push!(value, reterms[only(inds)])
        else
            trms = reterms[inds]
            trm1 = first(trms)
            trm = trm1.trm
            refs = refarray(trm1)
            levs = trm1.levels
            cnames =  foldl(vcat, rr.cnames for rr in trms)
            z = foldl(vcat, rr.z for rr in trms)

            Snew = size(z, 1)
            btemp = Matrix{Bool}(I, Snew, Snew)
            offset = 0
            for m in indmat.(trms)
                sz = size(m, 1)
                inds = (offset + 1):(offset + sz)
                view(btemp, inds, inds) .= m
                offset += sz
            end
            inds = (1:abs2(Snew))[vec(btemp)]

            λ = LowerTriangular(Matrix{T}(I, Snew, Snew))
            scratch =  foldl(vcat, rr.scratch for rr in trms)

            push!(
                value,
                ReMat{T,Snew}(trm, refs, levs, cnames, z, z, λ, inds, adjA(refs,z), scratch)
            )
        end
    end
    value
end

"""
    adjA(refs::AbstractVector, z::AbstractMatrix{T})

Returns the adjoint of an `ReMat` as a `SparseMatrixCSC{T,Int32}`
"""
function adjA(refs::AbstractVector, z::AbstractMatrix)
    S, n = size(z)
    length(refs) == n || throw(DimensionMismatch)
    J = Int32.(1:n)
    II = refs
    if S > 1
        J = repeat(J, inner=S)
        II = Int32.(vec([(r - 1)*S + j for j in 1:S, r in refs]))
    end
    sparse(II, J, vec(z))
end

Base.size(A::ReMat) = (length(A.refs), length(A.scratch))

SparseArrays.sparse(A::ReMat) = adjoint(A.adjA)

Base.getindex(A::ReMat, i::Integer, j::Integer) = getindex(A.adjA, j, i)

"""
    nranef(A::ReMat)

Return the number of random effects represented by `A`.  Zero unless `A` is an `ReMat`.
"""
nranef(A::ReMat) = size(A.adjA, 1)

LinearAlgebra.cond(A::ReMat) = cond(A.λ)

"""
    fname(A::ReMat)

Return the name of the grouping factor as a `Symbol`
"""
fname(A::ReMat) = fname(A.trm)
fname(A::CategoricalTerm) = A.sym
fname(A::InteractionTerm) = Symbol(join(fname.(A.terms), " & "))

getθ(A::ReMat{T}) where {T} = getθ!(Vector{T}(undef, nθ(A)), A)

"""
    getθ!(v::AbstractVector{T}, A::ReMat{T}) where {T}

Overwrite `v` with the elements of the blocks in the lower triangle of `A.Λ` (column-major ordering)
"""
function getθ!(v::AbstractVector{T}, A::ReMat{T}) where {T}
    length(v) == length(A.inds) || throw(DimensionMismatch("length(v) ≠ length(A.inds)"))
    m = A.λ.data
    @inbounds for (j, ind) in enumerate(A.inds)
        v[j] = m[ind]
    end
    v
end

function DataAPI.levels(A::ReMat)
    # These checks are for cases where unused levels are present.
    # Such cases may never occur b/c of the way an ReMat is constructed.
    pool = A.levels
    present = falses(size(pool))
    @inbounds for i in A.refs
        present[i] = true
        all(present) && return pool
    end
    pool[present]
end

"""
    indmat(A::ReMat)

Return a `Bool` indicator matrix of the potential non-zeros in `A.λ`
"""
function indmat end

indmat(rt::ReMat{T,1}) where {T} = ones(Bool, 1, 1)
indmat(rt::ReMat{T,S}) where {T,S} = reshape([i in rt.inds for i in 1:abs2(S)], S, S)

nlevs(A::ReMat) = length(A.levels)

"""
    nθ(A::ReMat)

Return the number of free parameters in the relative covariance matrix λ
"""
nθ(A::ReMat) = length(A.inds)

"""
    lowerbd{T}(A::ReMat{T})

Return the vector of lower bounds on the parameters, `θ` associated with `A`

These are the elements in the lower triangle of `A.λ` in column-major ordering.
Diagonals have a lower bound of `0`.  Off-diagonals have a lower-bound of `-Inf`.
"""
lowerbd(A::ReMat{T}) where {T} =
    T[x ∈ diagind(A.λ.data) ? zero(T) : T(-Inf) for x in A.inds]

"""
    isnested(A::ReMat, B::ReMat)

Is the grouping factor for `A` nested in the grouping factor for `B`?

That is, does each value of `A` occur with just one value of B?
"""
function isnested(A::ReMat, B::ReMat)
    size(A, 1) == size(B, 1) || throw(DimensionMismatch("must have size(A,1) == size(B,1)"))
    bins = zeros(Int32, nlevs(A))
    @inbounds for (a, b) in zip(A.refs, B.refs)
        bba = bins[a]
        if iszero(bba)    # bins[a] not yet set?
            bins[a] = b   # set it
        elseif bba ≠ b    # set to another value?
            return false
        end
    end
    true
end

function lmulΛ!(adjA::Adjoint{T,ReMat{T,1}}, B::Matrix{T}) where {T}
    lmul!(only(adjA.parent.λ.data), B)
end

function lmulΛ!(adjA::Adjoint{T,ReMat{T,1}}, B::SparseMatrixCSC{T}) where {T}
    lmul!(only(adjA.parent.λ.data), nonzeros(B))
    B
end

function lmulΛ!(adjA::Adjoint{T,ReMat{T,1}}, B::M) where{M<:AbstractMatrix{T}} where {T}
    lmul!(only(adjA.parent.λ.data), B)
end

function lmulΛ!(adjA::Adjoint{T,ReMat{T,S}}, B::VecOrMat{T}) where {T,S}
    lmul!(adjoint(adjA.parent.λ), reshape(B, S, :))
    B
end

function lmulΛ!(adjA::Adjoint{T,<:ReMat{T,S}}, B::BlockedSparse{T}) where {T,S}
    lmulΛ!(adjA, nonzeros(B.cscmat))
    B
end

function lmulΛ!(adjA::Adjoint{T,ReMat{T,1}}, B::BlockedSparse{T,1,P}) where {T,P}
    lmul!(only(adjA.parent.λ.data), nonzeros(B.cscmat))
    B
end

function lmulΛ!(adjA::Adjoint{T,<:ReMat{T,S}}, B::SparseMatrixCSC{T}) where {T,S}
    lmulΛ!(adjA, nonzeros(B))
    B
end

LinearAlgebra.Matrix(A::ReMat) = Matrix(sparse(A))

function LinearAlgebra.mul!(C::Diagonal{T}, adjA::Adjoint{T,<:ReMat{T,1}},
        B::ReMat{T,1}) where {T}
    A = adjA.parent
    @assert A === B
    d = C.diag
    fill!(d, zero(T))
    @inbounds for (ri, Azi) in zip(A.refs, A.wtz)
        d[ri] += abs2(Azi)
    end
    C
end

function *(adjA::Adjoint{T,<:ReMat{T,1}}, B::ReMat{T,1}) where {T}
    A = adjA.parent
    A === B ? mul!(Diagonal(Vector{T}(undef, size(B, 2))), adjA, B) :
    sparse(Int32.(A.refs), Int32.(B.refs), vec(A.wtz .* B.wtz))
end

*(adjA::Adjoint{T,<:ReMat{T}}, B::ReMat{T}) where {T} = adjA.parent.adjA * sparse(B)
*(adjA::Adjoint{T,<:FeMat{T}}, B::ReMat{T}) where {T} =
    mul!(Matrix{T}(undef, rank(adjA.parent), size(B, 2)), adjA, B)

function LinearAlgebra.mul!(C::Matrix{T}, adjA::Adjoint{T,<:FeMat{T}}, B::ReMat{T,1},
        α::Number, β::Number) where {T}
    A = adjA.parent
    Awt = A.wtx
    n, p = size(Awt)
    r = A.rank
    m, q = size(B)
    size(C) == (r, q) && m == n || throw(DimensionMismatch())
    isone(β) || rmul!(C, β)
    zz = B.wtz
    @inbounds for (j, rrj) in enumerate(B.refs)
        αzj = α * zz[j]
        for i in 1:r
            C[i, rrj] += αzj * Awt[j, i]
        end
    end
    C
end

function LinearAlgebra.mul!(C::Matrix{T}, adjA::Adjoint{T,<:FeMat{T}}, B::ReMat{T,S},
    α::Number, β::Number) where {T,S}
    A = adjA.parent
    Awt = A.wtx
    r = rank(A)
    rr = B.refs
    scr = B.scratch
    vscr = vec(scr)
    Bwt = B.wtz
    n = length(rr)
    q = length(scr)
    size(C) == (r, q) && size(Awt, 1) == n || throw(DimensionMismatch(""))
    isone(β) || rmul!(C, β)
    @inbounds for i in 1:r
        fill!(scr, 0)
        for k in 1:n
            aki = α * Awt[k,i]
            kk = Int(rr[k])
            for ii in 1:S
                scr[ii, kk] += aki * Bwt[ii, k]
            end
        end
        for j in 1:q
            C[i, j] += vscr[j]
        end
    end
    C
end

function LinearAlgebra.mul!(C::SparseMatrixCSC{T}, adjA::Adjoint{T,<:ReMat{T,1}},
        B::ReMat{T,1}) where {T}
    A = adjA.parent
    m, n = size(B)
    size(C, 1) == size(A, 2) && n == size(C, 2) && size(A, 1) == m || throw(DimensionMismatch)
    Ar = A.refs
    Br = B.refs
    Az = A.wtz
    Bz = B.wtz
    nz = nonzeros(C)
    rv = rowvals(C)
    fill!(nz, zero(T))
    for k in 1:m       # iterate over rows of A and B
        i = Ar[k]      # [i,j] are Cartesian indices in C - find and verify corresponding position K in rv and nz
        j = Br[k]
        coljlast = Int(C.colptr[j + 1] - 1)
        K = searchsortedfirst(rv, i, Int(C.colptr[j]), coljlast, Base.Order.Forward)
        if K ≤ coljlast && rv[K] == i
            nz[K] += Az[k] * Bz[k]
        else
            throw(ArgumentError("C does not have the nonzero pattern of A'B"))
        end
    end
    C
end

function LinearAlgebra.mul!(C::UniformBlockDiagonal{T}, adjA::Adjoint{T,ReMat{T,S}},
        B::ReMat{T,S}) where {T,S}
    A = adjA.parent
    @assert A === B
    Cd = C.data
    size(Cd) == (S, S, nlevs(B)) || throw(DimensionMismatch(""))
    fill!(Cd, zero(T))
    Awtz = A.wtz
    for (j, r) in enumerate(A.refs)
        @inbounds for i in 1:S
            zij = Awtz[i,j]
            for k in 1:S
                Cd[k, i, r] += zij * Awtz[k,j]
            end
        end
    end
    C
end

function LinearAlgebra.mul!(C::Matrix{T}, adjA::Adjoint{T,ReMat{T,S}},
        B::ReMat{T,P}) where {T,S,P}
    A = adjA.parent
    m, n = size(A)
    p, q = size(B)
    m == p && size(C, 1) == n && size(C, 2) == q || throw(DimensionMismatch(""))
    fill!(C, zero(T))

    Ar = A.refs
    Br = B.refs
    if isone(S) && isone(P)
        for (ar, az, br, bz) in zip(Ar, vec(A.wtz), Br, vec(B.wtz))
            C[ar, br] += az * bz
        end
        return C
    end
    ab = S * P
    Az = A.wtz
    Bz = B.wtz
    for i in 1:m
        Ari = Ar[i]
        Bri = Br[i]
        ioffset = (Ari - 1) * S
        joffset = (Bri - 1) * P
        for jj in 1:P
            jjo = jj + joffset
            Bzijj = Bz[jj, i]
            for ii in 1:S
                C[ii + ioffset, jjo] += Az[ii, i] * Bzijj
            end
        end
    end
    C
end

function *(adjA::Adjoint{T,<:ReMat{T,S}}, B::ReMat{T,P}) where {T,S,P}
    A = adjA.parent
    if A === B
        return mul!(UniformBlockDiagonal(Array{T}(undef, S, S, nlevs(A))), adjA, A)
    end
    cscmat = A.adjA * adjoint(B.adjA)
    if nnz(cscmat) > *(0.25, size(cscmat)...)
        return Matrix(cscmat)
    end

    BlockedSparse{T,S,P}(cscmat, reshape(cscmat.nzval, S, :), cscmat.colptr[1:P:(cscmat.n + 1)])
end

function PCA(A::ReMat{T,1}; corr::Bool=true) where {T}
    val = ones(T, 1, 1)
    # TODO: use DataAPI
    PCA(corr ? val : abs(only(A.λ)) * val, A.cnames; corr=corr)
end

# TODO: use DataAPI
PCA(A::ReMat{T,S}; corr::Bool=true) where {T,S} = PCA(A.λ, A.cnames; corr=corr)

refarray(A::ReMat) = A.refs

refpool(A::ReMat) = A.levels

refvalue(A::ReMat, i::Integer) = A.levels[i]

function reweight!(A::ReMat, sqrtwts::Vector)
    if length(sqrtwts) > 0
        if A.z === A.wtz
            A.wtz = similar(A.z)
        end
        rmul!(copyto!(A.wtz, A.z), Diagonal(sqrtwts))
    end
    A
end

rmulΛ!(A::Matrix{T}, B::ReMat{T,1}) where{T} = rmul!(A, only(B.λ.data))

function rmulΛ!(A::SparseMatrixCSC{T}, B::ReMat{T,1}) where {T}
    rmul!(nonzeros(A), only(B.λ.data))
    A
end

function rmulΛ!(A::Matrix{T}, B::ReMat{T,S}) where {T,S}
    m, n = size(A)
    q, r = divrem(n, S)
    iszero(r) || throw(DimensionMismatch("size(A, 2) is not a multiple of block size"))
    λ = B.λ
    for k = 1:q
        coloffset = (k - 1) * S
        rmul!(view(A, :, coloffset+1:coloffset+S), λ)
    end
    A
end

function rmulΛ!(A::BlockedSparse{T,S,P}, B::ReMat{T,P}) where {T,S,P}
    cbpt = A.colblkptr
    csc = A.cscmat
    nzv = csc.nzval
    for j in 1:div(csc.n, P)
        rmul!(reshape(view(nzv, cbpt[j]:(cbpt[j + 1] - 1)), :, P), B.λ)
    end
    A
end

rowlengths(A::ReMat{T,1}) where {T} = vec(abs.(A.λ.data))

function rowlengths(A::ReMat)
    ld = A.λ.data
    [norm(view(ld, i, 1:i)) for i in 1:size(ld, 1)]
end

"""
    scaleinflate!(L::AbstractMatrix, Λ::ReMat)

Overwrite L with `Λ'LΛ + I`
"""
function scaleinflate! end

function scaleinflate!(Ljj::Diagonal{T}, Λj::ReMat{T,1}) where {T}
    Ljjd = Ljj.diag
    broadcast!((x, λsqr) -> x * λsqr + 1, Ljjd, Ljjd, abs2(only(Λj.λ.data)))
    Ljj
end

function scaleinflate!(Ljj::Matrix{T}, Λj::ReMat{T,1}) where {T}
    lambsq = abs2(only(Λj.λ.data))
    @inbounds for i in diagind(Ljj)
        Ljj[i] *= lambsq
        Ljj[i] += one(T)
    end
    Ljj
end

function scaleinflate!(Ljj::UniformBlockDiagonal{T}, Λj::ReMat{T,S}) where {T,S}
    λ = Λj.λ
    dind = diagind(S, S)
    Ldat = Ljj.data
    for k in axes(Ldat, 3)
        f = view(Ldat, :, :, k)
        lmul!(λ', rmul!(f, λ))
        for i in dind
            f[i] += one(T)  # inflate diagonal
        end
    end
    Ljj
end

function scaleinflate!(Ljj::Matrix{T}, Λj::ReMat{T,S}) where{T,S}
    n = LinearAlgebra.checksquare(Ljj)
    q, r = divrem(n, S)
    iszero(r) || throw(DimensionMismatch("size(Ljj, 1) is not a multiple of S"))
    λ = Λj.λ
    offset = 0
    @inbounds for k in 1:q
        inds = (offset + 1):(offset + S)
        tmp = view(Ljj, inds, inds)
        lmul!(adjoint(λ), rmul!(tmp, λ))
        offset += S
    end
    for k in diagind(Ljj)
        Ljj[k] += 1
    end
    Ljj
end

function setθ!(A::ReMat{T}, v::AbstractVector{T}) where {T}
    A.λ.data[A.inds] = v
    A
end

function σs(A::ReMat{T,1}, sc::T) where {T}
    NamedTuple{(Symbol(only(A.cnames)),)}(sc*abs(only(A.λ.data)),)
end

function σs(A::ReMat{T}, sc::T) where {T}
    λ = A.λ.data
    NamedTuple{(Symbol.(A.cnames)...,)}(ntuple(i -> sc*norm(view(λ,i,1:i)), size(λ, 1)))
end

function σρs(A::ReMat{T,1}, sc::T) where {T}
    NamedTuple{(:σ,:ρ)}(
        (NamedTuple{(Symbol(only(A.cnames)),)}((sc*abs(only(A.λ.data)),)), ())
    )
end

function ρ(i, λ::AbstractMatrix{T}, im::Matrix{Bool}, indpairs, σs, sc::T)::T where {T}
    row, col = indpairs[i]
    if iszero(dot(view(im, row, :), view(im, col, :)))
        -zero(T)
    else
        dot(view(λ, row, :), view(λ, col, :)) * abs2(sc) / (σs[row] * σs[col])
    end
end

function σρs(A::ReMat{T}, sc::T) where {T}
    λ = A.λ.data
    k = size(λ, 1)
    im = indmat(A)
    indpairs = checkindprsk(k)
    σs = NamedTuple{(Symbol.(A.cnames)...,)}(ntuple(i -> sc*norm(view(λ,i,1:i)), k))
    NamedTuple{(:σ,:ρ)}((σs, ntuple(i -> ρ(i,λ,im,indpairs,σs,sc), (k * (k - 1)) >> 1)))
end

"""
    corrmat(A::ReMat)

Return the estimated correlation matrix for `A`.  The diagonal elements are 1
and the off-diagonal elements are the correlations between those random effect
terms

# Example

Note that trailing digits may vary slightly depending on the local platform.

```julia-repl
julia> using MixedModels

julia> mod = fit(MixedModel,
                 @formula(rt_trunc ~ 1 + spkr + prec + load + (1 + spkr + prec | subj)),
                 MixedModels.dataset(:kb07));

julia> VarCorr(mod)
Variance components:
             Column      Variance  Std.Dev.  Corr.
subj     (Intercept)     136591.782 369.583
         spkr: old        22922.871 151.403 +0.21
         prec: maintain   32348.269 179.856 -0.98 -0.03
Residual                 642324.531 801.452

julia> MixedModels.corrmat(mod.reterms[1])
3×3 LinearAlgebra.Symmetric{Float64,Array{Float64,2}}:
  1.0        0.214816   -0.982948
  0.214816   1.0        -0.0315607
 -0.982948  -0.0315607   1.0
```
"""
function corrmat(A::ReMat{T}) where {T}
    λ = A.λ
    λnorm = rownormalize!(copy!(zeros(T, size(λ)), λ))
    Symmetric(λnorm * λnorm', :L)
end

vsize(A::ReMat{T,S}) where {T,S} = S

function zerocorr!(A::ReMat{T}) where {T}
    λ = A.λ
    # zero out all entries not on the diagonal
    λ[setdiff(A.inds, diagind(λ))] .= 0
    A.inds = intersect(A.inds, diagind(λ))
    A
end


text:
module Array2file

export writefile, readfile

"""
    writefile(A::Matrix{<:Number}, str::String)

Write a matrix on disc.

This function will write the matrix A in the file located in str.
Be aware that the file will write the colums of the matrix in the
rows of the data file, and the other way around. This is due to the
column majorness of julia, but the row majorness of a text file.

# Example
```jldoctest
julia> A = reshape([1, 2, 3, 4], 2, 2)
2×2 Array{Int64,2}:
 1  3
 2  4

julia> writefile(A, "test.txt")

shell> cat test.txt
1	2
3	4
```

See also: [`readfile`](@ref)
"""
function writefile(A::Matrix{<:Number}, str::String)
    # A     ...   the array that will be saved on disc
    # str   ...   the location of the text file

    # open the output stream
    open(str, "w") do io

        # iterate over the elements of the matrix
        for i=1:size(A, 2)
            for j=1:size(A, 1)

                # write the corresponding element on the file
                # and seperate the value with an tabulator
                write(io, string(A[j, i])*"\t")
            end

            #after each written column write a newline
            write(io, "\n")
        end

    # close the io stream
    end

    # return nothing
    return
end

"""
    writefile(A::Vector{<:Number}, str::String)

Write a vector to a file.

This function will write the Vector A to the file
located at str.

# Example
```jldoctest
julia> writefile([1, 2, 3, 4], "test.txt")

shell> cat test.txt
1	2	3	4
```

See also: [`readfile`](@ref)
"""
function writefile(A::Vector{<:Number}, str::String)
    # A     ...   the vector that will be saved on disc
    # str   ...   the location of the text file

    # open the output stream
    open(str, "w") do io

        # iterate over the elements of the vector
        for i=1:size(A, 1)

            # write the elements on the file, add separator
            write(io, string(A[i])*"\t")
        end

        # signal the end of the file
        write(io, "\n")

    # close the file stream
    end

    # return nothing
    return
end

"""
    readfile(str::String, dtype::DataType)

Read in a matrix.

This function reads in a vector or a matrix that is save at the str
position. The input will be parsed into the DataType that is given
in the dtype variable.

# Example
```jldoctest
shell> cat test.txt
1	2
3	4

julia> readfile("test.txt", Int)
2×2 Array{Int64,2}:
 1  3
 2  4
```

See also: [`writefile`](@ref)
"""
function readfile(str::String, dtype::DataType)
    # str     ...   the location of the file that contains the matrix
    # dtype   ...   the datatype that will be parsed

    # load the data into the input string
    input = open(f->read(f, String), str)

    # set the starting index
    index_start = 1

    # start the row counter
    rows = 0

    # create the starting Array
    output = Array{dtype}(undef, 0)

    # create the endless loop that evaluates the input stream data
    while true

        # read out the position of the next tabulator or newline command
        index_tab = findnext(isequal('\t'), input, index_start)
        index_nl  = findnext(isequal('\n'), input, index_start)

        # if the matrix is already read in, break the loop
        if index_tab == nothing
            rows += 1
            break

        # if the next seperator is a tab, parse the chars to this point,
        # and fill the output array with it. Set a new index_start
        elseif index_tab < index_nl
            index_end = index_tab
            append!(output, parse(dtype, input[index_start:index_end-1]))
            index_start = index_end + 1

        # if the next seperator is a newline, increase the row counter,
        # and set the new start index
        else
            rows += 1
            index_start += 1
        end
    end

    # reshape the array to the wanted outputsize
    output = reshape(output, round(Int, length(output)/rows), rows)

    return output
end

# end module
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iupim.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const libiupim_ = "iupim") : (const libiupim_ = "libiupim")  # Name of shared lib.

const Ihandle = Cvoid

function IupLoadImage(file_name::String)
    ccall((:IupLoadImage, libiupim_), Ptr{Ihandle}, (Ptr{UInt8},), file_name)
end

function IupSaveImage(ih::Ptr{Ihandle}, file_name::String, format::String)
    ccall((:IupSaveImage, libiupim_), Cint, (Ptr{Ihandle}, Ptr{UInt8}, Ptr{UInt8}), ih, file_name, format)
end

function IupGetNativeHandleImage(handle::Ptr{Cvoid})
    ccall((:IupGetNativeHandleImage, libiupim_), Ptr{imImage}, (Ptr{Cvoid},), handle)
end

function IupGetImageNativeHandle(image::Ptr{imImage})
    ccall((:IupGetImageNativeHandle, libiupim_), Ptr{Cvoid}, (Ptr{imImage},), image)
end

function IupImageFromImImage(image::Ptr{imImage})
    ccall((:IupImageFromImImage, libiupim_), Ptr{Ihandle}, (Ptr{imImage},), image)
end


text:
"""To implement a sub-type of this you need to implement the methods below."""
abstract type Results end
function get_base_power(r::T) where {T <: Results}
    error("get_base_power must be implemented for $T")
end

function get_variables(r::T) where {T <: Results}
    error("get_variables must be implemented for $T")
end

function get_parameters(r::T) where {T <: Results}
    error("get_parameters must be implemented for $T")
end

function get_total_cost(r::T) where {T <: Results}
    error("get_total_cost must be implemented for $T")
end

function get_optimizer_log(r::T) where {T <: Results}
    error("get_optimizer_log must be implemented for $T")
end

function get_time_stamp(r::T) where {T <: Results}
    error("get_time_stamp must be implemented for $T")
end

function write_results(r::T) where {T <: Results}
    error("write_results must be implemented for $T")
end


text:
"""
    save_at(worker, sym, val)

Saves value `val` to symbol `sym` at `worker`. `sym` should be quoted (or
contain a symbol). `val` gets unquoted in the processing and evaluated at the
worker, quote it if you want to pass exact command to the worker.

This is loosely based on the package ParallelDataTransfers, but made slightly
more flexible by omitting/delaying the explicit fetches etc. In particular,
`save_at` is roughly the same as `ParallelDataTransfers.sendto`, and
`get_val_from` works very much like `ParallelDataTransfers.getfrom`.

# Return value

A future with Nothing that can be fetched to see that the operation has
finished.

# Examples

    addprocs(1)
    save_at(2,:x,123)       # saves 123
    save_at(2,:x,myid())    # saves 1
    save_at(2,:x,:(myid())) # saves 2
    save_at(2,:x,:(:x))     # saves the symbol :x
                            # (just :x won't work because of unquoting)

# Note: Symbol scope

The symbols are saved in Main module on the corresponding worker. For example,
`save_at(1, :x, nothing)` _will_ erase your local `x` variable. Beware of name
collisions.
"""
function save_at(worker, sym::Symbol, val; mod = Main)
    remotecall(() -> Base.eval(mod, :(
        begin
            $sym = $val
            nothing
        end
    )), worker)
end

"""
    get_from(worker,val)

Get a value `val` from a remote `worker`; quoting of `val` works just as with
`save_at`. Returns a future with the requested value.
"""
function get_from(worker, val; mod = Main)
    remotecall(() -> Base.eval(mod, val), worker)
end

"""
    get_val_from(worker,val)

Shortcut for instantly fetching the future from `get_from`.
"""
function get_val_from(worker, val)
    fetch(get_from(worker, val))
end

"""
    remove_from(worker,sym)

Sets symbol `sym` on `worker` to `nothing`, effectively freeing the data.
"""
function remove_from(worker, sym::Symbol)
    save_at(worker, sym, nothing)
end

"""
    scatter_array(sym, x::Array, pids; dim=1)::Dinfo

Distribute roughly equal parts of array `x` separated on dimension `dim` among
`pids` into a worker-local variable `sym`.

Returns the `Dinfo` structure for the distributed data.
"""
function scatter_array(sym::Symbol, x::Array, pids; dim = 1)::Dinfo
    n = length(pids)
    dims = size(x)

    for f in [
        begin
            extent = [(1:s) for s in dims]
            extent[dim] = (1+div((wid - 1) * dims[dim], n)):div(wid * dims[dim], n)
            save_at(pid, sym, x[extent...])
        end for (wid, pid) in enumerate(pids)
    ]
        fetch(f)
    end

    return Dinfo(sym, pids)
end

"""
    unscatter(sym, workers)

Remove the loaded data from workers.
"""
function unscatter(sym::Symbol, workers)
    for f in [remove_from(pid, sym) for pid in workers]
        fetch(f)
    end
end

"""
    unscatter(dInfo::Dinfo)

Remove the loaded data described by `dInfo` from the corresponding workers.
"""
function unscatter(dInfo::Dinfo)
    unscatter(dInfo.val, dInfo.workers)
end

"""
    dexec(val, fn, workers)

Execute a function on workers, taking `val` as a parameter. Results are not
collected. This is optimal for various side-effect-causing computations that
are not easily expressible with `dtransform`.
"""
function dexec(val, fn, workers)
    for f in [get_from(pid, :(
        begin
            $fn($val)
            nothing
        end
    )) for pid in workers]
        fetch(f)
    end
end

"""
    dexec(dInfo::Dinfo, fn)

Variant of `dexec` that works with `Dinfo`.
"""
function dexec(dInfo::Dinfo, fn)
    dexec(dInfo.val, fn, dInfo.workers)
end

"""
    dtransform(val, fn, workers, tgt::Symbol=val)

Transform the worker-local distributed data available as `val` on `workers`
in-place, by a function `fn`. Store the result as `tgt` (default `val`)

# Example

    # multiply all saved data by 2
    dtransform(:myData, (d)->(2*d), workers())
"""
function dtransform(val, fn, workers, tgt::Symbol = val)::Dinfo
    for f in [save_at(pid, tgt, :($fn($val))) for pid in workers]
        fetch(f)
    end
    return Dinfo(tgt, workers)
end

"""
    dtransform(dInfo::Dinfo, fn, tgt::Symbol=dInfo.val)::Dinfo

Same as `dtransform`, but specialized for `Dinfo`.
"""
function dtransform(dInfo::Dinfo, fn, tgt::Symbol = dInfo.val)::Dinfo
    dtransform(dInfo.val, fn, dInfo.workers, tgt)
end

"""
    dmapreduce(val, map, fold, workers)

A distributed work-alike of the standard `mapreduce`: Take a function `map` (a
non-modifying transform on the data) and `fold` (2-to-1 reduction of the
results of `map`), systematically run them on the data described by `val`
distributed on `workers`, and return the final reduced result.

It is assumed that the fold operation is associative, but not commutative (as
in semigroups). If there are no workers, operation returns `nothing` (we don't
have a monoid to magically conjure zero elements :[ ).

In current version, the reduce step is a sequential left fold, executed in the
main process.

# Example
    # compute the mean of all distributed data
    sum,len = dmapreduce(:myData,
        (d) -> (sum(d),length(d)),
        ((s1, l1), (s2, l2)) -> (s1+s2, l1+l2),
        workers())
    println(sum/len)

# Processing multiple arguments (a.k.a. "zipWith")

The `val` here does not necessarily need to refer to a symbol, you can easily
pass in a quoted tuple, which will be unquoted in the function parameter. For
example, distributed values `:a` and `:b` can be joined as such:

    dmapreduce(:((a,b)),
        ((a,b)::Tuple) -> [a b],
        vcat,
        workers())
"""
function dmapreduce(val, map, fold, workers)
    if isempty(workers)
        return nothing
    end

    futures = [get_from(pid, :($map($val))) for pid in workers]
    res = fetch(futures[1])

    # replace the collected futures with new empty futures to allow them to be
    # GC'd and free memory for more incoming results
    futures[1] = Future()

    for i = 2:length(futures)
        res = fold(res, fetch(futures[i]))
        futures[i] = Future()
    end
    res
end

"""
    dmapreduce(dInfo::Dinfo, map, fold)

Distributed map/reduce (just as the other overload of `dmapreduce`)
that works with `Dinfo`.
"""
function dmapreduce(dInfo::Dinfo, map, fold)
    dmapreduce(dInfo.val, map, fold, dInfo.workers)
end

"""
    dmapreduce(vals::Vector, map, fold, workers)

Variant of `dmapreduce` that works with more distributed variables
at once.
"""
function dmapreduce(vals::Vector, map, fold, workers)
    return dmapreduce(Expr(:vect, vals...), vals -> map(vals...), fold, workers)
end

"""
    dmapreduce(dInfo1::Dinfo, dInfo2::Dinfo, map, fold)

Variant of `dmapreduce` that works with more `Dinfo`s at
once.  The data must be distributed on the same set of workers, in the same
order.
"""
function dmapreduce(dInfos::Vector{Dinfo}, map, fold)
    if (isempty(dInfos))
        return nothing
    end

    if any([dInfos[1].workers] .!= [di.workers for di in dInfos])
        @error "workers in Dinfo objects do not match" dInfos[1].workers
        error("data distribution mismatch")
    end

    return dmapreduce([di.val for di in dInfos], map, fold, dInfos[1].workers)
end

"""
    gather_array(val::Symbol, workers, dim=1; free=false)

Collect the arrays distributed on `workers` under value `val` into an array. The
individual arrays are pasted in the dimension specified by `dim`, i.e. `dim=1`
is roughly equivalent to using `vcat`, and `dim=2` to `hcat`.

`val` must be an Array-based type; the function will otherwise fail.

If `free` is true, the `val` is `unscatter`ed after being gathered.

This preallocates the array for results, and is thus more efficient than e.g.
using `dmapreduce` with `vcat` for folding.
"""
function gather_array(val::Symbol, workers, dim = 1; free = false)
    size0 = get_val_from(workers[1], :(size($val)))
    innerType = get_val_from(workers[1], :(typeof($val).parameters[1]))
    sizes = dmapreduce(val, d -> size(d, dim), vcat, workers)
    ressize = [size0[i] for i = 1:length(size0)]
    ressize[dim] = sum(sizes)
    result = zeros(innerType, ressize...)
    off = 0
    for (i, pid) in enumerate(workers)
        idx = [(1:ressize[j]) for j = 1:length(ressize)]
        idx[dim] = ((off+1):(off+sizes[i]))
        result[idx...] = get_val_from(pid, val)
        off += sizes[i]
    end
    if free
        unscatter(val, workers)
    end
    return result
end

"""
    gather_array(dInfo::Dinfo, dim=1; free=false)

Distributed gather_array (just as the other overload) that works with
`Dinfo`.
"""
function gather_array(dInfo::Dinfo, dim = 1; free = false)
    return gather_array(dInfo.val, dInfo.workers, dim, free = free)
end

"""
    dmap(arr::Vector, fn, workers)

Call a function `fn` on `workers`, with a single parameter arriving from the
corresponding position in `arr`.
"""
function dmap(arr::Vector, fn, workers)
    map(fetch, [get_from(w, :($fn($(arr[i])))) for (i, w) in enumerate(workers)])
end

"""
    dpmap(fn, args...; mod = Main, kwargs...)

"Distributed pool map."

A wrapper for `pmap` from `Distributed` package that executes the code in the
correct module, so that it can access the distributed variables at remote
workers. All arguments other than the first function `fn` are passed to `pmap`.

The function `fn` should return an expression that is going to get evaluated.

# Example

```julia
using Distributed
dpmap(x -> :(computeSomething(someData, \$x)), WorkerPool(workers), Vector(1:10))
```

```julia
di = distributeSomeData()
dpmap(x -> :(computeSomething(\$(di.val), \$x)), CachingPool(di.workers), Vector(1:10))
```
"""
function dpmap(fn, args...; mod = Main, kwargs...)
    return pmap(x -> Base.eval(mod, fn(x)), args...; kwargs...)
end

"""
    tmp_symbol(s::Symbol; prefix="", suffix="_tmp")

Decorate a symbol `s` with prefix and suffix, to create a good name for a
related temporary value.
"""
function tmp_symbol(s::Symbol; prefix = "", suffix = "_tmp")
    return Symbol(prefix * String(s) * suffix)
end

"""
    tmp_symbol(dInfo::Dinfo; prefix="", suffix="_tmp")

Decorate the symbol from `dInfo` with prefix and suffix.
"""
function tmp_symbol(dInfo::Dinfo; prefix = "", suffix = "_tmp")
    return tmp_symbol(dInfo.val, prefix = prefix, suffix = suffix)
end

"""
    @remote module expr

A version of [`@remote`](@ref) that adds additional choice of the module for
scope.
"""
macro remote(mod, x)
    :(Base.eval($mod, $(QuoteNode(x))))
end

"""
    @remote expr

In a function that will get evaluated on a remote worker, this ensures the
evaluation scope of the expression `expr` (usually a variable) is taken on
the remote side, preventing namespace clash with the local session.

This is mainly useful for making the functions from `Distributed` package (such
as `pmap` and `remotecall`) work with the data stored by `DistributedData`
package.

Internally, this is handled by wrapping in `eval`.

# Example
```
julia> save_at(2, :x, 321)
Future(2, 1, 162, nothing)

julia> let x=123
         remotecall_fetch(() -> x + (@remote x), 2)
       end
444
```
"""
macro remote(x)
    :(@remote Main $x)
end


text:
type Interpreter <: AbstractPlanOutMapper
    serialization::Dict
    inputs::Dict
    env::Assignment
    salt::AbstractString
    in_experiment::Bool
    evaluated::Bool
end

default_salt = "global_salt"
Interpreter(serialization::Dict) = Interpreter(serialization, Dict(), Assignment(default_salt), default_salt, false, false)
Interpreter(serialization::Dict, inputs::Dict) = Interpreter(serialization, inputs, Assignment(default_salt), default_salt, false, false)
Interpreter(serialization::Dict, salt::AbstractString) = Interpreter(serialization, Dict(), Assignment(salt), salt, false, false)
Interpreter(serialization::Dict, salt::AbstractString, inputs::Dict) = Interpreter(serialization, inputs, Assignment(salt), salt, false, false)

function str_to_op(s::AbstractString)
  if s == AbstractString("<")
    return AbstractString("lessthan")
  elseif s == AbstractString(">")
    return AbstractString("greaterthan")
  elseif s == AbstractString(">=")
    return AbstractString("greaterthanorequalto")
  elseif s == AbstractString("<=")
    return AbstractString("lessthanorequalto")
  elseif s == AbstractString("%")
    return AbstractString("mod")
  elseif s == AbstractString("/")
    return AbstractString("div")
  elseif s == AbstractString("length")
    return AbstractString("len")
  elseif s == AbstractString("return")
    return AbstractString("ret")
  end
  return s
end

setindex!(i::Interpreter, value, key) = setindex!(i.env, value, key)
getindex(i::Interpreter, key) = getindex!(i.env, key)
delete!(i::Interpreter, key) = delete!(i.env, key)

function get_params(i::Interpreter)
  if !i.evaluated
    try
      evaluate(i, i.serialization)
    catch e
      if isa(e, StopPlanOutException)
        i.in_experiment = e.in_experiment
      end
    end
    i.evaluated = true
  end
  return i.env
end

function evaluate(i::Interpreter, planout_code::Any)
  return planout_code
end

function evaluate(i::Interpreter, planout_code::AbstractArray)
  return [evaluate(i, j) for j in planout_code]
end

function evaluate(i::Interpreter, planout_code::Associative)
  if haskey(planout_code, "op")
    op = str_to_op(getindex(planout_code, "op"))
    fxn_call = "PlanOut.$op($planout_code)"
    return execute(eval(parse(fxn_call)), i)
  else
    return planout_code
  end
end


text:
using Documenter, EDFPlus

makedocs(sitename = "EDFPlus Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
# FlyToken type and associated methods

"maximal count of code units directly encoded in a FlyToken"
const MAX_DIRECT_SIZE = 7%UInt64

"maximal count of code units in any buffer based Token"
const MAX_TOKEN_SIZE = ((1<<27) - 1) %UInt64

"bitmask to check if any inline code unit is not ascii"
const NOTASCII_BITS = UInt64(0b10000000100000001000000010000000100000001000000010000000)

"bitmask to restrict to all code units in a DirectFly"
const CODEUNIT_BITS =((1<<56)-1)%UInt64

"bitmask to isolate size bits in DirectFly"
const DIRECT_SIZE_BITS = MAX_DIRECT_SIZE<<56

"bitmask to isolate size bits in BufferFly"
const BUFFER_SIZE_BITS = MAX_TOKEN_SIZE<<32

"bitmask to isolate splitted size bits in BufferFly"
const SPLIT_SIZE_BITS = (((1<<24)-1)<<32)%UInt64

"bitmask to restrict to all code units in a DirectFly"
const OFFSET_BITS =((1<<32)-1)%UInt64

"bitmask to isolate category and fly type bit in any FlyToken"
const CATEGORY_FLY_BITS = (31<<59)%UInt64

"bitmask to isolate category in any FlyToken"
const CATEGORY_BITS = (15%UInt64)<<59

"bitmask to test for tiny. Bit set -> (offset,size) pair stored"
const NOTDIRECT_BIT = (1%UInt64)<<63

"true: size bitfield is split in [`BufferFly`](@ref)"
const BUFFER_SPLIT_SIZE = false

"Known standard types usable in reinterpret for 64 bit types"
const Bits64 = Union{UInt64, Int64}

"""
A flyweight data structure for tokens.

All subtypes must be 64 bit primitive data types which implement a certain
bitmap layout and methods for processing it.

FlyToken implementations will throw an error if a Token API method requires
a buffer and no buffer is supplied to the method. Not all AbstractToken API
methods are implemented by all FlyToken subtypes, many FlyToken subtypes
can only be used as  part of a compound token, not "standalone".

"""
abstract type FlyToken <: AbstractToken end

"""
Flyweight token with direct ("inline") encoding of its code units.

This string data type directly stores very short strings with a length
of up to 7 code units in its value. It supports the AbstractString API and can be
used as a memory efficient substitute for String instances, avoiding
any pointer overhead. It implements the full AbstractToken API.

# performance considerations

DirectFly instances consume much less memory, compared to String
instances. They are true values (no heap allocation), which can improve
data locality and thus CPU cache use. Comparison is considerably faster than
String comparison.

Depending on the API used, there is a disadvantage because it is not possible
to obtain a C-stype pointer to the contents of a DirectFly - without
creating a temporary String copy. E. G. hashing suffers from that.

# implementation details

Design goal is to treat a DirectFly instance as a 64 bit primitive value.
To allow bit operations, memory layout is in host endianness like Int64.
Take that into account on binary serialization.
The sign bit (most significant bit) is used to distinguish a DirectFly
value from other FlyToken types. For a DirectFly, it must always be 0.

DirectFly encodes really tiny strings of up to 7 code units, directly
in its 64 bit value. If size is less than 7, unused code units must be 0.

The memory layout is similar to the julia memory layout for the Char type,
with the difference that the code unit count is explicitly stored in a separate
byte, together with category and the "isdirect flag" in the sign bit.

Order is chosen in a way that a lexikographically correct string compare
for DirectFly-s can be done with one single UInt64 compare operation.

Memory layout:

```
bit 63 63 62 60 59 58 57 56 byte6 byte5 byte4 byte3 byte2 byte1 byte0
    == =========== ========
    |  |           |        |     |     |     |     |     |     |
    |  |           |        cu[1] cu[2] cu[3] cu[4] cu[5] cu[6] cu[7]
    |  |           size: 0..7
    |  category: 0..15
    |
    =0: marker flag must be 0 in DirectFly instances
```

"""
primitive type DirectFly <: FlyToken 64 end

"alias name similar to HToken and BToken"
const DToken = DirectFly

"""
Flyweight token accomplished by a code unit buffer.

This token type is a flyweight data structure which needs an accompanying
code unit buffer. Any method which operates on the contents, needs access to
the associated buffer. Usually, the buffer is supplied as an additional
parameter.

BufferFly is considered a module-private type and thus not exported by Tokens
module. It does NOT implement the full AbstractString API! Only
AbstractBuffer and AbstractToken methods which do not need contents access
are implemented.

See [`BToken`](@ref) for an implementation with explicitly associated buffer.

Memory layout:

```
bit 63 63 62 60 59 58 57 56 byte6 byte5 byte4 byte3 byte2 byte1 byte0
    =1 =========== ========================== =======================
    |  |           |                          |
    |  |           |                          offset: 32 bits
    |  |           size: 26 bits
    |  category: 0..15
    |
    =1: must be 1 for any BToken instance if size>0 or offset>0
    code units are stored in buffer[offset+1] .. buffer[offset+size]
    buffer must be known in the processing context.
```

if [`BUFFER_SPLIT_SIZE`](@ref) is true, the following bit layout is used:

    ```
    bit 63 63 62 60 59 58 57 56 byte6 byte5 byte4 byte3 byte2 byte1 byte0
        =1 =========== ======== ================= =======================
        |  |           |        |                 |
        |  |           |        |                 offset: 0..2^32-1
        |  |           |        size>>3: bits 3..26 of size
        |  |           size&7: lowest 3 bits
        |  category: 0..15
        |
        =1: must be 1 for any BToken instance if size>0 or offset>0
        code units are stored in buffer[offset+1] .. buffer[offset+size]
        buffer must be known in the processing context.
    ```

"""
primitive type BufferFly <: FlyToken 64 end

"""
Flyweight/direct token union type.

This token type is an alternative implementation for the julia
construct *union {DirectFly, BufferFly}*.

The sign bit in the 64-bit value (interpreted as UInt64) distinguishes
between the two types that make up the union, instead of a separate
tag byte which would be used by a Julia union construct.

The type acts as a DirectFly or a BufferFly depending on its value,
which is tested at runtime, introducing the overhead of one test and
conditional jump on access.

HybridFly is considered a private type of module Tokens and thus not exported.
It does NOT implement the full AbstractString API! Only AbstractBuffer and
AbstractToken methods which do not need contents access are implemented.

Use HybridFly instead of BufferFly in cases where strings with up to
7 code units have a significant proportion. Benchmark both alternatives
if runtime performance is critical - it depends on your operation mix
whether HybridFly gives faster code or not.

[`Token`](@ref) and [`TokenVector`](@ref) are using HybridFly.

# implementation details

Because all values of DirectFly and BufferFly, seen as native 64-bit-chunk,
are not overlapping, the concrete type of a HybridFly is easily determined
at runtime (test the sign bit in Int64 interpretation).
Use [`isdirect`](@ref)f or this purpose.

usize can be implemented without conditional jumps, by tricky bit
operations. This is left to llvm optimization - benchmarks have shown that
llvm is capable of doing it and decides on code generation which implementation
is regarded faster.

category bits are identical for DirectFly and BufferFly, so no conditional
code is needed.

If offset and size bitfields are both 0, DirectFly and BufferFly encode
the same content: an empty string.

TODO unique nothing/missing encoding with nondirect flag cleared even for BufferFly?!

TODO we have unnormalized values which may NEVER occur with proper constructor calls:
length==0, but offset bytes !=0. Currently, not used. Could be used to store
an Int32 or a Float in binary format.
"""
primitive type HybridFly <: FlyToken 64 end

## basic helpers ##

"""
Private convenience converter to an UInt64 value.

All FlyToken operations are finally implemented with native 64 bit
operations defined for Int64 or UInt64. We need a short notation
to convert between (U)Int64 and FlyToken. No checks!!
"""
u64(t::FlyToken) =  reinterpret(UInt64,t)

"""
Private unsafe convenience converter to a DirectFly.

Argument must be a 64 bit primitive value.
No checks performed!
"""
df(bits) = reinterpret(DirectFly,bits)

"""
Private unsafe convenience converter to a HybridFly.

Argument must be a 64 bit primitive value.
No checks performed!
"""
hf(bits) = reinterpret(HybridFly,bits)

"""
Private unsafe convenience converter to a BufferFly.

Argument must be a 64 bit primitive value.
No checks performed!
"""
bf(bits) = reinterpret(BufferFly,bits)

"""set a codeunit within a DirectFly assuming current value is 0.

index must be 1..7, t mush have codeunit value 0 at index i.
titled unsafe, because no checks are performed.
"""
function unsafe_setcodeunit(t::DirectFly, index, codeunit::UInt8)
    t | ( UInt64(codeunit)<< ((7-index)<<3))
end

offset(t::BufferFly) = (u64(t) & OFFSET_BITS)%UInt32

offset(t::DirectFly) = UInt32(0)

offset(t::HybridFly) = isdirect(t) ? offset(df(t)) : offset(bf(t))

## constructors  ##

"empty token (offset, length are 0)"
function BufferFly(cat::Nibble)
    bf(NOTDIRECT_BIT | (cat%UInt64)<<59)
end

"""
token with given size, bounds-checked, offset 0
"""
Base.@propagate_inbounds function BufferFly(cat::Nibble, size::UInt64)
    @boundscheck checkulimit(size,MAX_TOKEN_SIZE)
    if BUFFER_SPLIT_SIZE
        BufferFly(cat) | ((size&7)<<56 | (size>>>3)<<32)
    else
        BufferFly(cat) | (size<<32)
    end
end

"""
raw incomplete token with category and size in packed format, offset 0
"""
function BufferFly(cat_size::Packed31)
    u = NOTDIRECT_BIT | (cat_size%UInt64) <<59 # tricky: all bits above category are shifted out or set by NOTTINYBIT
    s = bits4_30(cat_size)%UInt64
    u |= BUFFER_SPLIT_SIZE ? ((s&7)<<56) | ((s>>>3)<<32)  :  (s<<32)
    bf(u)
end

"raw DirectFly with all content bytes 0"
Base.@propagate_inbounds function DirectFly(cat_size::Packed31)
    u = ((cat_size%UInt64) <<59) & CATEGORY_BITS
    s = bits4_30(cat_size)%UInt64
    @boundscheck checkulimit(s,MAX_DIRECT_SIZE)
    u |= ((s&7)<<56)
    df(u)
end

"raw Hybridfly: direct if size <8"
function HybridFly(cat_size::Packed31)
    @inbounds hf( UInt32(cat_size)<= 0x0000007F ? DirectFly(cat_size) : BufferFly(cat_size))
end

"""
Lowlevel constructor for a token referencing some external buffer.

UNSAFE because the referenced buffer is not given.
Method cannot check whether (offset,size) is valid for the buffer which is
used with this token. You have to grant (or check) validity of (offset,size)
with respect to the referenced buffer in your code elsewere.

"""
Base.@propagate_inbounds function BufferFly(cat::Nibble, offset::UInt32, size::UInt64)
    BufferFly(cat,size)|offset
end

"""
Constructor of an empty token with given category.

All codeunit bytes are 0, length is 0.
"""
DirectFly(cat::Nibble) = df( UInt64(cat)<<59 )

"token with given size, bounds-checked, all code units are 0, prepared for unsafe_setcodeunit"
Base.@propagate_inbounds function DirectFly(cat::Nibble, size::UInt64)
    @boundscheck checkulimit(size,MAX_DIRECT_SIZE)
    df(DirectFly(cat)) | (size << 56)
end

"Construct a DirectFly for size below 8, else a BufferFly with offset 0"
Base.@propagate_inbounds function HybridFly(cat::Nibble, size::UInt64)
    size <= MAX_DIRECT_SIZE ?
        hf(DirectFly(cat,size)) :
        hf(BufferFly(cat,size))
end

"""
constructor for direct encoding from Utf8-encoded strings.

bounds checks performed
"""
Base.@propagate_inbounds function DirectFly(cat::Nibble, offset::UInt32, size::UInt64, s::Utf8String)
    size == 0 && return DirectFly(cat)
    @boundscheck checkbounds(s,offset+size)
    fly = DirectFly(cat,size)
    for i in 1:size
        fly = unsafe_setcodeunit(fly,i,codeunit(s,offset+i))
    end
    fly
end

"token from a string constant (requires size<8)"
function DirectFly(cat::Nibble, s::Utf8String)
    size = usize(s)
    @boundscheck checkulimit(size,MAX_DIRECT_SIZE)
    fly = DirectFly(cat,size)
    for i in 1:size
        fly = unsafe_setcodeunit(fly,i,codeunit(s,i))
    end
    fly
end

function DirectFly(cat::Nibble, c::Char)
    df(u64(DirectFly(cat,Base.codelen(c))) | UInt64(reinterpret(UInt32, c)) << 24)
end

## special values for missing and nothing

const DIRECT_NOTHING = DirectFly(Nibble(T_END))
const DIRECT_MISSING = DirectFly(Nibble(T_SPECIAL))
Base.isnothing(t::FlyToken) = t & (DIRECT_SIZE_BITS | CATEGORY_BITS) == DIRECT_NOTHING
Base.ismissing(t::FlyToken) = t & (DIRECT_SIZE_BITS | CATEGORY_BITS) == DIRECT_MISSING
DirectFly(::Nothing) = DIRECT_NOTHING
DirectFly(::Missing) = DIRECT_MISSING
HybridFly(::Nothing) = hf(DIRECT_NOTHING)
HybridFly(::Missing) = hf(DIRECT_MISSING)
BufferFly(::Nothing) = bf(DIRECT_NOTHING|NOTDIRECT_BIT)
BufferFly(::Missing) = bf(DIRECT_MISSING|NOTDIRECT_BIT)

Base.@propagate_inbounds DirectFly(offset::UInt32, size::UInt64,t::AbstractToken) = DirectFly(t.cat,offset,size,t)

## Token API methods ##

category(t::FlyToken) = reinterpret(Nibble, UInt8((u64(t) >>>59)&15))

category(t::DirectFly) = reinterpret(Nibble, UInt8(u64(t) >>>59))

isdirect(t::HybridFly) = reinterpret(Int64,t)>=0

isdirect(t::DirectFly) = true

isdirect(t::BufferFly) = false

usize(t::DirectFly) = (u64(t)>>56) & MAX_DIRECT_SIZE

function usize(t::BufferFly)
    if BUFFER_SPLIT_SIZE
        usize(df(t)) | (u64(t)&SPLIT_SIZE_BITS)>>29
    else
        (u64(t)&BUFFER_SIZE_BITS)>>32
    end
end

function usize(t::HybridFly)
    if BUFFER_SPLIT_SIZE
        # tricky code without jumps:
        # 3 lowest bits of the size are always stored at bit position 56..58.
        # if NONASCII_BIT is set, we have additional 24 bits for size at bits 32..55
        # We build a mask for bits 32..55 all 1 (bit 63 set, data in buffer)
        # or 0 (bit 63 clear, all data in token, size is 0..7)
        # by arithmetic shift of bit 63. ANDing with t.bits and shift by 29 gives
        # the high bits 3..26 of the size.
        masksize = (reinterpret(Int64,t)>>31) & ((1<<24-1)<<32)
        ((u64(t)&masksize)>>>29) | ((u64(t) >>>56) & 7)
    else
        if isdirect(t)
            usize(df(t))
        else
            usize(bf(t))
        end
    end
end

## helpers for IO
Packed31(f::FlyToken)=Packed31(category(f)%UInt8, f.len%UInt32)

## Base operators and functions overloading ##

function Base.String(d::DirectFly)
    size = d.len
    s = Base._string_n(size)
    p = pointer(s)
    GC.@preserve s begin
        for i in 1::size
            @inbounds unsafe_store!(p,codeunit(d,i))
            p += 1
        end
    end
    return s
end

"bitwise OR applied to a FlyToken. DANGEROUS: no validity checks!!"
Base.:|(f::F, orValue::T) where {F<:FlyToken,T<:Unsigned} = reinterpret(F,u64(f) | orValue)

"bitwise AND applied to a FlyToken. DANGEROUS: no validity checks!!"
Base.:&(f::F, andValue::T) where {F<:FlyToken,T<:Unsigned} = reinterpret(F,u64(f) & andValue)

"add an offset to a FlyToken (error if isdirect(f)). DANGEROUS: no validity or overflow checks!!"
function Base.:+(f::F, addValue::UInt32) where {F<:FlyToken}
    isdirect(f) && error("cannot add offset to DirectFly")
    reinterpret(F,u64(f) + addValue)
end

## interpretation

Base.@propagate_inbounds function Base.iterate(s::AbstractToken, i::Int=firstindex(s))
    i > ncodeunits(s) && return nothing
    b = codeunit(s, i)
    u = UInt32(b) << 24
    # code from String
    #Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    # does work for valid Utf8 but does not detect illegal code: bytes C0, C1, F5..F8 are all illegal,
    # according to wikipedia (referencing RFC3629), but are accepted here.
    # we can spare one comparison by this code:
    b >= 0x80 || return reinterpret(Char, u), i+1
    return Base.iterate_continued(s, i, u)
end

function Base.iterate_continued(s::AbstractToken, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret) # u<0x80000000 treated in iterate. all other u<0xc2000000 are illegal!!
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

## Serialization of DirectFly

function Base.write(io::IO, t::DirectFly)
    write(io,Packed31(t))
    for i in 1:t.len
        write(io, @inbounds codeunit(t,i))
    end
end

"read contents of a DirectFly for already given category and size"
function Base.read(io::IO, cat_size::Packed31,::Type{DirectFly})
    size = bits4_30(cat_size)
    t = DirectFly(cat_size)
    for i in 1:size
        t = unsafe_setcodeunit(t,i,read(io,UInt8))
    end
    t
end

"read complete DirectFly"
function Base.read(io::IO, ::Type{DirectFly})
    read(io, read(io,Packed31),DirectFly) # category and size
end

Base.@propagate_inbounds function Base.codeunit(t::DirectFly, index::Int)
    @boundscheck checkbounds(t,index)
    (255%UInt8) & (u64(t) >>> ((7-index)<<3))%UInt8
end

Base.@propagate_inbounds function byte(t::DirectFly, ofs::UInt32)
    @boundscheck checkbyteofs(ofs,t.len)
    (255%UInt8) & (u64(t) >>> ((6-index)<<3))%UInt8
end

function Base.codeunit(t::FlyToken, index::Int)
    isdirect(t) || error("code unit buffer not available in codeunit($t,$index)")
    codeunit(df(t),index)
end

function Base.show(io::IO,t::T) where T <: FlyToken
    print(io, string(typeof(t))[1] * Char(t.cat))
    if isdirect(t)
        Base.print_quoted(io, df(t))
    else
        print(io,'<',offset(t)%Int,',',t.len%Int,'>')
    end
end

"generic comparison by code units - used e.g. for DirectFly"
function cmp_codeunits(a::Utf8String, b::Utf8String)
    al, bl = sizeof(a), sizeof(b)
    ml = min(al, bl)
    i = 1
    @inbounds while i<=ml
        ai = codeunit(a,i)
        bi = codeunit(b,i)
        ai < bi && return -1
        ai > bi && return 1
        i += 1
    end
    al < bl && return -1
    al > bl && return 1
    0
end

function Base.cmp(a::DirectFly, b::DirectFly)
    # both tiny: compare all code units in one step
    (u64(a)&CODEUNIT_BITS) < (u64(b)&CODEUNIT_BITS) && return -1
    (u64(a)&CODEUNIT_BITS) > (u64(b)&CODEUNIT_BITS) && return 1
    0
end

function Base.:(==)(a::DirectFly, b::DirectFly)
    # both tiny: compare all code units in one step
    (u64(a)&CODEUNIT_BITS) == (u64(b)&CODEUNIT_BITS)
end

# other fly tokens have no hash
function Base.hash(t::DirectFly, h::UInt)
    h += Base.memhash_seed
    # note: use pointer(s) here (see #6058).
    mmhash(t, h % UInt32) + h
end

# MurmurHash3 up to 8 bytes
#
function mmhash(str::DirectFly, seed::UInt32)
    h1 = h2 = seed%UInt64
    len = usize(str)
    k1 = 0%UInt64
    if len != 0
        k1 = Base.bswap(u64(str)<<8)
    end
    h1 = mhtail1(h1, k1)
    mhfin(len, h1, h2)
end

#########################################################
################ check methods ##########################
#########################################################

#=

endianness
Base.lastindex(s::FlyToken) = Int(s.size_content & 0xf)
Base.iterate(s::FlyToken, i::Integer) = iterate(String(s), i)
Base.iterate(s::FlyToken) = iterate(String(s))
Base.print(s::FlyToken) = print(String(s))
Base.display(s::FlyToken) = display(String(s))
Base.firstindex(::FlyToken) = 1
Base.collect(s::FlyToken) = getindex.(s, 1:lastindex(s))

==(s::FlyToken, b::String) = begin
    String(s)  == b
end
=#


text:
module modelling

function step(x)
    return x >= 0 ? 1.0 : 0.0
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#*****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using TFTools
using Datasets
using Base.Test
using TensorFlow
import TensorFlow: DT_FLOAT32
import TensorFlow.API: l2_loss, AdamOptimizer, cast, round_

# Parameters
const LEARNING_RATE = 0.001
const TRAINING_EPOCHS = 200
const BATCH_SIZE = 2000
const DISPLAY_STEP = 1

# Network parameters
const HIDDEN_UNITS = [100, 100]

function test_softmux(datname::AbstractString="bin_synth",
    featfile::AbstractString="feats", labelfile::AbstractString="labels";
    labelfield::AbstractString="x1",
    b_debug::Bool=false)

    Dfeats = dataset(datname, featfile)
    Dlabels = dataset(datname, labelfile)
    n_feats = ncol(Dfeats)
    @assert nrow(Dfeats) == nrow(Dlabels)

    data_set = TFDataset(Dfeats, Dlabels[symbol(labelfield)])

    # Construct model
    muxselect = Placeholder(DT_FLOAT32, [-1, n_feats])
    muxin = Placeholder(DT_FLOAT32, [-1, n_feats])
    mux = Softmux(n_feats, n_feats, HIDDEN_UNITS, Tensor(muxin), Tensor(muxselect))
    pred = out(mux)
    labels = Placeholder(DT_FLOAT32, [-1])

    # Define loss and optimizer
    cost = l2_loss(pred - labels) # Squared loss
    optimizer = minimize(AdamOptimizer(LEARNING_RATE), cost) # Adam Optimizer

    # Initializing the variables
    init = initialize_all_variables()

    sess = Session()
    try
        run(sess, init)

        # Training cycle
        for epoch in 1:TRAINING_EPOCHS
            avg_cost = 0.0
            total_batch = div(num_examples(data_set), BATCH_SIZE)

            # Loop over all batches
            for i in 1:total_batch
                batch_xs, batch_ys = next_batch(data_set, BATCH_SIZE)
                # Fit training using batch data
                run(sess, optimizer, FeedDict(muxin => batch_xs, muxselect => batch_xs, labels => batch_ys))
                # Compute average loss
                batch_average_cost = run(sess, cost, FeedDict(muxin => batch_xs, muxselect => batch_xs, labels => batch_ys))
                avg_cost += batch_average_cost / (total_batch * BATCH_SIZE)
            end

            # Display logs per epoch step
            if epoch % DISPLAY_STEP == 0
                println("Epoch $(epoch)  cost=$(avg_cost)")
            end
        end
        println("Optimization Finished")

        # Test model
        correct_prediction = (round_(pred) == labels)
        # Calculate accuracy
        accuracy = mean(cast(correct_prediction, DT_FLOAT32))
        acc = run(sess, accuracy, FeedDict(muxin => data_set.X, muxselect => data_set.X, labels => data_set.Y))
        println("Accuracy:", acc)
        if b_debug
            db_x = data_set.X
            db_nnout = run(sess, mux.nnout, FeedDict(muxin => data_set.X, muxselect => data_set.X, labels => data_set.Y))
            db_ypred = run(sess, mux.muxout, FeedDict(muxin => data_set.X, muxselect => data_set.X, labels => data_set.Y))
            db_labels = data_set.Y

            NSHOW = 20
            @show db_x[1:NSHOW]
            @show db_nnout[1:NSHOW,:]
            @show db_ypred[1:NSHOW]
            @show db_labels[1:NSHOW]
        end
    finally
        close(sess)
    end

end #module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FileStatusProperties <: SwaggerModel
    accessTime::Any # spec type: Union{ Nothing, Int64 } # spec name: accessTime
    blockSize::Any # spec type: Union{ Nothing, Int64 } # spec name: blockSize
    msExpirationTime::Any # spec type: Union{ Nothing, Int64 } # spec name: msExpirationTime
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    length::Any # spec type: Union{ Nothing, Int64 } # spec name: length
    modificationTime::Any # spec type: Union{ Nothing, Int64 } # spec name: modificationTime
    owner::Any # spec type: Union{ Nothing, String } # spec name: owner
    pathSuffix::Any # spec type: Union{ Nothing, String } # spec name: pathSuffix
    permission::Any # spec type: Union{ Nothing, String } # spec name: permission
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    aclBit::Any # spec type: Union{ Nothing, Bool } # spec name: aclBit

    function FileStatusProperties(;accessTime=nothing, blockSize=nothing, msExpirationTime=nothing, group=nothing, length=nothing, modificationTime=nothing, owner=nothing, pathSuffix=nothing, permission=nothing, type=nothing, aclBit=nothing)
        o = new()
        validate_property(FileStatusProperties, Symbol("accessTime"), accessTime)
        setfield!(o, Symbol("accessTime"), accessTime)
        validate_property(FileStatusProperties, Symbol("blockSize"), blockSize)
        setfield!(o, Symbol("blockSize"), blockSize)
        validate_property(FileStatusProperties, Symbol("msExpirationTime"), msExpirationTime)
        setfield!(o, Symbol("msExpirationTime"), msExpirationTime)
        validate_property(FileStatusProperties, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(FileStatusProperties, Symbol("length"), length)
        setfield!(o, Symbol("length"), length)
        validate_property(FileStatusProperties, Symbol("modificationTime"), modificationTime)
        setfield!(o, Symbol("modificationTime"), modificationTime)
        validate_property(FileStatusProperties, Symbol("owner"), owner)
        setfield!(o, Symbol("owner"), owner)
        validate_property(FileStatusProperties, Symbol("pathSuffix"), pathSuffix)
        setfield!(o, Symbol("pathSuffix"), pathSuffix)
        validate_property(FileStatusProperties, Symbol("permission"), permission)
        setfield!(o, Symbol("permission"), permission)
        validate_property(FileStatusProperties, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(FileStatusProperties, Symbol("aclBit"), aclBit)
        setfield!(o, Symbol("aclBit"), aclBit)
        o
    end
end # type FileStatusProperties

const _property_map_FileStatusProperties = Dict{Symbol,Symbol}(Symbol("accessTime")=>Symbol("accessTime"), Symbol("blockSize")=>Symbol("blockSize"), Symbol("msExpirationTime")=>Symbol("msExpirationTime"), Symbol("group")=>Symbol("group"), Symbol("length")=>Symbol("length"), Symbol("modificationTime")=>Symbol("modificationTime"), Symbol("owner")=>Symbol("owner"), Symbol("pathSuffix")=>Symbol("pathSuffix"), Symbol("permission")=>Symbol("permission"), Symbol("type")=>Symbol("type"), Symbol("aclBit")=>Symbol("aclBit"))
const _property_types_FileStatusProperties = Dict{Symbol,String}(Symbol("accessTime")=>"Int64", Symbol("blockSize")=>"Int64", Symbol("msExpirationTime")=>"Int64", Symbol("group")=>"String", Symbol("length")=>"Int64", Symbol("modificationTime")=>"Int64", Symbol("owner")=>"String", Symbol("pathSuffix")=>"String", Symbol("permission")=>"String", Symbol("type")=>"String", Symbol("aclBit")=>"Bool")
Base.propertynames(::Type{ FileStatusProperties }) = collect(keys(_property_map_FileStatusProperties))
Swagger.property_type(::Type{ FileStatusProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_FileStatusProperties[name]))}
Swagger.field_name(::Type{ FileStatusProperties }, property_name::Symbol) =  _property_map_FileStatusProperties[property_name]

const _allowed_FileStatusProperties_type = ["FILE", "DIRECTORY"]

function check_required(o::FileStatusProperties)
    true
end

function validate_property(::Type{ FileStatusProperties }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "FileStatusProperties", :enum, val, _allowed_FileStatusProperties_type)
    end
end


text:
import Base.isempty

export ConvexHull, CH,
       convex_hull,
       convex_hull!,
       ConvexHull!,
       swap

"""
    ConvexHull{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}

Type that represents the convex hull of the union of two convex sets.

### Fields

- `X` -- convex set
- `Y` -- convex set

### Notes

The `EmptySet` is the neutral element for `ConvexHull`.

### Examples

Convex hull of two 100-dimensional Euclidean balls:

```jldoctest
julia> b1, b2 = Ball2(zeros(100), 0.1), Ball2(4*ones(100), 0.2);

julia> c = ConvexHull(b1, b2);

julia> typeof(c)
ConvexHull{Float64,Ball2{Float64,Array{Float64,1}},Ball2{Float64,Array{Float64,1}}}
```
"""
struct ConvexHull{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2

    # default constructor with dimension check
    function ConvexHull(X::S1, Y::S2) where {N<:Real, S1<:LazySet{N},
                                             S2<:LazySet{N}}
        @assert dim(X) == dim(Y) "sets in a convex hull must have the same " *
            "dimension"
        return new{N, S1, S2}(X, Y)
    end
end

isoperationtype(::Type{<:ConvexHull}) = true
isconvextype(::Type{<:ConvexHull}) = true

# EmptySet is the neutral element for ConvexHull
@neutral(ConvexHull, EmptySet)

# Universe is the absorbing element for ConvexHull
@absorbing(ConvexHull, Universe)

"""
    CH

Alias for `ConvexHull`.
"""
const CH = ConvexHull

"""
    swap(ch::ConvexHull)

Return a new `ConvexHull` object with the arguments swapped.

### Input

- `ch` -- convex hull of two convex sets

### Output

A new `ConvexHull` object with the arguments swapped.
"""
function swap(ch::ConvexHull)
    return ConvexHull(ch.Y, ch.X)
end

"""
    dim(ch::ConvexHull)

Return the dimension of a convex hull of two convex sets.

### Input

- `ch` -- convex hull of two convex sets

### Output

The ambient dimension of the convex hull of two convex sets.
"""
function dim(ch::ConvexHull)
    return dim(ch.X)
end

"""
    σ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}

Return the support vector of a convex hull of two convex sets in a given
direction.

### Input

- `d`  -- direction
- `ch` -- convex hull of two convex sets

### Output

The support vector of the convex hull in the given direction.
"""
function σ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}
    σ1 = σ(d, ch.X)
    σ2 = σ(d, ch.Y)
    ρ1 = dot(d, σ1)
    ρ2 = dot(d, σ2)
    return ρ1 >= ρ2 ? σ1 : σ2
end

"""
    ρ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}

Return the support function of a convex hull of two convex sets in a given
direction.

### Input

- `d`  -- direction
- `ch` -- convex hull of two convex sets

### Output

The support function of the convex hull in the given direction.
"""
function ρ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}
    return max(ρ(d, ch.X), ρ(d, ch.Y))
end

"""
    isbounded(ch::ConvexHull)

Determine whether a convex hull of two convex sets is bounded.

### Input

- `ch` -- convex hull of two convex sets

### Output

`true` iff both wrapped sets are bounded.
"""
function isbounded(ch::ConvexHull)
    return isbounded(ch.X) && isbounded(ch.Y)
end

"""
    isempty(ch::ConvexHull)

Return if a convex hull of two convex sets is empty or not.

### Input

- `ch` -- convex hull

### Output

`true` iff both wrapped sets are empty.
"""
function isempty(ch::ConvexHull)
    return isempty(ch.X) && isempty(ch.Y)
end


text:
export
    CuArray, CuVector, CuMatrix, CuVecOrMat,
    device, curand

type CuArray{T,N}
    ptr::CuPtr
    dims::NTuple{N,Int}
end

typealias CuVector{T} CuArray{T,1}
typealias CuMatrix{T} CuArray{T,2}
typealias CuVecOrMat{T} Union{CuVector{T},CuMatrix{T}}

(::Type{CuArray{T}}){T,N}(dims::NTuple{N,Int}) = CuArray{T,N}(alloc(T,prod(dims)), dims)
(::Type{CuArray{T}}){T}(dims::Int...) = CuArray{T}(dims)
@generated function (::Type{CuArray{T}}){T,U}(x::CuArray{U})
    f = CuFunction("""
    __global__ void f($T *y, $U *x, int length) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx < length) {
            y[idx] = x[idx];
        }
    }""")
    quote
        y = CuArray{T}(size(x))
        $f(y.ptr, x.ptr, length(y), dx=length(y))
        y
    end
end

CuArray{T,N}(x::Array{T,N}) = copy!(CuArray{T}(size(x)), x)
Base.Array{T,N}(x::CuArray{T,N}) = copy!(Array{T}(size(x)), x)

device(x::CuArray) = x.ptr.dev
Base.pointer{T}(x::CuArray{T}, index::Int=1) = Ptr{T}(x.ptr) + sizeof(T) * (index-1)

Base.length(x::CuArray) = prod(x.dims)
Base.size(x::CuArray) = x.dims
Base.size(x::CuArray, dim::Int) = x.dims[dim]
Base.ndims{T,N}(x::CuArray{T,N}) = N
Base.eltype{T}(x::CuArray{T}) = T
Base.isempty(x::CuArray) = length(x) == 0

Base.strides(x::CuVector) = (1,)
Base.strides(x::CuMatrix) = (1,size(x,1))
function Base.strides{T}(x::CuArray{T,3})
    s2 = size(x,1)
    s3 = s2 * size(x,2)
    (1,s2,s3)
end
function Base.strides{T}(x::CuArray{T,4})
    s2 = size(x,1)
    s3 = s2 * size(x,2)
    s4 = s3 * size(x,3)
    (1,s2,s3,s4)
end
function Base.stride{T,N}(x::CuArray{T,N}, dim::Int)
    d = 1
    for i = 1:dim-1
        d *= size(x,i)
    end
    d
end

Base.similar{T,N}(x::CuArray{T}, dims::NTuple{N,Int}) = CuArray{T}(dims)
Base.similar(x::CuArray) = similar(x, size(x))
Base.similar(x::CuArray, dims::Int...) = similar(x, dims)

Base.convert{T}(::Type{Ptr{T}}, x::CuArray{T}) = Ptr{T}(x.ptr)
Base.convert(::Type{CUdeviceptr}, x::CuArray) = CUdeviceptr(x.ptr)
Base.convert{T,N}(::Type{CuArray{T,N}}, x::Array{T,N}) = CuArray(x)
Base.convert{T,N}(::Type{Array{T,N}}, x::CuArray{T,N}) = Array(x)
Base.unsafe_convert{T}(::Type{Ptr{T}}, x::CuArray) = Ptr{T}(x.ptr)
Base.unsafe_convert(::Type{CUdeviceptr}, x::CuArray) = CUdeviceptr(x.ptr)

Base.zeros{T,N}(::Type{CuArray{T}}, dims::NTuple{N,Int}) = fill(CuArray, T(0), dims)
Base.zeros{T,N}(x::CuArray{T,N}) = zeros(CuArray{T}, size(x))
Base.zeros{T}(::Type{CuArray{T}}, dims::Int...) = zeros(CuArray{T}, dims)
Base.ones{T}(x::CuArray{T}) = ones(CuArray{T}, x.dims)
Base.ones{T}(::Type{CuArray{T}}, dims::Int...) = ones(CuArray{T}, dims)
Base.ones{T}(::Type{CuArray{T}}, dims) = fill(CuArray, T(1), dims)

function Base.copy!{T}(dest::Array{T}, src::CuArray{T}; stream=C_NULL)
    cuMemcpyDtoHAsync(dest, src, length(src)*sizeof(T), stream)
    dest
end
function Base.copy!{T}(dest::CuArray{T}, src::Array{T}; stream=C_NULL)
    cuMemcpyHtoDAsync(dest, src, length(src)*sizeof(T), stream)
    dest
end
function Base.copy!{T}(dest::CuArray{T}, src::CuArray{T}; stream=C_NULL)
    cuMemcpyDtoDAsync(dest, src, length(src)*sizeof(T), stream)
    dest
end
Base.copy(x::CuArray) = copy!(similar(x),x)

@generated function Base.fill!{T}(x::CuArray{T}, value)
    f = CuFunction("""
    __global__ void f($T *x, int length, $T value) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx < length) {
            x[idx] = value;
        }
    }""")
    quote
        $f(x.ptr, length(x), T(value), dx=length(x))
        x
    end
end
Base.fill{T}(::Type{CuArray}, value::T, dims::NTuple) = fill!(CuArray{T}(dims), value)

Base.reshape{T,N}(x::CuArray{T}, dims::NTuple{N,Int}) = CuArray{T,N}(x.ptr, dims)
Base.reshape{T}(x::CuArray{T}, dims::Int...) = reshape(x, dims)

curand{T}(::Type{T}, dims::NTuple) = CuArray(rand(T,dims))
curand{T}(::Type{T}, dims::Int...) = CuArray(rand(T,dims))
curand(dims::NTuple) = curand(Float64, dims)
curand(dims::Int...) = curand(Float64, dims)

@generated function Base.cat{T,N}(dim::Int, xs::CuArray{T,N}...)
    f = CuFunction("""
    __global__ void f(Array<$T,$N> y, Array<$T,$N> *xs) {
        int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
        int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
        if (idx_x < y.length()) {

            //int subs[$N];
            //getindex(subs);

            //$T x = xs(subs);
            //subs[] = cumdims[idx_y];
            //y(subs) = x;
        }
    }
    """)
    quote
        split = Array{Cint}(length(xs))
        split[1] = 1
        for i = 2:length(xs)
            split[i] = split[i-1] + size(xs[i-1],dim)
        end

        cumdim = split[end] + size(xs[end],dim)
        dims = ntuple(i -> i == dim ? cumdim : size(xs[1],i), ndims(xs[1]))
        y = CuArray{T}(dims)
        $f(y, (y,y), dx=length(dy))

        #x3ds = map(reshape3d, xs)
        #p = Ptr{Void}(map(x -> Ptr{T}(reshape3d(x)), xs))
        #$f(x3d, y3d, dx=size(x3d,1), dy=size(x3d,2), dz=size(x3d,3), bx=1, by=by, bz=1)

        #p = Ptr{Void}[Ptr{T}(xs[i].ptr) for i=1:length(xs)]
        #$f(p, dx=10)
    end
end


text:
# Translating
# https://proplot.readthedocs.io/en/latest/basics.html#Subplot-grids
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

f, axs = proplot.subplots(ncols=4, nrows=4, axwidth=1.2)
[ax.format(
    xlabel="xlabel", ylabel="ylabel", suptitle="Subplot grid demo",
    grid=false, xlim=(0, 50), ylim=(-4, 4)
) for ax in axs]

# Must reshape axs (and transpose it...)
axs = permutedims(reshape(axs, 4, 4), [2, 1])

# Various ways to select subplots in the subplot grid
[ax.format(color="gray7", facecolor="gray3", linewidth=1) for ax in axs[:, 1]]
[ax.format(color="red", facecolor="gray3", linewidth=1) for ax in axs[1, :]]
axs[1].format(color="black", facecolor="gray5", linewidth=1.4)
[ax.format(facecolor="gray1") for ax in axs[2:end, 2:end]]
for ax in axs[2:end, 2:end]
    ax.plot(cumsum(rand(50, 5) .- 0.5, dims=1), cycle="Grays", lw=2)
end
f.savefig("figs/Subplot-grids.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples